%%This section calculates the lambda and nu values at each of the four %%branch points in the generalised absolute spectrum. %%Copyright (c) 2008, Microsoft %%All rights reserved. %%Copyright details (BSD License) is included at the end of this code %%Written by Matthew Smith, Microsoft Research Ltd. for the paper %%"Explaining spatiotemporal dynamics in lambda-omega %%reaction-diffusion equations using absolute stability spectra of %%wavetrains" %%by Matthew J. Smith, Jens D.M. Rademacher and Jonathan A. Sherratt. SUBMITTED MANUSCRIPT. clear all w0=3.0d0; %Constant in lambda-omega equations w1=1.0d0; %Constant in lambda-omega equations syms n f r fd l gd rh th %symbolic parameters D=n^4+n^2*(4*f+r*fd-2*l)+2*n*r*f^(1/2)*gd+l^2-r*fd*l; %This is the general dispersion relation for the lambda omega equations D2=diff(D,n); %differentiate it s1=solve(D2,l); %solve it to get one value of lambda D3=subs(D,l,s1); %substitute it back into D to get a 4th power polynomial equation that can be solved for nu, this gives the values associated with the branch points amp=1./sqrt(0.5*(1+sqrt(1+(8/9)*w1.^2))); %amplitude of wave selected by dirichlet boundary %amp=sqrt((2/w1^2)*(sqrt(1+w1^2)-1)); %amplitude of wave selected by invasion fr=(1-amp^2); %our "lambda" function (1-R^2) fr2=-2*amp; %the derivative of the "lambda" function gdr=(-2*w1*amp); %The derivative of the omega function D4=subs(D3,{r,f,fd,gd},{amp,fr,fr2,gdr}); %sub the values into our polynomial for the branch points s2=double(solve(D4,n)); %solve the polynomial to get the nu values for the branch points s3=subs(s1,{r,f,fd,gd},{amp,fr,fr2,gdr}); %substitute these nu values in to get corresponding lambda values - step 1 s4=subs(s3,n,s2); %gives the 4 lambda values at the branch points lambdafile=zeros(4,2); %variable to save the 4 lambdas nufile=zeros(4,4,2); %variable to save the double root nu values associated with each lambda uvfile=zeros(4,4,4); %variable to save the 4 sets of 4 u and v values associated with each nu lambdafile(:,:)=[real(s4) imag(s4)]; for j=1:4 %for each lambda s5=subs(D,{r,f,fd,gd,l},{amp,fr,fr2,gdr,s4(j)}); %sub lambda back in to equation s6=single(solve(s5)); %%calculate the 4 nus at one branch point nufile(j,:,:)=[real(s6) imag(s6)]; %save the nu values for k=1:4 %for each nu within each lambda Dscript=[[s6(k)^2+amp*fr2-s4(j)],[-2*s6(k)*amp*fr^(1/2)];[gdr+2*s6(k)*fr^(1/2)/amp],[s6(k)^2-s4(j)]]; %write out the dispersion relation a=-Dscript(1,2); %work out the a coefficient b=Dscript(1,1); %work out the b coefficient rnew=a/sqrt(a^2+b^2); %work out one component of the eigenvector thnew=b/sqrt(a^2+b^2); %work out the other component of the eigenvector uvfile(j,k,:)=[real(rnew) real(thnew) imag(rnew) imag(thnew)]; %arrange the eigenvectors in the appropriate format end end %%This section reports absolute and essential stability stablim=sqrt((2+2*w1^2)/(3+2*w1^2)) if amp>stablim EssentialStability='stable' else EssentialStability='unstable' end AbsoluteStability='stable'; tol = 0.000001; for j=1:4 if lambdafile(j,1)>0 renus=nufile(j,:,1); renus=sortrows(renus'); if abs(renus(2)-renus(3))