JBDMV-matlab/seco/brsG.m 0000755 0001750 0001750 00000000415 14041401156 012773 0 ustar hu hu function r=brsG(p,u)
%f=nodalf(p,u);
par=u(p.nu+1:end); j0=par(1); al=par(2); tau=par(3); D=par(4);
u1=u(1:p.np); u2=u(p.np+1:2*p.np);
f1=(u2-u1)./((u2-u1).^2+1)-tau*u1; f2=al*(j0-(u2-u1));
f=[f1;f2];
Ks=p.mat.K; K=[Ks 0*Ks; 0*Ks D*Ks];
r=K*u(1:p.nu)-p.mat.M*f; JBDMV-matlab/seco/sGjac.m 0000755 0001750 0001750 00000000563 14117707476 013153 0 ustar hu hu function Gu=sGjac(p,u) % Jac for MWBS'97 semiconductor model
n=p.np; [f1u,f1v,f2u,f2v]=nodaljac(p,u); % Jac of 'nonlinearity'
Fu=[[spdiags(f1u,0,n,n),spdiags(f1v,0,n,n)];% put f_u into block matrix
[spdiags(f2u,0,n,n),spdiags(f2v,0,n,n)]];
D=u(p.nu+4); Ks=p.mat.K; K=[Ks 0*Ks; 0*Ks D*Ks];% diffusion matrix
Gu=K-p.mat.M*Fu; % Jacobian JBDMV-matlab/seco/secoinit.m 0000755 0001750 0001750 00000001725 14117707321 013727 0 ustar hu hu function p=secoinit(lx,nx,par) % init for MWBS'97 semiconductor model
p=stanparam(); p.nc.neq=2; % init with stanparam, 2-compo-system
p.sol.ds=-0.01; p.sol.dsmax=0.05; p.sw.bifcheck=2; % reset some pars
p.fuha.outfu=@secobra; p.plot.bpcmp=8; % output function handle, and compo
p=setbel(p,0,1e-6,10,@lss); % use bordered elim. lss (always good in 1D),
% with border width 0, tolerance 1e-6, and at most 10 iterations
pde=stanpdeo1D(lx,2*lx/nx); p.vol=2*lx; % standard 1D PDE object
p.np=pde.grid.nPoints; p.pdeo=pde; % store number of grid-points, and pdeo
p.nu=p.np*p.nc.neq; p.sol.xi=1/p.nu; % DoFs, and weight for arclength
p.nc.lammin=0; p.nc.lammax=5; p.nc.dsmax=0.1; % lam-range, and max stepsize
j0=par(1); tau=par(3); als=j0/(tau*(j0^2+1)); us=als+j0; % initial sol
u=als*ones(p.np,1);v=us*ones(p.np,1);p.u=[u;v;par]; % append pars and store
p.nc.ilam=1; % select the (initial) primary active par, here j0
p.sw.sfem=-1; p=oosetfemops(p); % set FEM matrices JBDMV-matlab/seco/sG.m 0000755 0001750 0001750 00000000701 14121235065 012450 0 ustar hu hu function r=sG(p,u) % rhs for MWBS'97 semiconductor model
% split u into pars and PDE fields u1 and u2:
par=u(p.nu+1:end); j0=par(1); al=par(2); tau=par(3); D=par(4);
u1=u(1:p.np); u2=u(p.np+1:2*p.np);
% compute nodal 'nonlinearity' (terms without spat.derivatives):
f1=(u2-u1)./((u2-u1).^2+1)-tau*u1; f2=al*(j0-(u2-u1)); f=[f1;f2];
Ks=p.mat.K; K=[Ks 0*Ks; 0*Ks D*Ks]; % the diffusion matrix
r=K*u(1:p.nu)-p.mat.M*f; % the residual JBDMV-matlab/seco/cmds2.m 0000664 0001750 0001750 00000004165 14131115275 013117 0 ustar hu hu %% spatially homogeneous primary Hopf; use tl=40 points in t;
ds=0.1; aux.tl=40; p=hoswibra('0','hpt1',ds,4,'H1',aux); p.nc.dsmax=2;
p.sw.bifcheck=0; p.file.smod=1;% switch off bif.detection, save every point
p.hopf.flcheck=1; p.hopf.fltol=1e-6; % switch on multiplier comp
p.sw.verb=0; p=setbel(p,2,1e-4,5,@lss); p=cont(p,40); % use bel, and go!
%% splice together Hopf and Turing: choose xcut and j0 by trial and error
p1=loadp('T1','pt25'); p=loadp('H1','pt25'); xcut=-160; j0guess=3.35;
x=getpte(p);idx=find(x>xcut);%find indizes where to replace Hopf by Turing
p.u(p.nu+1)=j0guess; p.hopf.lam=j0guess; % set j0 to the guess j0guess
for j=1:p.hopf.tl % put Turing-soln into Hopf-data
p.hopf.y(idx,j)=p1.u(idx); p.hopf.y(idx+p.np,j)=p1.u(idx+p.np);
end
for i=1:1; hoplot(p,1,i); end; pause % plot for checking
p.sw.verb=2; p.sol.ds=-0.01; p.nc.dsmin=0.01;
p.nc.imax=20; p.nc.tol=1; p.hopf.flcheck=0; p.nc.dsmax=5;
p=resetc(p); p=setfn(p,'lh1a'); p.hopf.nfloq=10;
p=cont(p,1); % 1 step to see if initial guess was OK
%% a few more steps to get on the branch
p=resetc(p); p.nc.imax=10; p.nc.tol=1e-2; p=cont(p,3); p.nc.tol=1e-4; p=cont(p,3);
%% continue !
p=resetc(p); huclean(p);
p.nc.imax=10; p.nc.tol=1e-6; p.hopf.flcheck=0; p=cont(p,1); p.file.smod=10; p=cont(p,50);
%% other direction
p=loadp('lh1a','pt1','lh1b'); p.sol.ds=-p.sol.ds; p=resetc(p); p=cont(p,200);
%% BD
fnr=3; figure(fnr); clf; cmp=8; % desired branch-compo,
ptli0=[16 22]; ptli1=[1 40]; ptli2=[20 90 110 170 230 290]; % label lists
plotbra('lh1a','pt40',fnr,cmp,'cl','r','lab',ptli1,'lp',40);
plotbra('lh1b','pt300',fnr,cmp,'cl','m','lab',ptli2,'fp',4);
plotbra('H1','pt50',fnr,cmp,'cl','b','lab',ptli0,'fp',16);
xlabel('j0'); ylabel('||U||');
%% soln plots, with Fl-multipl., uncomment for desired directory
v=[15 60]; aux.nfloq=30; % view, and number of Floquet-multipliers
dir='H1';ptli=ptli0;
%dir='lh1a';ptli=ptli1;
%dir='lh1b';ptli=ptli2;
for i=ptli;
pt=['pt' mat2str(i)]; hoplotf(dir,pt,1,1); xlabel('t');
figure(1); title([dir '/' pt]); ylabel('t'); zticks([4 8]);
view(v); shading interp; muv1=floqap(dir,pt,aux);
pause
end JBDMV-matlab/seco/cmds1.m 0000664 0001750 0001750 00000006507 14131113356 013116 0 ustar hu hu %% MWBS97 model, trivial and Turing branches, and bif from Turing
al=0.02; j0=3.3; tau=0.05; D=8; kc=(al*tau/D)^0.25; nw=8; % parameters,
lx=nw*pi/kc; nx=nw*40; par=[j0;al;tau;D]; dir='0'; % dom.size and nx
p=secoinit(lx,nx,par); p=setfn(p,dir); % init, set output dir for tr.branch
%% compute guess for Hopf-spectral shift, then cont trivial branch
p=initeig(p,0.5); p.nc.neig=[10,10]; p.nc.eigref(1)=-0.05;
p.nc.dsmax=0.05; p.sw.verb=1; p.file.smod=5; p=cont(p,30);
%% primary Turing mode
p=swibra('0','bpt1','T1'); p.nc.dsmax=0.05; p.sol.ds=-0.02; p=cont(p,60);
%% T-mode via cswibra, useful for eigenvector inspection
aux.m=6; aux.besw=0; p0=cswibra('0','bpt2',aux);
%% now generate bif directions and go
p=gentau(p0,[1],'T2'); p.nc.dsmax=0.1; p=cont(p,40);
p=gentau(p0,[0 0 1],'T3'); p.nc.dsmax=0.1; p=cont(p,40);
%% T4-mode via cswibra
aux.m=6; aux.besw=0; p0=cswibra('0','bpt3',aux); pause
p=gentau(p0,-1,'T4'); p.nc.dsmax=0.1; p=cont(p,40);
%% T-front snake
clf(1); p=swibra('T1','bpt1','T1-1'); p.nc.dsmax=0.05; p.sw.foldcheck=1;
p.sol.ds=-0.02; p.file.smod=10; p.nc.mu2=0.01; p=cont(p,300);
%% HP10 on front-snake
para=4; ds=0.1; aux.tl=30; aux.dlam=0;
p=hoswibra('T1-1','hpt10',ds,para,'T1-1-H10',aux); pause
p.hopf.flcheck=0; p.sw.verb=0;
p=setbel(p,2,1e-4,5,@lss); p.nc.dsmax=0.5; p=cont(p,50);
%% HP18
p=hoswibra('T1-1','hpt18',ds,para,'T1-1-H18',aux); pause
p.hopf.flcheck=0; p.sw.verb=0;
p=setbel(p,2,1e-4,5,@lss); p.nc.dsmax=0.5; p=cont(p,50);
%% loc.Turing snake
clf(1); p=swibra('T1','bpt2','T1-2'); p.nc.dsmax=0.05; p.sw.foldcheck=1;
p.file.smod=10; p.sol.ds=-0.02; p=cont(p,200);
%% HP8 on loc.Turing snake
para=4; ds=0.1; dsmax=1.1; xi=1e-3; aux.tl=20;
p=hoswibra('T1-2','hpt8',ds,para,'T1-2-H8',aux); pause
p.hopf.flcheck=0; p.sw.verb=0; p=setbel(p,2,1e-4,5,@lss); p.nc.dsmax=0.5; p.plot.bpcmp=8;
p.file.smod=10; p=cont(p,100);
%% BD-plot,
fnr=3; figure(fnr); clf; cmp=8; plotbra('0',fnr,cmp,'cl','k','lsw',0);
plotbra('T1','pt60',fnr,cmp,'cl','b','lab',20);
plotbra('T2','pt40',fnr,cmp,'cl',p2pc('b2'),'lsw',0);
plotbra('T3','pt40',fnr,cmp,'cl',p2pc('g1'),'lsw',0);
plotbra('T4','pt40',fnr,cmp,'cl',p2pc('g2'),'lab',25);
plotbra('T1-1','pt400',fnr,cmp,'cl',p2pc('o1'),'lab',350,'lsw',0,'lp',380,'hplab',[10 18]);
plotbra('T1-1-H10','pt100',fnr,cmp,'cl',p2pc('r1'),'lsw',0,'lab',[20 50]);
plotbra('T1-1-H18','pt100',fnr,cmp,'cl',p2pc('r2'),'lsw',0);
axis([3.05 3.7 0 1.8]); xlabel('j_0'); ylabel('||u_1||_*');
%% BD plot, zoom of loc T-snake
figure(fnr);clf;
plotbra('T1-2','pt100',fnr,cmp,'cl',p2pc('o2'),'lsw',0,'hplab',8,'lp',175);
plotbra('T1-2-H8','pt100',fnr,cmp,'cl',p2pc('r2'),'lsw',0,'lab',[20 80],'lp',90);
axis([3.24 3.35 1.15 1.45]); xlabel('j_0'); ylabel('||u_1||_*');
%% soln plots, steady
figure(1); clf;
plotsol('T1','pt30',1,1,1); noltix; pause;
plotsol('T4','pt35',1,1,1); noltix; pause;
plotsol('T1-1','hpt10',1,1,1); noltix; pause;
plotsol('T1-1','hpt14',1,1,1); noltix; pause;
plotsol('T1-1','hpt18',1,1,1); noltix; pause;
plotsol('T1-1','pt370',1,1,1); noltix
%% hoplots, with Fl-multipl., uncomment the desired directory
v=[15 60]; aux.nfloq=20;
dir='T1-1-H10'; ptli=[20 50];
dir='T1-2-H8'; ptli=[20 80];
for i=ptli; %120 % 10:10:130;
pt=['pt' mat2str(i)]; hoplotf(dir,pt,1,1); figure(1);
title([dir '/' pt]); view(v); shading interp; ylabel('t');
muv1=floqap(dir,pt,aux); pause
end
JBDMV-matlab/seco/nodaljac.m 0000664 0001750 0001750 00000000434 14117710225 013655 0 ustar hu hu function [f1u,f1v,f2u,f2v]=nodaljac(p,u) % for seco
n=p.np; u1=u(1:n); u2=u(n+1:2*n); den=(u2-u1).^2+1; % u, and a denominator
par=u(p.nu+1:end); al=par(2); tau=par(3); % parameters in nonlinearity
f1u=-1./den+2*(u2-u1).^2./(den.^2)-tau; f1v=-f1u-tau;
f2u=al*ones(n,1); f2v=-f2u; JBDMV-matlab/seco/hobra.m 0000664 0001750 0001750 00000001604 14041476663 013210 0 ustar hu hu function out=hobra(p,u)
% hobra: output for bifurcation diagram, setup for Hopf
% [paras, T,max(u),min(u),||u||]
switch p.sw.spcont
case 0; np=p.nu/p.nc.neq; n=p.nu; % use nu, i.e. all compos of u
case 1; np=p.nu/p.nc.neq; n=p.nu/2;
case 2; np=p.nu/p.nc.neq; n=p.nu/2;
case 3; np=(p.nu-1)/p.nc.neq; n=(p.nu-1)/3;
end
par=u(p.nu+1:end); % par', pause
[as,us]=getss(par);
if p.sw.para>2 % Hopf setup; get data from p.hopf
ho=p.hopf; m1=ho.T;
m2=max(max(ho.y(1:np,:)-as)); m3=min(min(ho.y(1:np,:)-as));
l2v=zeros(1,ho.tl);
for i=1:length(ho.t); l2v(i)=(ho.y(1:np,i)'-as)*(ho.tom.M(1:np,1:np)*(ho.y(1:np,i)-as)); end
l2=trapz(ho.t,l2v); m4=sqrt(l2/p.vol);
else % steady state continuation
m1=0; m2=max(u(1:np)-as); m3=min(u(1:np)-as);
l2=(u(1:np)-as)'*(p.mat.M(1:np,1:np)*(u(1:np)-as)); m4=sqrt(l2/p.vol);
end
out=[par; m1; m2; m3; m4];
JBDMV-matlab/seco/getss.m 0000664 0001750 0001750 00000000117 14035603760 013231 0 ustar hu hu function [as,us]=getss(par)
j0=par(1); T=par(3); as=j0/(T*(j0^2+1)); us=as+j0; JBDMV-matlab/seco/oosetfemops.m 0000644 0001750 0001750 00000000301 14121235007 014427 0 ustar hu hu function p=oosetfemops(p)
% generate FEM matrices, here just mass M and Neumann-Laplacian K
[K,M,~]=p.pdeo.fem.assema(p.pdeo.grid,1,1,1);
p.mat.M=[M 0*M;0*M M]; p.mat.K=K; % store matrices JBDMV-matlab/seco/spjac.m 0000664 0001750 0001750 00000000722 14041335064 013202 0 ustar hu hu function Guuph=spjac(p,u) % for BP and FP cont
n=p.np; u1=u(1:p.np); u2=u(p.np+1:2*p.np);
d1=(u2-u1).^2+1; d2=d1.^2; d3=d1.*d2;
f1uu=-6*(u2-u1)./d2+8*(u2-u1).^3./d3; f1uv=-f1uu; f1vv=f1uu;
f2uu=0*u1; f2uv=f2uu; f2vv=f2uu;
ph1=u(p.nu+1:p.nu+n); ph2=u(p.nu+n+1:2*p.nu);
M1=spdiags(f1uu.*ph1+f1uv.*ph2,0,n,n);
M2=spdiags(f1uv.*ph1+f1vv.*ph2,0,n,n);
M3=spdiags(f2uu.*ph1+f2uv.*ph2,0,n,n);
M4=spdiags(f2uv.*ph1+f2vv.*ph2,0,n,n);
Guuph=-p.mat.M*[[M1 M2]; [M3 M4]]; JBDMV-matlab/seco/spufu.m 0000644 0001750 0001750 00000004266 14042633632 013254 0 ustar hu hu %%
% spufu, adadption of STANUFU to plot spectrum of lin. around hom sol
function [p,cstop]=spufu(p,brout,ds)
p=ulamcheck(p); % test if a desired lambda value has been passed
brplot=brout(length(bradat(p))+p.plot.bpcmp); %y-axis value in bif-figure
fprintf('%4i %s %s %5.2e %4i %s %s ', ...
p.file.count,printcon(getlam(p)),printcon(brplot),p.sol.res,p.sol.iter, ...
p.sol.meth,printcon(ds));
if(p.sw.errcheck>0) fprintf('%5.2e ',p.sol.err);end
if(p.sw.spcalc==1) fprintf(' %2i ', p.sol.ineg); end;
npr=length(p.sw.bprint);
for i=1:npr; fprintf('%s ',printcon(brout(length(bradat(p))+p.sw.bprint(i)))); end;
fprintf('\n');
cstop=0;
if(getlam(p)
p.nc.lammax)
fprintf(' lam=%g > lammax=%g, stopping\n',getlam(p),p.nc.lammax); cstop=1;
end
% addition to STANUFU: plot the dispersion relation!
n=p.np;nu=p.nu; par=p.u(p.nu+1:end);u=[p.u(1);p.u(n+1)]; % hom-state 2-vector
u=[u;par];p.np=1;p.nu=2;[f1u,f1v,f2u,f2v]=nodaljac(p,u); J=[[f1u f1v];[f2u f2v]];
d=par(4); % diffusion param.; this and k-range usually only problem dep. things
kmax=0.2; kv=0:0.002:kmax; kl=length(kv); muv=zeros(2,kl); % provide wave-nrs and mem
x=getpte(p); lx=max(x); dk=pi/2/lx;
kv2=0:dk:kmax; kl2=length(kv2); muv2=zeros(2,kl2);
for i=1:kl % now loop over wave-nr and compute Evals
k=kv(i); K=[[k^2 0];[0 d*k^2]]; A=J-K; % Jac in Fourier space
mu=eig(A); [mus, ix]=sort(real(mu),'descend'); % sorted eigenvalues
for j=1:2; muv(j,i)=mu(ix(j)); end
end
for i=1:kl2 % now loop over wave-nr and compute Evals
k=kv2(i); K=[[k^2 0];[0 d*k^2]]; A=J-K; % Jac in Fourier space
mu=eig(A); [mus, ix]=sort(real(mu),'descend'); % sorted eigenvalues
for j=1:2; muv2(j,i)=mu(ix(j)); end
end
figure(10); clf; plot(kv,real(muv(1,:)),kv,imag(muv(1,:)),kv,0*kv,'--k'); % plot leading Eval
title(['par=' mat2str(par,3)],'fontsize',p.plot.fs); set(gca,'FontSize',p.plot.fs);
p.np=n; p.nu=nu; % reset # spatial points to full problem
hold on
plot(kv2,real(muv2(1,:)),'*b',kv2,imag(muv2(1,:)),'*r');
legend('real','imag'); xlabel('k'); title(''); axis tight;
JBDMV-matlab/seco/bpjac.m 0000664 0001750 0001750 00000000752 14041322370 013160 0 ustar hu hu function Guuph=bpjac(p,u) % for BP and FP cont
n=p.np; u1=u(1:p.np); u2=u(p.np+1:2*p.np);
ov=ones(n,1); den=(u2-u1).^2+1; d2=den.^2; d3=d2.*den;
f1uu=-6*(u2-u1)./d2+8*(u2-u1).^3./d3;
f1uv=-f1uu; f1vv=f1uu;
f2uu=0*u1; f2uv=f1uv; f2vv=f1vv;
ph1=u(p.nu+1:p.nu+p.np); ph2=u(p.nu+p.np+1:2*p.nu);
M1=spdiags(f1uu.*ph1+f1uv.*ph2,0,n,n);
M2=spdiags(f1uv.*ph1+f1vv.*ph2,0,n,n);
M3=spdiags(f2uu.*ph1+f2uv.*ph2,0,n,n);
M4=spdiags(f2uv.*ph1+f2vv.*ph2,0,n,n);
Guuph=-p.mat.M*[[M1 M2]; [M3 M4]]; JBDMV-matlab/seco/hpjac.m 0000664 0001750 0001750 00000001760 14041340005 013161 0 ustar hu hu function Guuph=hpjac(p,u) % for HP continuation
nu=p.nu; np=p.np; M=p.mat.M(1:np,1:np);
u1=u(1:np); u2=u(np+1:2*np);
ph1r=u(nu+1:nu+np); ph2r=u(nu+np+1:nu+2*np);
ph1i=u(2*nu+1:2*nu+np); ph2i=u(2*nu+np+1:2*nu+2*np);
Guuph=sparse(2*nu,nu);
d1=(u2-u1).^2+1; d2=d1.^2; d3=d1.*d2;
f1uu=-6*(u2-u1)./d2+8*(u2-u1).^3./d3; f1uv=-f1uu; f1vv=f1uu;
f2uu=0*u1; f2uv=f2uu; f2vv=f2uu;
M11=spdiags(f1uu.*ph1r+f1uv.*ph2r,0,np,np);
M12=spdiags(f1uv.*ph1r+f1vv.*ph2r,0,np,np);
M21=spdiags(f2uu.*ph1r+f2uv.*ph2r,0,np,np);
M22=spdiags(f2uv.*ph1r+f2vv.*ph2r,0,np,np);
Guuph(1:np,1:np)=-M*M11; Guuph(1:np,np+1:2*np)=-M*M12;
Guuph(np+1:2*np,1:np)=-M*M21; Guuph(np+1:2*np,np+1:2*np)=-M*M22;
M11=spdiags(f1uu.*ph1i+f1uv.*ph2i,0,np,np);
M12=spdiags(f1uv.*ph1i+f1vv.*ph2i,0,np,np);
M21=spdiags(f2uu.*ph1i+f2uv.*ph2i,0,np,np);
M22=spdiags(f2uv.*ph1i+f2vv.*ph2i,0,np,np);
Guuph(nu+1:nu+np,1:np)=-M*M11; Guuph(nu+1:nu+np,np+1:2*np)=-M*M12;
Guuph(nu+np+1:nu+2*np,1:np)=-M*M21; Guuph(nu+np+1:nu+2*np,np+1:2*np)=-M*M22; JBDMV-matlab/seco/secobra.m 0000664 0001750 0001750 00000001516 14041477101 013521 0 ustar hu hu function out=secobra(p,u) % mod of hobra.m to compute deviation from tr. state
switch p.sw.spcont
case 0; np=p.nu/p.nc.neq; n=p.nu; % use nu, i.e. all compos of u
case 1; np=p.nu/p.nc.neq; n=p.nu/2;
case 2; np=p.nu/p.nc.neq; n=p.nu/2;
case 3; np=(p.nu-1)/p.nc.neq; n=(p.nu-1)/3;
end
par=u(p.nu+1:end); [as,us]=getss(par);
if p.sw.para>2 % Hopf setup; get data from p.hopf
ho=p.hopf; m1=ho.T;
m2=max(max(ho.y(1:np,:)-as)); m3=min(min(ho.y(1:np,:)-as));
l2v=zeros(1,ho.tl);
for i=1:length(ho.t); l2v(i)=(ho.y(1:np,i)'-as)*(ho.tom.M(1:np,1:np)*(ho.y(1:np,i)-as)); end
l2=trapz(ho.t,l2v); m4=sqrt(l2/p.vol);
else % steady state continuation
m1=0; m2=max(u(1:np)-as); m3=min(u(1:np)-as);
l2=(u(1:np)-as)'*(p.mat.M(1:np,1:np)*(u(1:np)-as)); m4=sqrt(l2/p.vol);
end
out=[par; m1; m2; m3; m4]; JBDMV-matlab/seco/cmds1b.m 0000664 0001750 0001750 00000003307 14131115125 013247 0 ustar hu hu %% BP and HP cont, needs a larger delta for FDs for the param-derivatives
p=bpcontini('0','bpt1',2,'bpc1'); p.sol.ds=0.005; p.nc.dsmax=0.2; huclean(p);
p.nc.del=1e-3; p.plot.bpcmp=1; p.nc.tol=1e-6; p=cont(p,100);
%% reverse direction
p=loadp('bpc1','pt5','bpc1b'); p.sol.ds=-p.sol.ds; p=cont(p,20);
%% primary Hopf
p=hpcontini('0','hpt1',2,'hpc1'); p.sol.ds=0.05; p.nc.dsmax=0.2;
p.nc.del=1e-2; p.fuha.spjac=@hpjac; p.sw.spjac=1; %[Ga,Gn]=hpjaccheck(p);
p.plot.bpcmp=1; p.sol.ds=0.01; p=cont(p,50);
p=loadp('hpc1','pt5','hpc1b'); p.sol.ds=-p.sol.ds; p=cont(p,20);
%% standard branch plot
mclf(4); plotbra('bpc1',4,1,'cl','b'); plotbra('hpc1',4,1,'cl','r');
axis([0 0.1 1 3.5]); ylabel('j0'); xlabel('\alpha');
%% plotbradat
mclf(3); p=loadp('bpc1'); h1=plotbradat(p,3,7,8,'-b'); hold on;
set(h1,'linewidth',2); p=loadp('hpc1'); h3=plotbradat(p,3,7,8,'d-r');
p=loadp('bpc1b'); h1=plotbradat(p,3,7,8,'-b'); set(h1,'linewidth',2);
p=loadp('hpc1b'); h3=plotbradat(p,3,7,8,'d-r');
axis([1 3.6 0 0.1]); legend('T_1','H_1'); xlabel('j0'); ylabel('\alpha');
annotation('textarrow',[0.64 0.58],[0.62 0.62],'String','(c)','Fontsize',14);
annotation('textarrow',[0.84 0.84],[0.4 0.28],'String','(d)','Fontsize',14);
text(3.22,0.05,'(b)','Fontsize',14); set(gca,'fontsize',14);
%% MWBS, disp.relations, uncomment the desired par-values in the next line
al=0.04; j0=3.3; %al=0.06; j0=2.4; %al=0.01; j0=3.5;
tau=0.05; D=8; kc=(al*tau/D)^0.25; nw=8; % parameters,
lx=nw*pi/kc; nx=nw*40; par=[j0;al;tau;D]; dir='0b'; % dom.size and nx
p=secoinit(lx,nx,par); p=setfn(p,dir); % init, set output dir for tr.branch
p.fuha.ufu=@spufu; p.sol.ds=1e-6; p=cont(p,1); figure(10);
title(['(j0,\alpha)=(' mat2str(j0) ',' mat2str(al) ')']); JBDMV-matlab/acdc/mypsol.m 0000644 0001750 0001750 00000000411 14131116266 013360 0 ustar hu hu function mypsol(dir,pt);
p=loadp(dir,pt); y2=max(max(p.u(1:p.nu)),1);
plotsol(p); axis([0 1 0 y2]); hold on;
m=round(p.nu/2); % p.u(m-2:m+4)', min(p.u(1:p.nu)) % uncomment to see minima
idx=find(p.u(1:p.nu)<1e-8); x=getpte(p);
plot(x(idx),p.u(idx),'*m'); nola
JBDMV-matlab/acdc/nloopext.m 0000644 0001750 0001750 00000005074 14044240520 013711 0 ustar hu hu function [u,res,iter,Gu,Glam,p]=nloopext(p,u1,ds,varargin)
% NLOOPEXT: newton loop for arclength extended system.
%
% [u,res,iter,Gu,Glam]=nloopext(p,u1,ds) : use p.u and u1 for initial arclength step
% [u,res,iter,Gu,Glam]=nloopext(p,u1,ds,u) : use u and u1 for initial arclength step
%
% * u1 is initial guess, ds initial stepsize usually taken from p-structure.
% * Returns u, res=residuum, iter=number of iterations, Gu,Glam=derivatives.
% * Note: returns full vector u (including auxiliary variables)!
% * Important settings: p.fsol.fsol, p.nc.tol, p.nc.imax, p.sw.newt, p.fuha.blss
%
% See also fsolext, nloop, nlooppde, genamat, getder, stanparam
if (p.fsol.fsol==2 || p.fsol.fsol==3)
[u,res,iter,Gu,Glam]=fsolext(p,u1,ds,varargin); return; end
if ~isempty(varargin); au0=u2au(p,varargin{1},1);
else au0=u2au(p,p.u,1); % last point on branch, don't change!
end
iter=0; u=u1; idx=find(u(1:p.nu)<0); u(idx)=0; % idx, pause
alpha=1;
try almin=p.nc.almine; catch almin=0.2; end % minimal damping, try-catch for backward-comp.
r=resi(p,u); res0=norm(r,p.sw.norm); res=res0; % the residual
[Gu,Glam]=getder(p,u,r); % the derivatives
if(resp.nc.tol && iter0)
p1=p.tau'*[p.sol.xi*aud(1:p.nu);p.sol.xiq*aud(p.nu+1:p.nu+p.nc.nq);(1-(p.sol.xi+p.sol.xiq)/2)*aud(p.nu+p.nc.nq+1)]-ds;
else
p1=p.tau'*[p.sol.xi*aud(1:p.nu);(1-p.sol.xi)*aud(p.nu+p.nc.nq+1)]-ds;
end
[upd,p]=p.fuha.blss(amat,[r;p1],p); stepok=0;
while(stepok==0 && alpha>almin)
au1=au-alpha*upd; u1=au2u(p,au1,1); idx=u1<0; u1(idx==1)=0;
r=resi(p,u1); res=norm(r,p.sw.norm);
if(res0); % if chord, then now get derivatives at new point!
[Gu,Glam]=getder(p,u,r); end
if(p.sw.verb>1); if(alpha<1) fprintf('\nnloopext: damp alpha=%g, res=%g, ds=%g\n',...
alpha,res,ds); end; end; % inform user if damping was used
res=res0;
JBDMV-matlab/acdc/sGjac.m 0000755 0001750 0001750 00000000612 14122154224 013066 0 ustar hu hu function Gu=sGjac(p,u) % 'Jacobian' for AC dead core
a=afu(p,u); par=u(p.nu+1:end); lam=par(1); ga=par(2); c=par(3); n=p.np;
u=u(1:p.nu); up=max(u,p.Jdel); % keeping u away from 0 for singular part
fusing=ga*up.^(ga-1); % the singular part,
fu=(lam-a)+2*(a+1).*u-3*u.^2-lam*fusing; Fu=spdiags(fu,0,n,n);
Gu=c*p.mat.K-p.mat.M*Fu+p.nc.sf*p.mat.Q; % build Jac from bulk and BCs JBDMV-matlab/acdc/nloop.m 0000644 0001750 0001750 00000004633 14045264630 013201 0 ustar hu hu function [u,res,iter,Gu,Glam,p]=nloop(p,u1)
% NLOOP: Newton loop with damping according to settings,
% see in particular stanparam.
%
% [u,res,iter,Gu,Glam,p]=nloop(p,u1)
%
% * u1 is initial guess
% * Returns resulting u, res=residuum, iter=number of iterations, Gu,Glam
% derivatives of G.
% * Note: returns full vector u (including auxiliary variables)!
% * Important settings: p.fsol.fsol, p.nc.tol, p.nc.imax, p.sw.newt, p.fuha.lss
%
% See also au2u, fsol, nloopext, getder, resi, stanparam
if (p.fsol.fsol==1 || p.fsol.fsol==3)
[u,res,iter,Gu,Glam]=fsol(p,u1); return; end
if p.sw.newt==2; % stopping crit based on stepsize in u, not residual
[u,res,iter,Gu,Glam,p]=nloopnleq1(p,u1); return; end
try almin=p.nc.almin; catch almin=0.1; end % minimal damping, try-catch for backward-comp.
u=u1; idx=find(u(1:p.nu)<0); u(idx)=0; %idx, pause
alpha=1; iter=0;
r=resi(p,u); res0=norm(r,p.sw.norm); % (starting) residual
Gu=getGu(p,u,r); % derivatives (needed for next tangent, thus always computed)
if res0p.nc.tol && iteralmin)
au1=u2au(p,u,1) ... % au1=[upde,actaux,lam]
-alpha*[upd;0]; % the newton step, no change in primary param.
u1=au2u(p,au1,1); u1(1:p.nu)=max(u1(1:p.nu),0);
%idx=find(u1(1:p.nu)<0); u1(idx)=0; % idx, pause,
r=resi(p,u1); res=norm(r,p.sw.norm); % residual
if resp.nc.tol && stepok
% some postprocessing
if(p.sw.newt==0); Glam=getGlam(p,u,r); % Newton, only update Glam
else [Gu,Glam]=getder(p,u,r); end % chord, update Gu,Glam
if(p.sw.verb>1); if(alpha<1); % inform user about damping ...
fprintf('nloop: damp alpha=%g, res=%g\n', alpha,res0); end; end;
res=res0; % return res of u0=best u
JBDMV-matlab/acdc/mypsol2D.m 0000664 0001750 0001750 00000000377 14044527545 013574 0 ustar hu hu function mypsol2D(dir,pt,v);
p=loadp(dir,pt); y2=max(max(p.u(1:p.nu)),1);
plotsol(p); hold on;
m=round(p.nu/2); p.u(m-2:m+4)', min(p.u(1:p.nu))
idx=find(p.u(1:p.nu)<1e-6); x=getpte(p);
plot3(x(1,idx),x(2,idx),p.u(idx),'*m');
nola; view(v); getlam(p) JBDMV-matlab/acdc/nodalf.m 0000755 0001750 0001750 00000000164 14044240520 013302 0 ustar hu hu function f=nodalf(p,u) % nodal nonlinearity f, called in sG and e2rs
par=u(p.nu+1:end); u=u(1:p.nu); f=u.^par(2);
JBDMV-matlab/acdc/sG.m 0000755 0001750 0001750 00000000337 14122154213 012412 0 ustar hu hu function r=sG(p,u) % PDE rhs
a=afu(p,u); par=u(p.nu+1:end); lam=par(1); ga=par(2); c=par(3);
u=u(1:p.nu); up=max(u,0);
fti=-u.*(u-1).*(u-a); f=fti-lam*(up.^ga-u);
r=c*p.mat.K*u-p.mat.M*f+p.nc.sf*(p.mat.Q*u-p.mat.G); JBDMV-matlab/acdc/e2rs.m 0000755 0001750 0001750 00000000441 14045265307 012723 0 ustar hu hu function [p,idx]=e2rs(p,u) % classical elements2refine selector as in pdejmps
par=u(p.nu+1:end); u=u(1:p.nu); c=1; a=0; fv=0; %nodalf(p,u);
%E=p.pdeo.errorInd(u,c,a,fv);
E=abs(p.mat.K*p.u(1:p.nu)); E=p.mat.p2c*E;
p.sol.err=max(max(E));
idx=p.pdeo.selectElements2Refine(E,p.nc.sig); JBDMV-matlab/acdc/cmds1D.m 0000664 0001750 0001750 00000003042 14122155312 013150 0 ustar hu hu %% 1D Allen-Cahn dead core
close all; keep pphome;
%% init
p=[]; par=[0.1 0.85 0.01 1 0.2 8 0.6];
% lam, ga=1/m, c, cubic, pa, pk pm (ampl., wave-nr, mean)
dim=1; sw=1; p=acinit(p,1,200,par,dim,sw); p=setfn(p,'0'); plotsol(p);
p.usrlam=[18 25]; p.Jdel=1e-8; % u_+=max(u,Jdel) in sGjac
p=cont(p,30);
%% center DC
p=swibra('0','bpt1','b1',-0.1); p.nc.dsmax=5; p=cont(p,100);
% p=loadp('b1','pt20','b1r'); p=oomeshadac(p,'ngen',5,'sig',0.5); % meshada
%% further DC branches
p=swibra('0','bpt2','b2',0.1); p.nc.almin=0.1; p.sw.jac=1; p=cont(p,100); % asymetric
p=swibra('0','bpt3','b3a',0.1); p=cont(p,80); % center DC, 2 humps, unstable
p=swibra('0','bpt3','b3b',-0.1); p=cont(p,60); % 2 DCs, center hump, partly stable
%% BD, min
f=3; c=9; figure(f); clf; plotbra('0',f,c,'cl','k','lsw',0);
plotbra('b1',f,c,'cl','b');
axis([0 8 0 1.1]); ylabel('min(u)');
%% BD, L^2
f=3; c=0; figure(f); clf; plotbra('0',f,c,'cl','k','lsw',0);
plotbra('b1',f,c,'cl','b','lab', [60]); plotbra('b2',f,c,'cl','m','lab',[25 60]);
plotbra('b3a',f,c,'cl',p2pc('r1'),'lab',[10 40]);
plotbra('b3b',f,c,'cl',p2pc('r2'),'lab',55); ylabel('||u||_2');
axis([0 40 0.1 1.8]);
%% sol plot;
mypsol('b1','pt20'); pause; mypsol('b1','pt60'); pause; mypsol('b2','pt25'); pause;
mypsol('b2','pt60'); pause; mypsol('b3a','pt10'); pause; mypsol('b3a','pt40'); pause;
mypsol('b3b','pt55');
%% cont in c: DCs lift up
p=swiparf('b3b','pt40','b3b-c',3); p.sol.ds=0.01; p.nc.dsmax=0.1; p=cont(p,15);
%% 2ndary bifs,
p=swibra('b3b','bpt1','b3-1'); p=cont(p,50); JBDMV-matlab/acdc/diskpdeo.m 0000644 0001750 0001750 00000000700 14045270220 013633 0 ustar hu hu % mod of diskpdeo to ellipse; r=radius, nr=#of boundary points, del=ellipticity
classdef diskpdeo< pde
methods(Access = public)
function o=diskpdeo(r,nr,del) % constructor
o.grid=grid2D; s=linspace(0,2*pi,nr);
o.grid.freeGeometry([r*cos(s)*del;r*sin(s)]); o.fem=lagrange12D;
end
end
methods(Access = protected) % only here since required by pde-class
function r=df(~,~,~); r=0; end % rather use p.fuha.sG in pderesi
end
end
JBDMV-matlab/acdc/oosetfemops.m 0000755 0001750 0001750 00000000765 14044240520 014411 0 ustar hu hu function p=oosetfemops(p)
gr=p.pdeo.grid; fem=p.pdeo.fem; % just introduced as shorthands
[K,M,~]=fem.assema(gr,1,1,1); p.mat.K=K; p.mat.M=M; % indep. of BC
bcl=gr.robinBC(1,1);
gr.makeBoundaryMatrix(bcl); % intermediate step before assembling BC matr.
[Q,G,~,~]=fem.assemb(gr); p.mat.Q=Q; p.mat.G=G; % the BC matrices
p.nc.sf=1e3; % stiff spring constant for DBC via Robin-BC
p.mat.E=center2PointMatrix(gr); % to map element different. matrices to nodal ones
p.mat.p2c=point2CenterMatrix(gr); JBDMV-matlab/acdc/acinit.m 0000755 0001750 0001750 00000002135 14045266042 013316 0 ustar hu hu function p=acinit(p,lx,nx,par,dim,sw,del) % init with prep. of mesh-adaption
p=stanparam(p); screenlayout(p); p.nc.neq=1;
p.sw.sfem=-1; p.fuha.sG=@sG; p.fuha.sGjac=@sGjac; p.fuha.e2rs=@e2rs;
switch dim;
case 1; pde=stanpdeo1Db(0,lx,lx/nx); p.plot.axis=[0 1 0 1];
case 2; pde=diskpdeo(lx,nx,del); % refine all longest edges for more symm.
rlongs=pde.grid.rlong; pde.grid.rlong=1; idx=1:size(pde.grid.t,2);
pde.grid.refineMesh(idx); pde.grid.rlong=rlongs;
case 3; pde=stanpdeo2D(lx,lx/2,lx/nx);
end
p.pdeo=pde; p.np=pde.grid.nPoints; p.nu=p.np; p.sol.xi=1/(p.nu);
p.mesh.bp=pde.grid.p; p.mesh.bt=pde.grid.t; p.mesh.be=pde.grid.e;
[po,t,e]=getpte(p); x=po(1,:)'; x0=0.5; a=0.0001; b=16*(1-a); u0=a+b*(x-x0).^4;
switch sw;
case 0; u0=0*u0;
case 1; u0=1+0*x;
case 2; y=po(2,:)'; r=x.^2+y.^2; b=1-a; u0=a+b*r.^2;
end
p.u=[u0; par']; % initial guess, and parameters
p.usrlam=[]; p=setfemops(p); p.mesh.maxt=100;
p.plot.auxdict={'\lambda','\gamma'};
p.nc.ilam=1; p.nc.lammin=-1; p.nc.lammax=50;
p.sol.ds=0.1; p.nc.dsmax=2; p.nc.sf=1000; p.nc.del=1e-2; JBDMV-matlab/acdc/afu.m 0000664 0001750 0001750 00000000405 14122134576 012621 0 ustar hu hu function a=afu(p,u) % function a(x) for acdc
par=u(p.nu+1:end);pa=par(5);pk=par(6);pm=par(7); % amplitude, wave-nr, mean
po=getpte(p); if 1; x=po(1,:)'; else x=po'; end
if size(po,1)==2; y=po(2,:)'; x=sqrt(x.^2+y.^2); end
x2=max(x); a=pm+pa*cos(pk*pi*x/x2); JBDMV-matlab/acdc/cmds2D.m 0000664 0001750 0001750 00000001633 14122155723 013163 0 ustar hu hu %% 2D dead core
close all; keep pphome;
%% init,
p=[]; par=[-0.1 0.85 0.05 1 0.3 8 0.6]; dim=2; sw=1; del=1.2;
p=acinit(p,1,120,par,dim,sw,del); p=setfn(p,'2'); p.Jdel=1e-8; plotsol(p);
p.nc.ngen=1; p.nc.sig=0.6; p.nc.maxt=12000; % mesh-adapt.settings
p.nc.lammax=40; p.nc.tol=1e-6; % slightly relax tol
%% first 4 BPs, then cont with larger ds
p=findbif(p,4); p.nc.dsmax=5; p.nc.dlammax=2; p=cont(p,20);
%% center DC branch, with meshada
p=swibra('2','bpt1','c1',-0.1); p.nc.amod=10; p=cont(p,50); % center DC
%% 2nd branch
p=swibra('2','bpt2','c2',0.01); p.nc.sig=0.8; p.nc.amod=10; p=cont(p,100);
%% BD, L^2
f=3; c=0; figure(f); clf;
plotbra('2',f,c,'cl','k'); plotbra('c1',f,c,'cl','b','lab',[49]);
plotbra('c2',f,c,'cl','m','lab',[100]); ylabel('||u||_2');
%% sol plots;
v=[0 90];
mypsol2D('c1','pt49',v); pause; mypsol2D('c2','pt100',v); pause;
plotsol('c2','pt100',1,1,3); nola; colormap parula; JBDMV-matlab/schnakdc/mypsol.m 0000644 0001750 0001750 00000000420 14044766511 014253 0 ustar hu hu function mypsol(dir,pt);
p=loadp(dir,pt); y2=max(max(p.u(1:p.nu)),1);
plotsol(p); %axis([0 1 0 y2]);
hold on;
m=round(p.nu/2); p.u(m-2:m+4)', min(p.u(1:p.nu))
idx=find(p.u(1:p.nu)<1e-8); x=getpte(p);
plot(x(idx),p.u(idx),'*m'); nola
title([dir '/' pt]); xticks([]);
JBDMV-matlab/schnakdc/nloopext.m 0000644 0001750 0001750 00000005167 14045302410 014576 0 ustar hu hu function [u,res,iter,Gu,Glam,p]=nloopext(p,u1,ds,varargin)
% NLOOPEXT: newton loop for arclength extended system.
%
% [u,res,iter,Gu,Glam]=nloopext(p,u1,ds) : use p.u and u1 for initial arclength step
% [u,res,iter,Gu,Glam]=nloopext(p,u1,ds,u) : use u and u1 for initial arclength step
%
% * u1 is initial guess, ds initial stepsize usually taken from p-structure.
% * Returns u, res=residuum, iter=number of iterations, Gu,Glam=derivatives.
% * Note: returns full vector u (including auxiliary variables)!
% * Important settings: p.fsol.fsol, p.nc.tol, p.nc.imax, p.sw.newt, p.fuha.blss
%
% See also fsolext, nloop, nlooppde, genamat, getder, stanparam
if (p.fsol.fsol==2 || p.fsol.fsol==3)
[u,res,iter,Gu,Glam]=fsolext(p,u1,ds,varargin); return; end
if ~isempty(varargin); au0=u2au(p,varargin{1},1);
else au0=u2au(p,p.u,1); % last point on branch, don't change!
end
iter=0; alpha=1; u=u1; %if p.zc; u(1:p.np)=max(u(1:p.np),0); end
try almin=p.nc.almine; catch almin=0.2; end % minimal damping, try-catch for backward-comp.
try alpha=p.nc.almax; catch alpha=1; end
r=resi(p,u); res0=norm(r,p.sw.norm); res=res0; % the residual
[Gu,Glam]=getder(p,u,r); % the derivatives
if(resp.nc.tol && iter0)
p1=p.tau'*[p.sol.xi*aud(1:p.nu);p.sol.xiq*aud(p.nu+1:p.nu+p.nc.nq);(1-(p.sol.xi+p.sol.xiq)/2)*aud(p.nu+p.nc.nq+1)]-ds;
else
p1=p.tau'*[p.sol.xi*aud(1:p.nu);(1-p.sol.xi)*aud(p.nu+p.nc.nq+1)]-ds;
end
[upd,p]=p.fuha.blss(amat,[r;p1],p); stepok=0; upd=real(upd);
while(stepok==0 && alpha>almin)
au1=au-alpha*upd; u1=au2u(p,au1,1); u1(1:p.np)=max(u1(1:p.np),0);
r=resi(p,u1); res=norm(r,p.sw.norm);
if(res0); % if chord, then now get derivatives at new point!
[Gu,Glam]=getder(p,u,r); end
if(p.sw.verb>1); if(alpha<1) fprintf('\nnloopext: damp alpha=%g, res=%g, ds=%g\n',...
alpha,res,ds); end; end; % inform user if damping was used
res=res0;
JBDMV-matlab/schnakdc/sGjac.m 0000644 0001750 0001750 00000000425 14117716372 013765 0 ustar hu hu function Gu=sGjac(p,u)
par=u(p.nu+1:end); n=p.np;
[f1u,f1v,f2u,f2v]=nodaljac(p,u); % (nodal) Jacobian of 'nonlin'
Fu=[[spdiags(f1u,0,n,n),spdiags(f1v,0,n,n)];
[spdiags(f2u,0,n,n),spdiags(f2v,0,n,n)]];
Ks=p.mat.K; c=par(2); d=par(3);
Gu=[c*Ks 0*Ks; 0*Ks d*Ks]-p.mat.M*Fu; JBDMV-matlab/schnakdc/nloop.m 0000644 0001750 0001750 00000004625 14045302354 014062 0 ustar hu hu function [u,res,iter,Gu,Glam,p]=nloop(p,u1)
% NLOOP: Newton loop with damping according to settings,
% see in particular stanparam.
%
% [u,res,iter,Gu,Glam,p]=nloop(p,u1)
%
% * u1 is initial guess
% * Returns resulting u, res=residuum, iter=number of iterations, Gu,Glam
% derivatives of G.
% * Note: returns full vector u (including auxiliary variables)!
% * Important settings: p.fsol.fsol, p.nc.tol, p.nc.imax, p.sw.newt, p.fuha.lss
%
% See also au2u, fsol, nloopext, getder, resi, stanparam
np=p.np;
if (p.fsol.fsol==1 || p.fsol.fsol==3)
[u,res,iter,Gu,Glam]=fsol(p,u1); return; end
if p.sw.newt==2; % stopping crit based on stepsize in u, not residual
[u,res,iter,Gu,Glam,p]=nloopnleq1(p,u1); return; end
try almin=p.nc.almin; catch almin=0.05; end % minimal damping, try-catch for backward-comp.
try alpha=p.nc.almax; catch alpha=1; end
u=u1; iter=0; %if p.zc; u(1:np)=max(u1(1:np),0); end
r=resi(p,u); res0=norm(r,p.sw.norm); % (starting) residual
Gu=getGu(p,u,r); % derivatives (needed for next tangent, thus always computed)
if res0p.nc.tol && iteralmin)
au1=u2au(p,u,1) ... % au1=[upde,actaux,lam]
-alpha*[upd;0]; % the newton step, no change in primary param.
u1=au2u(p,au1,1); u1(1:np)=max(u1(1:np),0);
r=resi(p,u1); res=norm(r,p.sw.norm); % residual
if resp.nc.tol && stepok
% some postprocessing
if(p.sw.newt==0); Glam=getGlam(p,u,r); % Newton, only update Glam
else [Gu,Glam]=getder(p,u,r); end % chord, update Gu,Glam
if(p.sw.verb>1); if(alpha<1); % inform user about damping ...
fprintf('nloop: damp alpha=%g, res=%g\n', alpha,res0); end; end;
res=res0; % return res of u0=best u
JBDMV-matlab/schnakdc/mypsol2D.m 0000664 0001750 0001750 00000000417 14044723775 014456 0 ustar hu hu function mypsol2D(dir,pt,v);
p=loadp(dir,pt); y2=max(max(p.u(1:p.nu)),1);
plotsol(p); colormap parula; hold on;
m=round(p.nu/2); p.u(m-2:m+4)', min(p.u(1:p.nu))
idx=find(p.u(1:p.nu)<1e-8); x=getpte(p);
plot3(x(1,idx),x(2,idx),p.u(idx),'*m');
nola; view(v); getlam(p) JBDMV-matlab/schnakdc/myhopl.m 0000664 0001750 0001750 00000000760 14045002370 014234 0 ustar hu hu function myhopl(dir,pt,varargin)
hoplotf(dir,pt,1,1); xlabel('t');
p=loadp(dir,pt); x=getpte(p); x0=x(p.x0i);
legend(['u|_{x=' mat2str(x0) '}']);
figure(1); colormap parula;
shading interp; view([40 60]);
xlabel('x'); ylabel('t'); title([dir '/' pt]);
if nargin>2; zticks(varargin{1}); end
return
tl=p.hopf.tl; n=p.np; hold on
for i=1:tl;
u=p.hopf.y(1:n,1); t=p.hopf.t(i)*p.hopf.T;
idx=find(u(1:n)<1e-6); ndc=length(idx); idx
plot3(x(idx),t*ones(1,ndc),u(idx)','*m');
end JBDMV-matlab/schnakdc/sG.m 0000664 0001750 0001750 00000000444 14117717573 013316 0 ustar hu hu function r=sG(p,u) % Schnakenberg
par=u(p.nu+1:end); lam=par(1); c=par(2); d=par(3); m=par(4);
np=p.np; u1=u(1:np); u2=u(np+1:2*np); if p.zc2; u1=max(u1,0); end
uo=u1.^(1/m); f1=-uo+u1.^2.*u2; f2=lam-u1.^2.*u2; f=[f1; f2];
K1= p.mat.K; K=[c*K1 0*K1; 0*K1 d*K1]; r=K*u(1:p.nu)-p.mat.M*f; JBDMV-matlab/schnakdc/e2rs.m 0000644 0001750 0001750 00000000337 13665510751 013613 0 ustar hu hu function [p,idx]=e2rs(p,u) % classical elements2refine selector as in pdejmps
u=u(1:p.np); c=1; a=0; fv=0; %nodalf(p,u);
E=p.pdeo.errorInd(u,c,a,fv); p.sol.err=max(max(E));
idx=p.pdeo.selectElements2Refine(E,p.nc.sig); JBDMV-matlab/schnakdc/cmds1D.m 0000664 0001750 0001750 00000003710 14122156216 014042 0 ustar hu hu close all; keep pphome;
%% DC-Schnakenberg; don't choose m too large (non-conv at small lam)
p=[]; par=[3.2 1 100 1.5]; nx=300; lx=12; % [lam c d m]
p=schnakinit(p,lx,nx,par); p=setfn(p,'tr');
p.zc=1; p.zc2=0; % zero-correction in Newton and in sG
p=findbif(p,6); p.sw.bifcheck=0; p.nc.dsmax=1; p=cont(p,20);
%%
p=swibra('tr','bpt1','T1',0.1); p=cont(p,50);
p=swibra('tr','bpt2','T2',0.1); p=cont(p,50);
p=swibra('tr','bpt3','T3',0.1); p=cont(p,50);
%% steady from branch with DC
p=swibra('T1','bpt1','T1-1',0.05); p.nc.dsmax=0.1; p.nc.tol=2e-8; p=cont(p,120);
%% Hopf
aux.dlam=0; aux.tl=40; p=hoswibra('T1','hpt1',0.1,4,'T1-h1',aux); p.file.smod=2;
p.nc.tol=1e-6; p.hopf.flcheck=0; p.fuha.outfu=@hobra; p.usrlam=[];
p.x0i=p.np; p=setbel(p,1,1e-6,10,@lss); p=cont(p,10);
%% Hopf
aux.dlam=0; aux.tl=50; p=hoswibra('T3','hpt1',0.1,4,'T3-h1',aux); p.file.smod=2;
p.nc.tol=1e-6; p.hopf.flcheck=0; p.fuha.outfu=@hobra; p.usrlam=[];
p.x0i=p.np; p=setbel(p,1,1e-6,10,@lss); p=cont(p,10);
%% plot BD
fnr=3; figure(fnr); clf; c=6; plotbra('tr','pt50',fnr,c,'cl','k','fp',1);
plotbra('T1',fnr,c,'cl','b', 'lab',25);
plotbra('T1-H1','pt8',fnr,c,'cl',p2pc('o1'),'lab',[2 8]);
plotbra('T2',fnr,c,'cl',p2pc('b1'),'lab',30);
plotbra('T3',fnr,c,'cl',p2pc('g1'),'lab',40);
plotbra('T3-2',fnr,c,'cl',p2pc('r2'),'lab',15);
plotbra('T1-1','pt140',fnr,c,'cl',p2pc('r3'), 'lab',[40 100],'lp',105);
plotbra('T3-H1','pt8',fnr,c,'cl',p2pc('o2'),'lab',[8]);
ylabel('max(u)');
%% soln plot, steady
mclf(1);
mypsol('T1','bpt1'); pause; mypsol('T1','hpt1'); pause; mypsol('T1','pt25'); pause;
mypsol('T1-1','pt40'); pause; mypsol('T1-1','pt100'); pause;
mypsol('T3','bpt2'); pause; mypsol('T3','hpt1'); pause; mypsol('T3','pt40'); pause;
mypsol('T2','bpt1'); pause; mypsol('T2','hpt1'); pause; mypsol('T2','pt30'); pause;
mypsol('T3-2','pt15');
%% soln plots, POs
myhopl('T1-H1','pt2',[1 2]); pause; myhopl('T1-H1','pt8'); pause
myhopl('T3-H1','pt8'); JBDMV-matlab/schnakdc/nodaljac.m 0000664 0001750 0001750 00000000537 14117716357 014522 0 ustar hu hu function [f1u,f1v,f2u,f2v]=nodaljac(p,u) % Jacobian for Schnakenberg
par=u(p.nu+1:end); m=par(4); ga=1/m;
u1=u(1:p.nu/2); u2=u(p.nu/2+1:p.nu); del=p.nc.del;
iz=find(u1<2*p.nc.del);
fsd=ga*u1.^(ga-1);
fsd(iz)=((u1(iz)+del).^ga-u1(iz).^ga)/del; % singular derivative, FD-like
f1u=-fsd+2*u1.*u2;
f1v=u1.^2; f2u=-2*u1.*u2; f2v=-f1v; % regular part JBDMV-matlab/schnakdc/oosetfemops.m 0000644 0001750 0001750 00000000307 13665503537 015304 0 ustar hu hu function p=oosetfemops(p)
%gr=p.pdeo.grid; fem=p.pdeo.fem;
[p.mat.K,M,~]=p.pdeo.fem.assema(p.pdeo.grid,1,1,1); % FEM/mass matrices
p.mat.M=[M 0*M; 0*M M]; p.mat.p2c=point2CenterMatrix(p.pdeo.grid);
JBDMV-matlab/schnakdc/honloopext.m 0000664 0001750 0001750 00000010274 14045303162 015127 0 ustar hu hu function [y,T,lam,a,res,iter,A,p,jac]=honloopext(p,y1,T1,lam1,a1,ds,varargin)
% honloopext: Newton loop for Hopf in arclength
%
% updated for cases
% without t-transl. phase-condition: p.hopf.pc=0 (default 1)
% with fixed period T: p.hopf.freeT=0 (default 1)
% for compatibility with old data, these, and associated mods (in particular
% relating to p.hopf.ilam), are read out via try-catch
if p.fsol.fsol==1; % use fsolve
[y,T,lam,res,iter,A,p,jac]=fsolext(p,y1,T1,lam1,ds); return; end
try; freeT=p.hopf.freeT; catch; freeT=1; end % check if T is free (default)
try; pcsw=p.hopf.pc; catch; pcsw=1; end % check if t-PC is used (default)
try; nqh=p.hopf.nqh; catch; nqh=0; end
try; nqh2=length(p.hopf.ilam); catch; nqh2=0; end
try; a=p.u(p.nu+p.hopf.ilam); catch; a=0; end;
iter=0; y=max(y1,0); T=T1; lam=lam1; tol=p.nc.tol;
p=setlam(p,lam); % set lam to predictor
nu=p.nu; tl=p.hopf.tl;
alpha=1; try; almin=p.nc.almin; catch; almin=0.5; end % minimal damping
[f,jac,f_T,f_lam,f_a]=tomassempbc(p,y,T,lam);
pc=[]; pc_y=[]; if pcsw; [pc,pc_y]=hopc(p,y,T,lam); end
qf=[]; qhder=[]; if nqh>0; qf=p.hopf.qfh(p,y); qhder=getqhder(p,y); end;
arcl=hoarcl(p,y,T,lam,ds,a);
F=[f; pc; arcl; qf]; res0=norm(F,p.sw.norm); res=res0;
if p.sw.verb>1;
if pcsw; fprintf('initial (||f||, pc, arcl)=(%g,%g,%g)\n',norm(f,p.sw.norm),pc,arcl);
else; fprintf('initial (||f||, arcl)=(%g,%g)\n',norm(f,p.sw.norm),arcl); end
if nqh>0; fprintf('qh=%g\n',qf); end
end
% size(jac), size(qhder), size(p.hopf.tau), pause
A=gethoA(p,jac,f_T,f_lam,pc_y,p.hopf.tau,f_a,qhder);
if(restol && iter=almin)
y1=yv-alpha*upd(1:nu*tl);
y1(1:nu*tl)=max(y1(1:nu*tl),0);
if freeT; tupd=upd(nu*tl+1); T1=T-alpha*tupd;
lam1=lam-alpha*upd(nu*tl+2);
if nqh2>0; a=p.u(p.nu+p.hopf.ilam);
a1=a-alpha*upd(nu*tl+3:end);
p.u(p.nu+p.hopf.ilam)=a1;
end
else; T1=T; lam1=lam-alpha*upd(nu*tl+1); %nqh2, pause
if nqh2>0; a=p.u(p.nu+p.hopf.ilam);
a1=a-alpha*upd(nu*tl+2:end); p.u(p.nu+p.hopf.ilam)=a1;
end
end
ytom=reshape(y1(1:nu*tl),nu,tl); % convert to tom-input
f=tomassempbc(p,ytom,T1,lam1); % setlam(p,lam1) done in tomassempbc!
if pcsw; [pc,pc_y]=hopc(p,ytom,T1,lam1); end
if freeT;
arcl=hoarcl(p,ytom,T1,lam1,ds,a1);
if nqh>0; qf=p.hopf.qfh(p,ytom); end
else
if nqh2>0; p.u(p.nu+p.hopf.ilam)=a; arcl=hoarcl(p,ytom,T1,lam1,ds,a1);
p.u(p.nu+p.hopf.ilam)=a1;
else; arcl=hoarcl(p,ytom,T1,lam1,ds,a1);
end
end
F=[f; pc; arcl]; if nqh>0; F=[F; qf]; end
if p.sw.verb>1;
if nqh==0;
if pcsw; fprintf('al=%g, (||f||, pc, arcl)=(%g,%g,%g)\n',alpha,norm(f,p.sw.norm),pc,arcl);
else; fprintf('al=%g, (||f||, arcl)=(%g,%g)\n',alpha,norm(f,p.sw.norm),arcl);
end
else
if pcsw; fprintf('al=%g, (||f||, pc, arcl,qh)=(%g,%g,%g,%g)\n',alpha,norm(f,p.sw.norm),pc,arcl,qf);
else; fprintf('al=%g, (||f||, arcl, qh)=(%g,%g,%g)\n',alpha,norm(f,p.sw.norm),arcl,qf);
end
end
%if nqh>0; fprintf('qh=%g\n',qf); end
end
res=norm(F,p.sw.norm); % check new tom-resi
if(res1); if(alpha<1) fprintf('\nhonloopext: damp alpha=%g, res=%g, ds=%g\n',...
alpha,res,ds); end; end; % inform user if damping was used
res=res0; try; y=ytom; catch; y=y1; end
JBDMV-matlab/schnakdc/schnakinit.m 0000644 0001750 0001750 00000001763 14045302244 015064 0 ustar hu hu function p=schnakinit(p,dom,nx,par,varargin)
p=stanparam(p); screenlayout(p); p.nc.neq=2; p.sw.sfem=-1;
p.plot.auxdict={'\lambda','c','d','m'}; p.fuha.outfu=@hobra;
switch length(dom)
case 1; lx=dom; p.pdeo=stanpdeo1D(lx,2*lx/nx); p.vol=2*lx;
case 2; ny=round(dom(2)/dom(1)*nx); lx=dom(1); ly=dom(2); p.vol=4*lx*ly;
pde=stanpdeo2D(lx,ly,nx,ny,varargin{1}); p.pdeo=pde; p.plot.pstyle=2;
end
p.np=p.pdeo.grid.nPoints; p.nu=p.np*p.nc.neq; p.nc.neig=10; p.nc.nsteps=50;
p=setfemops(p); p.sol.xi=1/p.nu; p.file.smod=5; p.sw.para=2; p.sw.foldcheck=1;
p.nc.ilam=1; p.sol.ds=-0.1; p.nc.dsmax=0.5; p.nc.lammin=0.1; p.sw.bifcheck=2;
lam=par(1); m=par(4); al=m; bet=1-2*m;
u=lam^al*ones(p.np,1); v=lam^bet*ones(p.np,1);
p.u=[u;v;par']; % initial solution guess with parameters
[po,tr,ed]=getpte(p); p.mesh.bp=po; p.mesh.be=ed; p.mesh.bt=tr; p.nc.ngen=1;
p.plot.bpcmp=6; p.plot.axis='image'; p.plot.cm='parula';
p.nc.resfac=1e-3; p.pm.mst=4; p.nc.del=1e-8; % small delta helpful for FDs JBDMV-matlab/schnakdc/cmds2D.m 0000664 0001750 0001750 00000002304 14122157312 014037 0 ustar hu hu close all; keep pphome;
%% DC-Schnakenberg in semi-lin setting
p=[]; par=[3.2 1 100 1.5]; lx=12; ly=lx/sqrt(3); nx=80; sw.sym=1;
p=schnakinit(p,[lx ly],nx,par,sw); p.nc.dsmax=0.5; p.sol.ds=-0.1;
p=setfn(p,'tr2'); p.nc.lammin=0.1; p.fuha.outfu=@hobra;
p.zc=1; % zero-correction in Newton
p.zc2=0; % zero-correction in sGdc
p.sw.verb=2; p.file.smod=5; p.nc.neig=30; p.nc.dsmax=0.1; p.nc.del=1e-8;
p.plot.bpcmp=6; p.nc.neig=10;
p=cont(p,15); p.sw.bifcheck=0; p.nc.dsmax=1; p=cont(p,15); % continue for plotting d
%% hex via qswibra, continue in both directions
aux=[]; aux.m=3; p0=qswibra('tr2','bpt1',aux); p0.nc.dsmin=0.01; p0.sw.bifcheck=0;
p0.sw.foldcheck=0; p0.nc.dsmax=0.5;
%% select directions and go
p=seltau(p0,2,'h1+',2); p.sol.ds=-0.2; p.nc.tol=1e-6; p=cont(p,50);
p=seltau(p0,2,'h1-',2); p.sol.ds=0.1; p.nc.tol=1e-6; p=cont(p,50);
%% plot BD
fnr=3; figure(fnr); clf; c=6; plotbra('tr2',fnr,c,'cl','k','fp',1);
plotbra('h1+',fnr,c,'cl','b', 'lab',[10 27]);
plotbra('h1-',fnr,c,'cl',p2pc('b3'), 'lab',[5 27]);
ylabel('max(u)');
%%
mypsol2D('h1+','pt10',[0,90]); pause; mypsol2D('h1+','pt19',[0,90]); pause
mypsol2D('h1-','pt5',[0,90]); pause; mypsol2D('h1-','pt20',[0,90]); JBDMV-matlab/sh35disk/qxon.m 0000664 0001750 0001750 00000000161 14040302422 013557 0 ustar hu hu function p=qxon(p) % switch on x-phase condition
p.nc.nq=1; p.fuha.qf=@qx; p.fuha.qfder=@qxder; p.nc.ilam=[1 4]; JBDMV-matlab/sh35disk/qyder.m 0000664 0001750 0001750 00000000177 14040311762 013734 0 ustar hu hu function qu=qyder(p,u) % derivative of transl.phase condition in y
qy=(p.mat.Dy(1:p.np,1:p.np)*p.u(1:p.np))'; qu=[qy, 0*qy]; JBDMV-matlab/sh35disk/qxyron.m 0000664 0001750 0001750 00000000173 14040302513 014136 0 ustar hu hu function p=qxyron(p) % switch on x-phase condition
p.nc.nq=3; p.fuha.qf=@qxyr; p.fuha.qfder=@qxyrder; p.nc.ilam=[1 4 5 6]; JBDMV-matlab/sh35disk/qxyon.m 0000664 0001750 0001750 00000000166 14040302455 013763 0 ustar hu hu function p=qxyon(p) % switch on x-phase condition
p.nc.nq=2; p.fuha.qf=@qxy; p.fuha.qfder=@qxyder; p.nc.ilam=[1 4 5]; JBDMV-matlab/sh35disk/sGjac.m 0000755 0001750 0001750 00000001133 14040315767 013642 0 ustar hu hu function Gu=sGjac(p,u)
[f1u,f1v,f2u,f2v]=njac(p,u); n=p.nu/2; % the 'local' derivatives
par=u(p.nu+1:end); u=u(1:p.nu); sx=par(4); sy=par(5); s=par(6);
Fu=[[spdiags(f1u,0,n,n),spdiags(f1v,0,n,n)];
[spdiags(f2u,0,n,n),spdiags(f2v,0,n,n)]];
K=p.mat.Ks; Ms=p.mat.M(1:p.np, 1:p.np); Ksys=[[0*K -K];[K Ms]];
Gu=Ksys+sx*p.mat.Dx+sy*p.mat.Dy+s*p.mat.Krot-p.mat.M*Fu;
end
function [f1u,f1v,f2u,f2v]=njac(p,u)
n=p.nu/2; u1=u(1:n); par=u(p.nu+1:end); lam=par(1); nup=par(2);
try c2=par(7); catch; c2=0; end
ov=ones(n,1);
f1u=(lam-1)*ov+3*nup*u1.^2-5*u1.^4+2*c2*u1; f1v=-2*ov;
f2u=0*ov; f2v=0*ov;
end JBDMV-matlab/sh35disk/qrot.m 0000664 0001750 0001750 00000000145 14040312210 013555 0 ustar hu hu function q=qrot(p,u) % rotational phase condition
n=p.np; q=(p.mat.Krot(1:n,1:n)*p.u(1:n))'*u(1:n); JBDMV-matlab/sh35disk/qxyrder.m 0000664 0001750 0001750 00000000355 14040302311 014272 0 ustar hu hu function qu=qxyrder(p,u) % derivative of all three phase conditions
np=p.np; uold=p.u(1:np); u0x=p.mat.Dx(1:np,1:np)*uold; u0y=p.mat.Dy(1:np,1:np)*uold;
u0phi=p.mat.Krot(1:np,1:np)*uold;
qu=[u0x' 0*u0x';u0y' 0*u0y';u0phi' 0*u0phi']; JBDMV-matlab/sh35disk/shJ.m 0000664 0001750 0001750 00000000312 14035423621 013325 0 ustar hu hu function J=shJ(p,u) % energy for 2-3 SH
par=u(p.nu+1:end); lam=par(1); nup=par(2); n=p.nu/2; u1=u(1:n); u2=u(n+1:2*n);
dens=(u1+u2).^2/2-lam*u1.^2/2-nup*u1.^3/3+u1.^4/4; J=sum(p.mat.M(1:n,1:n)*dens); JBDMV-matlab/sh35disk/qyon.m 0000664 0001750 0001750 00000000161 14040302414 013561 0 ustar hu hu function p=qyon(p) % switch on y-phase condition
p.nc.nq=1; p.fuha.qf=@qy; p.fuha.qfder=@qyder; p.nc.ilam=[1 5]; JBDMV-matlab/sh35disk/sG.m 0000755 0001750 0001750 00000001025 14040317655 013162 0 ustar hu hu function r=sG(p,u) % rhs for SH; first split u into parameters and PDE part
n=p.np; par=u(p.nu+1:end); u=u(1:p.nu); lam=par(1); nup=par(2);
sx=par(4); sy=par(5); s=par(6); % Lagrange-pars for the phase-conditions
u1=u(1:n); u2=u(n+1:2*n); % the 2 components
f1=(lam-1)*u1-2*u2+nup*u1.^3-u1.^5; f2= 0*u2; f=[f1;f2]; % the 'nonlin'
K=p.mat.Ks; Ms=p.mat.M(1:n, 1:n); % scalar Lapl. and mass matrices
Ksys=[[0*K -K];[K Ms]]; % for linear part of rhs
r=Ksys*u-p.mat.M*f+s*p.mat.Krot*u+sx*p.mat.Dx*u+sy*p.mat.Dy*u; % rhs, +PCs JBDMV-matlab/sh35disk/cmds1.m 0000664 0001750 0001750 00000007760 14131115736 013630 0 ustar hu hu %% sh35 on (half) disk; init, and saving of u(eps=0) as 0h/pt1
% uses quadratic FEM, hence first run 'sethofem.m' in demos/hofem
p=[]; dswitch=2; lam=-0.01; nu=1.4; q=1; par=[lam nu q 0 0 0]; rad=14; nref=5;
p=shinit(p,dswitch,par,rad,nref); p=setfn(p,'0h'); % init, and set out-dir
p.sol.ds=0.1; p=cont(p,1); % one dummy-step (to eps=0)
%% bif directions
aux.m=6; aux.besw=0; p0=cswibra('0h','pt1',aux); p0.nc.eigref=-1;
%% radial
p=gentau(p0,[0 0 0 0 1],'r'); p.sw.bifcheck=0; p.u(p.nu+2)=2;
p.sol.ds=-0.01; p.nc.dsmax=0.02; p=cont(p,5); p=qyon(p); p=cont(p,75);
p.nc.dsmax=0.04; p=cont(p,50); % increase max stepsize and cont further
%% D_4^- mode
p=gentau(p0,[0 0 0 0 0 1],'d4m'); p.sw.bifcheck=0; p.u(p.nu+2)=2;
p.sol.ds=-0.01; p.nc.dsmax=0.02; p=cont(p,5); p=qyon(p); p=cont(p,50);
p.nc.dsmax=0.04; p=cont(p,40);
%% branch plot
f=11; figure(f); clf; plotbra('r','pt130','cmp',11,'cl','k','lab',[10 120],'wnr',f,'fms',0,'ms',0);
plotbra('d4m','pt140','cmp',11,'cl','r','lab',[20 60 90 130],'wnr',f,'fms',0,'ms',0); ylabel('||u||_2');
%% soln plot
plotsol('r','pt10'); solstyle('10'); pause; plotsol('r','pt120'); solstyle('120'); pause
plotsol('d4m','pt20'); solstyle('20','r'); pause; plotsol('d4m','pt60'); solstyle('60','r'); pause;
plotsol('d4m','pt90'); solstyle('90','r'); pause; plotsol('d4m','pt130'); solstyle('130','r');
%% ------- wall mode (daisy) with nu=1.4; initially with bifdetec, then without ------------
p=gentau(p0,[0 0 1],'wall'); p.nc.eigref=-1; p.sw.bifcheck=1; p.sol.ds=0.01; p.nc.dsmax=0.02;
p=cont(p,10); p.sw.bifcheck=0; p=cont(p,20);
%% daisysnake
clf(2); p=swibra('wall','bpt1','dsnake'); p.nc.dsmax=0.01; p.sw.bifcheck=0; p=cont(p,200);
%% branch plot
f=11; figure(f); clf; plotbra('wall','cmp',11,'cl','b','lab',[20],'wnr',f,'lp',20,'fms',0,'ms',0);
plotbra('dsnake','cmp',11,'cl',p2pc('g1'),'fplab',[4 12 20],'wnr',f); ylabel('||u||_2');
%% soln plots
plotsol('wall','pt20'); solstyle('20'); pause; plotsol('dsnake','fpt4'); solstyle('FP4'); pause
plotsol('dsnake','fpt12'); solstyle('FP12'); pause; plotsol('dsnake','fpt20'); solstyle('FP20');
%% -------- FP-continuation to illustrate daisy-snake breakup at nu=1.5 ---------------
figure(2); clf; p=spcontini('dsnake','fpt12',2,'f12'); p.sol.ds=0.01; p.sol.dsmax=0.01;
p.plot.bpcmp=1; p.sw.bifcheck=0; p.usrlam=1.5:0.1:2; p.file.smod=1; p.sw.para=2;
p.fuha.spjac=@spjac; p.fuha.lss=@lss; p.fuha.blss=@blss; p=cont(p,20);
%% FP-cont plots
f=7; c=1; figure(f); clf; plotbra('f12','pt20',f,c,'cl','m', 'lab',[12]);
xlabel('\nu'); ylabel('\epsilon'); plotsol('f12','pt12'); nolti;
%% return to eps-cont
p=spcontexit('f12','pt12','w1b'); p.sol.ds=1e-2; p.nc.tol=1e-4; clf(2); p.sw.spcalc=1;
p.plot.bpcmp=11; p.file.smod=5; p.sw.foldcheck=0; p=cont(p,100);
%% other direction
p=loadp('w1b','pt5','w1br'); p.sol.ds=-p.sol.ds; p=cont(p,100);
%% branch plots
f=11; c=11; figure(f); clf; plotbra('w1b','pt100',f,c,'cl',p2pc('g1'), 'lab',50);
plotbra('w1br','pt100',f,c,'cl',p2pc('r1'), 'lab',100);
xlabel('\epsilon'); ylabel('||u||_2');
%% soln plots
plotsol('w1b','pt50'); solstyle('50'); pause; plotsol('w1br','pt100'); solstyle('100');
%% ---------- Solutions on full disk: create structure for full disk, then mirror
% half-disk soln to full disk, switch on pertinent phase conditions, then cont
p0=[]; dsw=1; rad=14; nref=5; p0=shinit(p0,dsw,par,rad,nref);
%% daisies, need only rotational PC
p=h2fdisk('dsnake','pt30',p0,'dsnakefull');
p=resetc(p); p.sw.bifcheck=0; p=qroton(p); p.sol.ds=-p.sol.ds; p=cont(p,30);
%% radial, needs x and y phase-conditions:
p=h2fdisk('r','pt10',p0,'rfull');
p=resetc(p); p.sw.bifcheck=0; p=qxyon(p); p.sol.ds=-0.02; p.nc.dsmax=0.05;
p.nc.eigref=-1; p=cont(p,20);
%% D_4^-, needs all 3 phase-conditions:
p=h2fdisk('d4m','pt20',p0,'d4mfull');
p=resetc(p); p.sw.bifcheck=0; p=qxyron(p); p.sol.ds=-0.02; p.nc.dsmax=0.05;
p.nc.eigref=-1; p=cont(p,20);
%% compare half-disk and full-disk branches
figure(3), clf;
plotbra('dsnake','cmp',11); plotbra('dsnakefull','cmp',11,'cl','m'); ylabel('||u||_2'); JBDMV-matlab/sh35disk/qxder.m 0000664 0001750 0001750 00000000211 14040302166 013716 0 ustar hu hu function qu=qxder(p,u) % derivative of standard transl. phase condition in x
qx=(p.mat.Dx(1:p.np,1:p.np)*p.u(1:p.np))'; qu=[qx, 0*qx]; JBDMV-matlab/sh35disk/shbra.m 0000755 0001750 0001750 00000000743 14035423621 013711 0 ustar hu hu function out=shbra(p,u)
M=getM(p); n=p.np; J=shJ(p,u); po=getpte(p);
v=p.u(1:p.nu/2); radr=find(abs(po(2,:))<1e-2); xy=[po(1,radr); v(radr)'];
xy=xy'; xy=sortrows(xy,1); xmax=max(xy(:,1));
%figure(10); clf; plot(xy(:,1),xy(:,2),'-')
out=[u(p.nu+1:end);
sqrt(trapz(xy(:,1),xy(:,2).^2))/sqrt(xmax);
max(v);
min(v);
J; sqrt(u(1:n)'*(M(1:n,1:n)*u(1:n)))/sqrt(p.Om)];
%sqrt(trapz(xy(:,1),xy(:,2).^2/30))];
% out=[max(v);sqrt(trapz(xy(:,1),xy(:,2).^2/30))]; JBDMV-matlab/sh35disk/qxyder.m 0000664 0001750 0001750 00000000266 14040302237 014120 0 ustar hu hu function qu=qxyder(p,u) % derivative of all x-y phase conditions
np=p.np; uold=p.u(1:np); u0x=p.mat.Dx(1:np,1:np)*uold; u0y=p.mat.Dy(1:np,1:np)*uold;
qu=[u0x' 0*u0x';u0y' 0*u0y']; JBDMV-matlab/sh35disk/oosetfemops.m 0000755 0001750 0001750 00000001167 14040305077 015156 0 ustar hu hu function p=oosetfemops(p) % for SH as 2nd order system, hence singular p.mat.M
fem=p.pdeo.fem; gr=p.pdeo.grid; [K,M]=assem6(p); % using 6-node triangles for Lapl.and mass
p.mat.M=[[M 0*M];[0*M 0*M]]; % system mass matrix (here singular)
p.mat.Ks=K; p.mat.Ms=M; % save SCALAR Laplacian, system K composed in sG
x=gr.p(1,:); y=gr.p(2,:); % for setting up PCs (in 3-node discretization)
[Dx,Dy]=fem.gradientMatrices(gr); E=center2PointMatrix(gr); Dx=E*Dx; Dy=E*Dy;
n=p.np; Krot=spdiags(-y',0,n,n)*Dx+spdiags(x',0,n,n)*Dy;
p.mat.Dx=[Dx 0*Dx; 0*Dx 0*Dx]; p.mat.Dy=[Dy 0*Dy; 0*Dy 0*Dy]; p.mat.Krot=[Krot 0*Krot; 0*Krot 0*Krot]; JBDMV-matlab/sh35disk/spjac.m 0000664 0001750 0001750 00000000675 14035423621 013715 0 ustar hu hu function Guuph=spjac(p,u)
u1=u(1:p.np); u2=u(p.np+1:2*p.np);
par=u(2*p.nu+1:end); nup=par(2);
n=p.np; ov=ones(n,1);
f1uu=6*nup*u1-20*u1.^3; f1uv=0*u1;
f1vv=0*u1; f2uu=0*u1; f2uv=f1uv; f2vv=f1vv;
ph1=u(p.nu+1:p.nu+p.np); ph2=u(p.nu+p.np+1:2*p.nu);
M1=spdiags(f1uu.*ph1+f1uv.*ph2,0,n,n);
M2=spdiags(f1uv.*ph1+f1vv.*ph2,0,n,n);
M3=spdiags(f2uu.*ph1+f2uv.*ph2,0,n,n);
M4=spdiags(f2uv.*ph1+f2vv.*ph2,0,n,n);
Guuph=-p.mat.M*[[M1 M2]; [M3 M4]]; JBDMV-matlab/sh35disk/qx.m 0000664 0001750 0001750 00000000215 14040302156 013226 0 ustar hu hu function q=qx(p,u) % standard phase condition for translational invariance in 'x'
u0x=p.mat.Dx(1:p.np,1:p.np)*p.u(1:p.np); q=u0x'*u(1:p.np); JBDMV-matlab/sh35disk/qrotder.m 0000664 0001750 0001750 00000000173 14040312163 014260 0 ustar hu hu function qu=qrotder(p,u) % derivative of rotational phase-cond
n=p.np; qr=(p.mat.Krot(1:n,1:n)*p.u(1:n))'; qu=[qr, 0*qr]; JBDMV-matlab/sh35disk/qroton.m 0000664 0001750 0001750 00000000200 14040302705 014112 0 ustar hu hu function p=qroton(p) % switch on rotational phase-condition
p.nc.nq=1; p.fuha.qf=@qrot; p.fuha.qfder=@qrotder; p.nc.ilam=[1 6]; JBDMV-matlab/sh35disk/solstyle.m 0000664 0001750 0001750 00000000472 14040235717 014471 0 ustar hu hu function solstyle(txt,varargin) % modify soln plots
nolti; colormap parula; title(''); colorbar off; h=colorbar('southoutside');
%h.Ticks=[-0.5 0.5]; % for low ampl.soln plots, adapt ticks of colorbar by hand
if nargin==1; text(9,12,txt,'fontsize',16);
else text(9,12,txt,'fontsize',16,'color',varargin{1});
end JBDMV-matlab/sh35disk/qy.m 0000664 0001750 0001750 00000000167 14040320064 013233 0 ustar hu hu function q=qy(p,u) % phase condition for transl.invariance in y
n=p.np; u0y=p.mat.Dy(1:n,1:n)*p.u(1:n); q=u0y'*u(1:n); JBDMV-matlab/sh35disk/h2fdisk0.m 0000664 0001750 0001750 00000001232 14040301401 014200 0 ustar hu hu function p=h2fdisk0(dir,pt,p,newdir) % extend half-disk soln to full disk by 0.
% Input p must contain mesh etc for full disk
q=loadp(dir,pt,newdir); par=q.u(q.nu+1:end)'; % half-disk-soln
np=p.np; hp=q.pdeo.grid.p; nph=q.np; % half points
for i=1:np;
mz=0;
po=p.pdeo.grid.p(:,i); x=po(1); y=po(2);
if x<0; x=-x; mz=1; end
dx=hp-repmat([x;y],1,nph); dxv=vecnorm(dx);
[mi,in]=min(dxv); % find point in p closest to point in q
if mz==1; p.u(i)=0; p.u(np+i)=0;
else p.u(i)=q.u(in); p.u(np+i)=q.u(nph+in);
end
end
p.u(p.nu+1:end)=par; p.file=q.file; p.nc.dsmax=q.nc.dsmax; % cp stuff from half disk
p=setfn(p,newdir); JBDMV-matlab/sh35disk/qxy.m 0000664 0001750 0001750 00000000246 14040302251 013417 0 ustar hu hu function q=qxy(p,u) % x and y phase conditions
np=p.np; u=u(1:np); uold=p.u(1:np); u0x=p.mat.Dx(1:np,1:np)*uold;
u0y=p.mat.Dy(1:np,1:np)*uold; q=[u0x'*u; u0y'*u];
JBDMV-matlab/sh35disk/qxyr.m 0000664 0001750 0001750 00000000343 14040302276 013606 0 ustar hu hu function q=qxyr(p,u) % 3 phase conditions: x,y, and rotational
np=p.np; u=u(1:np); uold=p.u(1:np); u0x=p.mat.Dx(1:np,1:np)*uold;
u0y=p.mat.Dy(1:np,1:np)*uold; u0phi=p.mat.Krot(1:np,1:np)*uold;
q=[u0x'*u; u0y'*u; u0phi'*u];
JBDMV-matlab/sh35disk/shinit.m 0000755 0001750 0001750 00000004044 14040322076 014104 0 ustar hu hu function p=shinit(p,domsw,par,rad,nref,varargin) % Swift-Hohenberg as 2 component system
p=stanparam(p); p.nc.neq=2; p.nc.neig=4; % # evals to compute (few, for speed)
p.fuha.outfu=@shbra; p.sol.ds=0.01; p.sol.dsmax=0.1;
p.sw.bifcheck=2; p.sw.foldcheck=1; p.sw.spcalc=1; p.sw.sfem=-1;
pde=stanpdeo2D(1,1,10); % dummy pdeo, next filled by 6-node triangulations
switch domsw
case 1; % disk
[p.nt,p.np,p.po,tri,p.efl,p.gfl]=trgl6_disk(nref,rad); % 6-node-triangulation
c3=c3fu(tri,p.nt); size(c3), c3=[c3, ones(size(c3,1),1)]; % convert to 3-node triangle-mesh
pde.grid.p=p.po'; pde.grid.t=c3'; p.pdeo=pde; % store 3-node-mesh (for plotting and Krot)
p.hofem.tri=tri;
case 2; % half disk
[p.nt,p.np,p.po,tri,p.efl,p.gfl]=trgl6_hdisk(nref,rad); % 6-node-triangulation
c3=c3fu(tri,p.nt); size(c3), c3=[c3, ones(size(c3,1),1)]; % convert to 3-node triangle-mesh
pde.grid.p=p.po'; pde.grid.t=c3'; p.pdeo=pde; % store 3-node-mesh (for plotting and Krot)
p.hofem.tri=tri;
case 3; % quarter disk
[p.nt,p.np,p.po,tri,p.efl,p.gfl]=trgl6_qdisk(nref,rad); % 6-node-triangulation
c3=c3fu(tri,p.nt); size(c3), c3=[c3, ones(size(c3,1),1)]; % convert to 3-node triangle-mesh
pde.grid.p=p.po'; pde.grid.t=c3'; p.pdeo=pde; % store 3-node-mesh (for plotting and Krot)
p.hofem.tri=tri;
case 4; % disk sector
aux=varargin{2}; p.pdeo=secpdeo2(rad,aux.nr,aux.nphi,aux.al,aux.phi,aux.jf);
p.np=p.pdeo.grid.nPoints; p.u=zeros(2*p.np,1); % prepare conversion to 6-nodes-triang.
p.t2sinterpol=0; % u re-initialized below, hence no need to interpolate u (or tau) to 6-nodes mesh,
p=tri2six(p.pdeo.grid.t',p.pdeo.grid.p',p); % convert to 6-node mesh
end
p.plot.pstyle=2; p.plot.cm='parula'; p.plot.axis='image'; p.plot.bpcmp=11;
p.nu=p.np*p.nc.neq; p.sol.xi=1/p.nu; p.nc.lammin=-4; p.nc.lammax=2;
u=0*ones(p.np,1); v=u; u0=[u v]; p.u=u0(:);
p.u=[p.u; par']; p.nc.ilam=1; p.file.smod=5; p=setfemops(p); huclean(p);
p.plot.auxdict={'\epsilon','\nu','q','sx','sy','srot'};
p.Om=sum(sum(p.mat.M(1:p.np,1:p.np),1)); % domain size, for normalization of L2-norm JBDMV-matlab/sh35disk/h2fdisk.m 0000664 0001750 0001750 00000001075 14040301250 014127 0 ustar hu hu function p=h2fdisk(dir,pt,p,newdir) % mirror half-disk to full
% input p must contain mesh etc for full disk
q=loadp(dir,pt,newdir); par=q.u(q.nu+1:end)'; % half-disk-soln
np=p.np; hp=q.pdeo.grid.p; nph=q.np; % half points
for i=1:np;
po=p.pdeo.grid.p(:,i); x=po(1); y=po(2);
if x<0; x=-x; end
dx=hp-repmat([x;y],1,nph); dxv=vecnorm(dx);
[mi,in]=min(dxv); % find point in p closest to point in q
p.u(i)=q.u(in); p.u(np+i)=q.u(nph+in);
end
p.u(p.nu+1:end)=par; p.file=q.file; p.nc.dsmax=q.nc.dsmax; % cp stuff from half disk
p=setfn(p,newdir);