gpinit

 init-routine

0001 function p=gpinit(p) 
0002 % init-routine
0003 p=stanparam(p); p.pstyle=3; p.cm=cool;
0004 p.neq=2; p.f=@gpf; p.jac=@gpjac; p.outfu=@gpbra; 
0005 pre=sprintf('%s',inputname(1)); p=setfn(p,pre);
0006 lx=5; ly=1*lx; p.tol=1e-8; [p.geo,p.bc]=recnbc2(lx,ly); 
0007 p.bc=gnbc(2,4,[[100 0];[0 100]],[0;0]); p.bcf=@(p,u,lam) p.bc;
0008 nx=30; ny=nx; p=stanmesh(p,nx,ny); p=setbmesh(p); 
0009 p.dsmin=0.001; p.dsmax=0.5; p.dlammax=0.1; p.nsteps=20; 
0010 p.jsw=0; p.vsw=2; p.amod=0; p.bifchecksw=0; p.spcalcsw=0; % 0/1
0011 p.pcmp=10; p.imax=50; p.bpcmp=3; p.cm='cool'; 
0012 p.lammax=1.5; p.smod=5; p.errbound=0; p.parasw=0; p.smod=5; 
0013 %%%%%%%%%%%%%%%%%%%%%%%  starting point %%%%%%%%%%%%%%%%%%%%%%%%%%%%
0014 p.lam=0;p.ds=0.1;p.pa=1;p.mu=2; % pa & mu are parameters in gpf.
0015 p.xi=1/p.np; p.tau=zeros(p.neq*p.np+1,1); p.tau(p.neq*p.np+1)=1; 
0016 x=p.points(1,:)'; y=p.points(2,:)'; p.r=x.^2+y.^2; 
0017 u=2./cosh(p.r); % monopole (to work with m=1,2 first meshref the monopole)
0018 v=0.0*u; u0=[u v]; p.u=reshape(u0,p.neq*p.np,1); 
0019 res=norm(resi(p,p.u,p.lam), p.normsw); fprintf('initial res=%g\n',res); %return;
0020 p.jsw=1; p.imax=10;
0021 [p.u,res]=nlooppde(p,p.u,p.lam); fprintf('first res=%g\n',res);
0022 plotsol(p,1,1,3); p.tau=zeros(p.neq*p.np+1,1); p.tau(1)=1; 
0023 p=meshref(p,'maxt',2500); 
0024 % now initial quadpole guess on the already refined mesh
0025 x=p.points(1,:)'; y=p.points(2,:)'; p.r=x.^2+y.^2; p.phi=angle(x+1i*y); 
0026 p.tau=zeros(p.neq*p.np+1,1); p.tau(p.neq*p.np+1)=1; 
0027 m=2; n=0; p.lam=0;
0028 am=pi*factorial(m+1)/2; gm=am; bm=pi*factorial(2*m)/2^(2*m+5); 
0029 dm=pi*factorial(m)/2; a2=(sqrt(dm^2*p.mu^2+12*am*gm)+dm*p.mu)/(6*gm); 
0030 a=sqrt(a2)*(n+1); A2=(2*a2*gm-dm*p.mu)/(bm*3); A=sqrt(A2); 
0031 u=A*tfu(p.r/a,m,n).*(cos(m*p.phi))*(n+1); v=0.0*u; u0=[u v]; 
0032 p.u=reshape(u0,p.neq*p.np,1); plotsol(p,1,1,3); 
0033 p.ngen=10;p.adafac=1.5;p.maxt=p.adafac*p.nt; p.xi=1/p.np; 
0034 res=norm(resi(p,p.u,p.lam), p.normsw); fprintf('initial res=%g\n',res);
0035 p.imax=10;[p.u,res]=nlooppde(p,p.u,p.lam); fprintf('first res=%g\n',res);
0036 plotsol(p,1,1,3);p=meshref(p,'maxt',5000); p=setbmesh(p);
0037 
0038