bratudemo

 demo-script for Bratu. Note that the basic bif-diagram and 
 solution plots only require the first 30 lines

0001 % demo-script for Bratu. Note that the basic bif-diagram and
0002 % solution plots only require the first 30 lines
0003 close all; clear all; 
0004 p=[];
0005 fprintf('Init p for Bratu´s problem\n'); 
0006 p=bratuinit(p); % init
0007 mypause('Ready to call p=cont(p)'); 
0008 p=cont(p); % homogeneous branch
0009 fprintf('Ready to call q=swibra(´p´,´bp1´,´q´)\n'); 
0010 mypause('(branch bifurcating from homogeneous branch)'); 
0011 q=swibra('p','bp1','q'); q.lammin=0.1;q.nsteps=20;
0012 fprintf('See figure(6) for tangent of bifurcating branch.\n');
0013 mypause('Ready to call q=cont(q)'); 
0014 q=cont(q); 
0015 fprintf('Ready to call r=swibra(´q´,´bp1´,´r´) and r=cont(r)\n'); 
0016 mypause('(branch bifurcating from q-branch)'); 
0017 r=swibra('q','bp1','r'); r.lammin=0.05;r.nsteps=20;r=cont(r); 
0018 % ----------- postprocessing, plot BifDiagram
0019 mypause('Ready to plot BifDiagram to figure(3)'); 
0020 figure(3);clf(3);
0021 plotbra(p,3,2,'ms',12,'lwst',5,'lwun',2,'fs',16,'cl','k');
0022 plotbra(q,3,2,'ms',12,'lwst',6,'lwun',2,'lab',[20],'fs',16,'cl','b');
0023 plotbra(r,3,2,'ms',12,'lwst',6,'lwun',2,'lab',[20],'fs',16,'cl','r');
0024 axis([0.05 0.37 0.5 3.5]);
0025 % ----------- plot solutions
0026 mypause('Ready to plot some solutions'); 
0027 plotsolf('q','p20',7,1,1); % mesh plot
0028 figure(7); view(-30,50); xlabel('x');ylabel('y');
0029 plotsolf('r','p20',8,1,3); % rendered 3D plot
0030 figure(8); view(-30,50); xlabel('x');ylabel('y');
0031 % ----------- meshcheck
0032 mypause('Assess mesh at end of q-branch via mesh refinement');
0033 q0=loadp('q','p20','q0');   % load point to perform various checks
0034 [q0r,ud]=meshcheck(q0,1);   % first 'ad hoc meshcheck': refine mesh
0035 fprintf('See fig(6) for difference between old and new solution,\n');
0036 fprintf('and fig(1) for the new solution.\n');
0037 % ----------- mesh refinement
0038 mypause('An example for refining mesh aiming at error-bound 0.02'); 
0039 q0r=meshref(q0,'eb',0.02);  % refine aiming at error-bound
0040 plotsol(q0r,9,1,1);view(-30,50); xlabel('x');ylabel('y');
0041 fprintf('See fig(6) for new mesh and fig(9) for new solution\n'); 
0042 % ----------- errcheck
0043 fprintf('An error-estimate is calculated via err=errcheck(q0);\n'); 
0044 err=errcheck(q0r); fprintf('error-estimate %g\n',err); 
0045 % ----------- test of different adaption strategies
0046 fprintf('Test different adaption strategies; first load p10 on q branch,\n');
0047 fprintf('then run without adaption, with adaption every 10th step, and with\n'); 
0048 mypause('adaption if error-estimate > error-bound=0.1'); 
0049 q0=loadp('q','p10','q0');q0.lammin=0.05;q0=resetc(q0); q0.nsteps=30; 
0050 q0.neig=100; % need many eigenvalues since the negative one becomes quite large
0051 q0.errbound=0.1; q0b=q0; % save the point for restarts
0052 q0=cont(q0); % no refinement
0053 q1=q0b;q1.amod=10; q1=cont(q1);     % adapt every 10 steps
0054 q2=q0b;q2.errchecksw=2;q2=cont(q2); % adapt when err-est>errbound
0055 mypause('plot the error-estimates for the different adaption strategies to fig(5)'); 
0056 % ---------- error plot for different refinements
0057 figure(5);clf;hold on;
0058 ploterr(q0,5,3);ploterr(q1,5,3);ploterr(q2,5,3); legend('q0','q1','q2');
0059 mypause('plot the branches for the different adaption strategies to fig(3)'); 
0060 % ---------- branch plots for different refinements
0061 clf(3);cmp=2; 
0062 plotbra(q0,3,cmp,'labi',10,'cl','k','ms',8);
0063 plotbra(q1,3,cmp,'labi',10,'cl','b','ms',8);
0064 plotbra(q2,3,cmp,'labi',10,'cl','r','ms',8);
0065 legend('q0','q1','q2');
0066 % ---------- check how many triangles needed to fulfill error-estimate
0067 mypause('Finally, check how much refinement is needed for error-est < 0.0025:'); 
0068 q0r=meshref(q0,'eb',0.0025,'maxt',10^5, 'ngen',20);