%near_bifurcation.m %Code used to generate data in figure 10 in Shemilt et al. (2022). %Solves the primal and adjoint problems defined in Appendix C using bvp4c %BVP solver, then calculates the coefficients in (C11) then solves (C11) %using ode45. %Calls the script TF_bifurdiagram.m %Written by James D Shemilt, 02/2022. Bstar = .16296; %bifurcation point L = pi/k; %domain length C = 1; %set value of phi_1(L) (this is arbitrary - final solution does not depend on C) %Run TF_bifurdiagram to find static solution at B=Bstar, bifurcation. TF_bifurdiagram %--% Solve the phi_1 problem using the H_0 computed above %--% %Construct BVP xmesh = linspace(0,L,150); solinit = bvpinit(xmesh, @(x)guess(x,C,L)); options = bvpset('RelTol',1e-9,'AbsTol',1e-9); sol = bvp4c(@(x,y)bvpfn(x,y,sol_Bbif,Bstar), @(ya,yb)bcfn(ya,yb,C), solinit,options); %Define phi_1 phi_1 = sol.y(1,:); %Print error in volume of solution V_err = trapz(sol.x,sol.y(1,:)) %Integrate to determine phi_2 phi_2 = cumtrapz(sol.x,sol.y(1,:)); %--% Solve the Adjoint problem for phi^dagger %--% %Construct BVP xmesh = linspace(0,L,150); solinitAdj = bvpinit(xmesh, @(x) guessAdj(x,L)); optionsAdj = bvpset('RelTol',1e-9,'AbsTol',1e-9); solAdj = bvp4c(@(x,y)bvpfnAdj(x,y,sol_Bbif,Bstar), @bcfnAdj, solinit,options); %Define phi_1_dag phi_1_dag = solAdj.y(1,:); %--% Evaluate integrals for coefficients in the A(T) evolution equation %-% % Spline everything onto the adjoint solution mesh % phi_1_spl = spline(sol.x,phi_1,solAdj.x); phi_2_spl = spline(sol.x,phi_2,solAdj.x); H0_spl = spline(sol_Bbif.x,sol_Bbif.y(1,:),solAdj.x); %Determine the coefficients of the evolution equation for A(T), (C11). a = real(trapz(solAdj.x,phi_1_dag*(sqrt(2*Bstar)).*sqrt(-phi_2_spl)./(H0_spl.^2))); b = real(trapz(solAdj.x,-phi_1_dag./H0_spl)); c = real(trapz(solAdj.x,phi_1_dag.*(phi_1_spl.^2).*Bstar./(H0_spl.^3))); %--% Solve (C11) for A(T) %--% T = 200; %end-time of integration A0 = -2.19; %initial condition options = odeset('RelTol',1e-8); %Set error tolerance solA = ode45(@(T,A) ((b+c*A^2)/a)^2, [0 T], A0,options);%solve ODE %--% Functions for Primal Problem %--% function dydx = bvpfn(x,y,solgu2,Bstar) H0 = deval(solgu2,x); H0 = H0(1); dydx = [y(2) y(3) -y(2)-Bstar*y(1)/(H0^2)]; end function res = bcfn(ya,yb,C) res = [ya(2) yb(1)-C yb(2) ]; end function g = guess(x,C,L) g = zeros(3,1); g(1) = 0.5*C*(1-cos(pi*x/L)); g(2) = 0.5*C*pi*sin(pi*x/L)/L; g(3) = 0.5*C*pi*pi*cos(pi*x/L)/(L^2); end %--% Functions for Adjoint Problem %--% function dydx = bvpfnAdj(x,y,solgu2,Bstar) H0 = deval(solgu2,x); H0 = H0(1); dydx = [y(2) y(3) -y(2)+Bstar*y(1)/(H0^2)+1]; end function res = bcfnAdj(ya,yb) res = [ya(1) yb(3) yb(1) ]; end function g = guessAdj(x,L) a = L/pi; g = zeros(3,1); g(1) = (a/(1-(1/(a^2))))*sin(x/a); g(2) = (1/(1-(1/(a^2))))*cos(x/a); g(3) = -(1/(a*(1-(1/(a^2)))))*cos(x/a); end