Implements an implicit multistep method to solve y'=f. If the analytic solution is given, the error is calculated. Input arguments: f, function handle; f=@(t,y)() t0, y0 initial conditions h, step size T, end time rho, sigma, column vectors prescribing the method. The first element should contain the coefficient associated with l=0. analytic (optional), function handle to the analytic solution; analytic=@(t)() Output arguments: y, vector of solution at times defined by t t, vector of time steps err, pointwise error. If the analytic solution is not provided, err will be empty
0001 function [ y,t,err ] = implicit_multistep( f,t0,y0,h,T,rho, sigma, analytic ) 0002 % Implements an implicit multistep method to solve y'=f. 0003 % If the analytic solution is given, the error is calculated. 0004 % Input arguments: 0005 % f, function handle; f=@(t,y)() 0006 % t0, y0 initial conditions 0007 % h, step size 0008 % T, end time 0009 % rho, sigma, column vectors prescribing the method. 0010 % The first element should contain the coefficient associated with l=0. 0011 % analytic (optional), function handle to the analytic solution; 0012 % analytic=@(t)() 0013 % Output arguments: 0014 % y, vector of solution at times defined by t 0015 % t, vector of time steps 0016 % err, pointwise error. If the analytic solution is not provided, err will be empty 0017 0018 % check user input 0019 [n,m]=size(rho); % finding the size of rho 0020 [s,t]=size(sigma); % finding the size of sigma 0021 if m~=1 || t~=1; 0022 error('inputs must be column vectors'); 0023 elseif n~=s; 0024 error('inputs must be the same length'); 0025 % the method is an (s-1)-step method 0026 elseif isa(f,'function_handle')==0 0027 error('input must be a function handle'); 0028 elseif h<=0; 0029 error('step size must be positive'); 0030 elseif T<=t0; 0031 error('invalid interval'); 0032 end 0033 0034 s1 = s-1; % number of steps in the method 0035 n=ceil((T-t0)/h); % calculate the ceiling of (T-t0)/h. 0036 % This gives the number of steps 0037 t=linspace(t0,t0+h*n,n+1); t=t'; % initialize vector of time steps 0038 y=zeros(n+1,1); % initialize vector to hold y values 0039 fy=zeros(s1,1); % holds function values at s time steps 0040 0041 % Take first s-1 steps with a one step method. 0042 % This is not necessary if it is a one step method. 0043 if s>1; 0044 [a,b,~]=backward_euler(f,t0,y0,h,t0+s1*h,analytic); 0045 y(1:s1)=a(1:s1); 0046 for i=1:s1; 0047 fy(i)=f(b(i),a(i)); 0048 end 0049 else 0050 y(1)=y0; 0051 fy(1)=f(t0,y0); 0052 end 0053 0054 for k=0:n-s; 0055 % find next iterate as zero of a non-linear function 0056 g = @(x) rho(s)*x + (rho(1:s1)'*y(k+1:s1+k)) - ... 0057 h*(sigma(1:s1)'*fy) - h*sigma(s)*f(t(s+k),x); 0058 y(s+k)= fzero(g,y(s+k-1)); 0059 fy(1:s1-1)=fy(2:s1); 0060 fy(s1)=f(t(s+k),y(s+k)); 0061 end 0062 0063 % calculate error if possible 0064 err=[]; % initialize err 0065 if exist('analytic')==1; 0066 if isa(analytic,'function_handle')==0; 0067 disp('analytic must be a function handle') 0068 return; 0069 else 0070 true=analytic(t); 0071 err=abs(true-y); 0072 end 0073 end 0074 0075 end 0076