rbconvcmds

 command templates for Rayleigh-Benard; run cell-by-cell; 
 note that these commands already incorporate some knowledge 
 like approximate location of Bif. points and occurrence 
 of imperfect pitchfork, which has been found before (by toying 
 around with pde2path)

0001 % command templates for Rayleigh-Benard; run cell-by-cell;
0002 % note that these commands already incorporate some knowledge
0003 % like approximate location of Bif. points and occurrence
0004 % of imperfect pitchfork, which has been found before (by toying
0005 % around with pde2path)
0006 close all; clear all; 
0007 %% some steps to get first Bif-point
0008 ns=[]; ns=rbconvinit(ns); ns=cont(ns);
0009 % for speedup, 'jump' to near next Bif-point
0010 ns.lam=771; ns.ds=0.1; ns=cont(ns); 
0011 %% bifurcating branches, use pmcont for speedup
0012 ns1=swibra('ns','bp1','ns1');ns1.ds=2;ns1=pmcont(ns1);
0013 ns2=swibra('ns','bp2','ns2');ns2.ds=2;ns2=pmcont(ns2);
0014 %% Plot Bif. diagram
0015 figure(3);clf;cmp=1;
0016 plotbra(ns,3,cmp,'lw',6,'ms',10,'cl','b');
0017 plotbra(ns1,3,cmp,'ms',10,'cl','r','lab',10);
0018 plotbra(ns2,3,cmp,'ms',10,'cl','k','lab',8);
0019 axis([739 940 0 3]); xlabel('R'); ylabel('max\psi');
0020 %% plot solutions
0021 p1=loadp('ns1','p10','p1'); 
0022 p2=loadp('ns2','p8','p1'); 
0023 arrowplot(p1,8,1); arrowplot(p2,9,1); 
0024 %%  stress-free boundary conditions, trivial branch
0025 sf=[];sf=rbconvinit_stressfree(sf);
0026 sf.lam=665.5;sf=cont(sf); sf.lam=761;sf=cont(sf); 
0027 %% first bifurcating branch,
0028 sf1=swibra('sf','bp1','sf1');sf1.ds=2;df1.nsteps=10; sf1=pmcont(sf1); 
0029 %% second bifurcating branch
0030 sf2=swibra('sf','bp2','sf2');
0031 sf2.ds=2;sf2.lammax=845;sf2=pmcont(sf2); % stop before pitchfork
0032 sf2.ds=1;sf2.dsmax=2;sf2.lammax=900;sf2.nsteps=50;sf2=cont(sf2);
0033 %% Switch to disconnected branch by time-integration:
0034 p=sf2;p=resetc(p); r=resi(p,p.u,p.lam); Gu=getGu(p,p.u,p.lam,r); 
0035 [K,MId]=assema(p.points,p.tria,0,1,zeros(p.neq,1)); 
0036 [phi1,mu1]=eigs(Gu,MId,1,0,p.evopts);
0037 p1=p;p1.u=p.u+0.01*phi1/norm(phi1);
0038 p1=tint(p1,0.1,400,50); 
0039 %% Now continue from p1
0040 sf3=p1; sf3.restart=1; sf3.ds=-sf3.ds; sf3.lammax=sf3.lammax+2; 
0041 sf3.bifchecksw=0;sf3.dlammax=10;sf3.nsteps=100;sf3=cont(sf3);
0042 %% Plot Bif Diagram
0043 clf(3);cmp=1;
0044 plotbra(sf,3,cmp,'lw',5,'ms',10,'cl','b');
0045 plotbra(sf1,3,cmp,'ms',10,'cl','r','lab',40,'inc',40); 
0046 plotbra(sf2,3,cmp,'ms',10,'cl','k','lab',[30 60]);
0047 plotbra(sf3,3,cmp,'ms',10,'cl','b','lab',65,'inc',65);
0048 axis([660 920 0 3]); xlabel('R'); ylabel('max\psi');
0049 %% Plot solutions
0050 p1=loadp('sf1','p40','p1'); p2=loadp('sf2','p20','p2'); 
0051 p3=loadp('sf2','p60','p3'); p4=loadp('sf3','p65','p4'); 
0052 arrowplot(p1,8,1); arrowplot(p2,9,1); 
0053 arrowplot(p3,10,1); arrowplot(p4,11,1); 
0054 %% other useful stuff
0055 sf=loadp('sf','p12','sf');
0056 sf1=loadp('sf1','p20','sf1');
0057 sf2=loadp('sf2','p71','sf2');
0058 sf3=loadp('sf3','p85','sf3');
0059