Contents

% Skiba example, continues cpdemo.m
% find paths from the uv0 initial states from cpdemo.m to FSM
js=10; je=30; jl=js-je+1; % alpha-range; now set the target and Psi to FSM:
s1=loadp('f2','pt11'); u1=s1.u(1:s1.nu); [Psi,muv,d,t1]=getPsi(s1);
a0l=length(alv0); tva=zeros(jl,opt.Nmax+1); % some prep. and fields to hold paths
uva=zeros(jl,n+1,opt.Nmax+1); opt.start=0; alva=[]; vva=[]; tavl=[]; sol=[];
alvin=[0.1 0.25 0.5 0.75 1]; % we run from uv0(j,:) to FSM with iscnat
tv=linspace(0,opt.t1,opt.nti); se=2; opt.tv=tv.^se./opt.t1^(se-1); doplot=1;
opt.msw=0; opt.Stats_step='off'; v=[50,8]; % switch off stats
for j=js:je;
   fprintf('j=%i, al=%g\n', j, alv0(j)); s0.u(1:n)=uv0(j,1:n,1)';
   [alv,vv,sol,udat]=iscnat(alvin,[],[],opt,fn);
   if alv(end)==1; Jd=vv0(j)-vv(end); fprintf('J1-J2=%g\n',Jd); % cont. success
      alva=[alva alv0(j)]; vva=[vva vv(end)]; tl=length(sol.x); % put vals in vector
      talv=[tavl tl]; tva(j,1:tl)=sol.x; uva(j,1:n,1:tl)=sol.y;
      if abs(Jd)<0.05; % (approximate) Skiba point(s) found
        if doplot==1; sol0=[]; alp=alv0(j); % plot the paths to FSC and FSM
          sol0.x=tv0(j,1:tlv0(j));sol0.y=squeeze(uv0(j,1:n,1:tlv0(j)));
          psol3Dm(s1,sol0,sol,1,1,[]); zlabel('P');
          psol3Dm(s1,sol0,sol,2,0,[]); zlabel('k'); pause
        end
      end
   end
end

getting Psi, done in 0.044582 sn/2=102, d=0, suggested T=32.7303
j=10, al=0.424976
al=0.1, flag=0 
al=0.25, flag=0
al=0.5, flag=0
al=0.75, flag=0
al=1, flag=0
J1-J2=0.136151
j=11, al=0.437192
al=0.1, flag=0 
al=0.25, flag=0
al=0.5, flag=0
al=0.75, flag=0
al=1, flag=0
J1-J2=0.0849527
j=12, al=0.449595
al=0.1, flag=0 
al=0.25, flag=0
al=0.5, flag=0
al=0.75, flag=0
...

plot value diagram

jep=je; figure(6); clf; plot(alv0(js:jep),vv0(js:jep),'-*b');hold on;
plot(alva,vva,'-*r');set(gca,'FontSize',s1.plot.fs); axis tight;
xlabel('\alpha','FontSize',s1.plot.fs); ylabel('J_{a}','FontSize',s1.plot.fs);

evaluate selected path from continuation

j=22;sol0=[];v=[25,15]; n=s1.nu;
tl=tlv0(j); sol0.x=tv0(j,1:tlv0(j));sol0.y=squeeze(uv0(j,1:n,1:tlv0(j)));
al=alv0(j), u0=al*s0.u(1:n)+(1-al)*s1.u(1:n); u1=s1.u(1:s1.nu); slsolplot(sol0,v);
psol3D(s1,sol,3,1,v,[]); zlabel('P'); psol3D(s1,sol,4,0,v,[]); zlabel('k');

al =
    0.6035

Published with MATLAB® 7.14