Home > K25104 > ODEs > trapezoidal_rule.m

trapezoidal_rule

PURPOSE ^

Implements the trapezoidal rule to solve y'=f.

SYNOPSIS ^

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

DESCRIPTION ^

 Implements the trapezoidal rule 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
   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 ] = trapezoidal_rule( f,t0,y0,h,T,analytic )
0002 % Implements the trapezoidal rule 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 %   analytic (optional), function handle to the analytic solution;
0010 %     analytic=@(t)()
0011 % Output arguments:
0012 %   y, vector of solution at times defined by t
0013 %   t, vector of time steps
0014 %   err, pointwise error. If the analytic solution is not provided,
0015 %     err will be empty
0016 
0017 % check user input
0018 if isa(f,'function_handle')==0
0019     error('input must be a function handle');
0020 elseif h<=0;    
0021     error('step size must be positive');
0022 elseif T<=t0;
0023     error('invalid interval');
0024 end
0025 
0026 n=ceil((T-t0)/h);   % calculate the ceiling of (T-t0)/h.
0027                     % This gives the number of steps
0028 t=linspace(t0,t0+h*n,n+1); t=t';    % initialize vector of time steps
0029 y=zeros(n+1,1);     % initialize vector to hold y values
0030 y(1)=y0;            % set initial y value
0031 
0032 for k=2:n+1;        % move forward in time with trapezoidal rule
0033     % find next iterate as zero of a non-linear function
0034     g = @(x) (x - y(k-1) - h/2 * (f(t(k-1),y(k-1)) + f(t(k),x)));
0035     y(k)=fzero(g,y(k-1));
0036 end
0037 
0038 % calculate error if possible
0039 err=[]; % initialize err
0040 if exist('analytic')==1;
0041     if isa(analytic,'function_handle')==0;
0042         disp('analytic must be a function handle')
0043         return;
0044     else
0045         true=analytic(t);
0046         err=abs(true-y);
0047     end
0048 end
0049 
0050 end
0051

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