acdemo

 demo script for Allen-Cahn equation with Dirchlet BC over rectangle

0001 % demo script for Allen-Cahn equation with Dirchlet BC over rectangle
0002 close all; clear all; 
0003 mypause('Init and continue trivial branch of Allen-Cahn'); 
0004 p=[];p=acinit(p);p=cont(p);
0005 mypause('Switch branches at bifurcation points and continue'); 
0006 q=swibra('p','bp1','q',0.1);
0007 % 4th argument of swibra is ds (aux argument)
0008 q.amod=5;q.nsteps=10;q.lammax=3;q.smod=3;q=cont(q);
0009 r=swibra('p','bp2','r',0.1);
0010 r.amod=5;r.nsteps=15;r.lammax=4; r=cont(r);
0011 % ----------------------------------------------------
0012 mypause('Plot bifurcation diagram (L2-norm over lam)'); 
0013 figure(3);clf; cmp=2; 
0014 plotbra(p,3,cmp,'ms',8,'lw',6); 
0015 plotbra(q,3,cmp,'ms',8,'lab',10,'lw',4, 'cl', 'b');
0016 plotbra(r,3,cmp,'ms',8,'lab',10,'lw',4, 'cl', 'r');
0017 xlabel('\lambda');ylabel('||u||_2');
0018 mypause('Plot some solutions from file'); 
0019 plotsolf('q','p9',5,1,2);plotsolf('r','p10',7,1,2); 
0020 % ----------------------------------------------------
0021 mypause('switch to contuation in diffusion constant'); 
0022 clf(2); w=q; w=setfn(w); w.up1=w.lam; % save the old lam as 'user-parameter'
0023 w.lam=0.25; % reset the new parameter lamda (=diff. coeff.)
0024 w.lammin=0.1;w.nsteps=20; w.restart=1; w.ds=-0.01; 
0025 w.xi=1e-6; % small xi since problem quite sensitive in diff coeff.
0026 w.f=@acfmu;w.jac=@acjacmu;  % new f, jac since new continuation param,
0027 w.jsw=1; % jsw=1 since Glam by finite differences
0028 %w.f=@acmuf; w.jsw=3; %alternative to above 2 lines; slower on fine meshes
0029 w=cont(w);plotsol(w,8,1,2);title('u at \mu=0.1, \lambda=2');
0030 % ----------------------------------------------------
0031 fprintf('Example for tint: perturb point 3 from q branch into \n');
0032 mypause('unstable direction (in two directions) and run tint'); 
0033 q=loadp('q','p3','q'); 
0034 % calculate unstable direction by finding leading Eval/Evec of Jacobian
0035 Gu=getGu(q,q.u,q.lam,r); 
0036 [K,MId]=assema(q.points,q.tria,0,1,0); 
0037 q.evopts=[]; q.evopts.disp=0;[phi1,mu1]=eigs(Gu,MId,1,0,q.evopts);
0038 mypause('Perturb in one direction');
0039 qa=q;qa.u=q.u-0.2*phi1;qa=tint(qa,0.2,200,40);plotsol(qa,9,1,1);
0040 mypause('Perturb in the other direction');
0041 qb=q;qb.u=q.u+0.2*phi1;qb=tint(qb,0.3,600,60);