Home > K25104 > ODEs > explicit_multistep.m

explicit_multistep

PURPOSE ^

Implements an explicit multistep method to solve y'=f.

SYNOPSIS ^

function [ y,t,err ] = explicit_multistep( f,t0,y0,h,T,rho, sigma, analytic )

DESCRIPTION ^

 Implements an explicit 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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [ y,t,err ] = explicit_multistep( f,t0,y0,h,T,rho, sigma, analytic )
0002 % Implements an explicit 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 sigma(s)~=0;
0027     error('given method is not explicit');
0028 elseif rho(s)==0;
0029     error('something wrong with the number of steps');
0030 elseif isa(f,'function_handle')==0
0031     error('input must be a function handle');
0032 elseif h<=0;    
0033     error('step size must be positive');
0034 elseif T<=t0;
0035     error('invalid interval');
0036 end
0037 
0038 s1 = s-1; % number of steps in the method
0039 n=ceil((T-t0)/h);   % calculate the ceiling of (T-t0)/h.
0040                     % This gives the number of steps
0041 t=linspace(t0,t0+h*n,n+1); t=t';    % initialize vector of time steps
0042 y=zeros(n+1,1);     % initialize vector to hold y values
0043 fy=zeros(s1,1);      % holds function values at s time steps
0044 
0045 % Take first s-1 steps with a one step method.
0046 % This is not necessary if it is a one step method.
0047 if s>1;
0048     [a,b,~]=forward_euler(f,t0,y0,h,t0+s1*h,analytic);
0049     y(1:s1)=a(1:s1);
0050     for i=1:s1;
0051         fy(i)=f(b(i),a(i));
0052     end
0053 else
0054     y(1)=y0;
0055     fy(1)=f(t0,y0);
0056 end
0057 
0058 for k=0:n-s;
0059     % calculate next value according to prescribed method
0060     y(s+k)= (h*(sigma(1:s1)'*fy)-(rho(1:s1)'*y(k+1:s1+k)))/rho(s); 
0061     fy(1:s1-1)=fy(2:s1);
0062     fy(s1)=f(t(s+k),y(s+k));
0063 end
0064 
0065 % calculate error if possible
0066 err=[]; % initialize err
0067 if exist('analytic')==1;
0068     if isa(analytic,'function_handle')==0;
0069         disp('analytic must be a function handle')
0070         return;
0071     else
0072         true=analytic(t);
0073         err=abs(true-y);
0074     end
0075 end
0076 
0077 end
0078

Generated on Mon 18-Jan-2016 10:25:49 by m2html © 2005