%latetimefig3.m %Code used to generate the data in figure 3 of Shemilt et al. (2022) %Calls function thinfilm_evolution to find late time numerical solution, %then solve the leading-order (static) problem for H_0, then solves the %O(1/t) problem to find H_1 and Y_1. Uses bvp4c to solve BVPs. %Written by James D Shemilt, 02/2022 B = 0.05; L = sqrt(2)*pi; A = 0.4; N = 81; dz = L/(N-1); X = 0:dz:L; Tend = 1e5; %get long time solution from simulation [sol, H, Y, pz, tauw] = thinfilm_evolution(A,B,1e-6,Tend,N,L,1e-11,[0 Tend]); %Approximate H_0 and its derivatives, Y, H1 and its derivatives by finite %differences. These will be used as guess functions for the BVP solvers. H = sol.y(:,end); Hz = zeros(N,1); Hzz = zeros(N,1); Hzzz = zeros(N,1); Y = zeros(N,1); H1 = sol.y(:,end); H1z = zeros(N,1); H1zz = zeros(N,1); Hz(2:N-1) = (H(3:N)-H(1:N-2))./(2*dz); Hzz(1) = 2*(H(2)-H(1))/(dz^2); Hzz(2:N-1) = (H(3:N)-2.*H(2:N-1)+H(1:N-2))./(dz^2); Hzz(N) = 2.*(H(N-1)-H(N))/(dz^2); Hzzz(2) = (-0.5*H(2)+ H(1)- H(3)+ 0.5*H(4))/(dz^3); Hzzz(N-1) = (-0.5*H(N-3)+ H(N-2)- H(N)+ 0.5 *H(N-1))/(dz^3); Hzzz(3:N-2) = (-0.5.*H(1:N-4) +H(2:N-3) -H(4:N-1) + 0.5 .*H(5:N))./(dz^3); Y(2:N-1) = B.*Tend.*(H(2:N-1) - B./abs(Hz(2:N-1)+Hzzz(2:N-1))); Y2 = Y.^2; H1(2:N-1) = (1/2).*(Y2(3:N)-Y2(1:N-2))./(2*dz); H1(1) = H1(2); H1(N) = H1(N-1); H1z(2:N-1) = (H1(3:N)-H1(1:N-2))./(2*dz); H1zz(2:N-1) = (H1z(3:N)-H1z(1:N-2))./(2*dz); %--% Leading order solution %--% %Calculate H(z=L) Cf = fzero(@(c)fun(c,L,X,H,Hz,Hzz,B),H(end)); %Construct BVP for leading order xmesh = linspace(0,L,150); solinit0 = bvpinit(xmesh, @(x)guess(x,X,H,Hz,Hzz)); options = bvpset('RelTol',1e-9,'AbsTol',1e-9); sol0 = bvp4c(@(x,y)bvpfn(x,y,B), @(ya,yb)bcfn(ya,yb,Cf), solinit0,options);%leading order soln %--% O(1/t) problem %--% %Approximate H_1 for the guess for bvp4c H1 = (H' - H0).*B.*Tend; %Construct BVProblem solinit1 = bvpinit(xmesh,@(x) guess1(x,X,H1,H1z,H1zz,Y2)); sol1 = bvp4c(@(x,y)bvpfn1(x,y,sol0,B),@bcfn1,solinit1,options);%next order soln %--% functions %--% %function to shoot for fixed volume V = L (leading-order problem) function Verror = fun(c,L,X,H,Hz,Hzz,B) %Construct BVP xmesh = linspace(0,L,150); solinit = bvpinit(xmesh, @(x)guess(x,X,H,Hz,Hzz)); options = bvpset('RelTol',1e-9,'AbsTol',1e-9); sol = bvp4c(@(x,y)bvpfn(x,y,B), @(ya,yb)bcfn(ya,yb,c), solinit,options); %Integrate for volume V = trapz(sol.x,sol.y(1,:)); Verror = V - L; end %--% Functions to solve the leading-order problem %--% function dydx = bvpfn(x,y,B) %y(1)-y(3) are H_0 and its derivatives dydx = [y(2) y(3) (B/y(1))-y(2)]; end function res = bcfn(ya,yb,C) res = [ya(2) yb(1)-C yb(2) ]; end function g = guess(z,X,H,Hz,Hzz) %Spline the finite differenced approximation to H_0 and its derivatives %from the numerical solution. g = zeros(3,1); g(1) = spline(X,H,z); g(2) = spline(X,Hz,z); g(3) = spline(X,Hzz,z); end %--% Functions to solve the O(1/t) problem %--% function dydx = bvpfn1(x,y,sol,B) %y(1)-y(3) are H_1 and derivatives %y(4) is Y_1^2 Z = deval(sol,x); dydx = [y(2) y(3) -y(2)-(B/(Z(1)^2))*(y(1)-sqrt(y(4))) 2*y(1)]; end function res = bcfn1(ya,yb) res = [ya(2) yb(2) ya(4) yb(4)]; end function g = guess1(z,X,H1,H1z,H1zz,Y2) %Spline the finite differenced approximation to H1 and its derivatives and %Y_1^2 calculated from the numerical solution. g = zeros(4,1); g(1) = spline(X,H1,z); g(2) = spline(X,H1z,z); g(3) = spline(X,H1zz,z); g(4) = spline(X,Y2,z); end