Home > K25104 > ODEs > forward_euler.m

forward_euler

PURPOSE ^

Implements forward Euler method to solve y'=f.

SYNOPSIS ^

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

DESCRIPTION ^

 Implements forward Euler 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
   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 ] = forward_euler( f,t0,y0,h,T,analytic )
0002 % Implements forward Euler 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 %   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, err will be empty
0015 
0016 % check user input
0017 if isa(f,'function_handle')==0
0018     error('input must be a function handle');
0019 elseif h<=0;    
0020     error('step size must be positive');
0021 elseif T<=t0;
0022     error('invalid interval');
0023 end
0024 
0025 n=ceil((T-t0)/h);   % calculate the ceiling of (T-t0)/h.
0026                     % This gives the number of steps
0027 t=linspace(t0,t0+h*n,n+1); t=t';    % initialize vector of time steps
0028 y=zeros(n+1,1);     % initialize vector to hold y values
0029 y(1)=y0;            % set initial y value
0030 
0031 for k=2:n+1;        % move forward in time with Forward Euler scheme
0032     y(k)=y(k-1)+h*f(t(k-1),y(k-1));
0033 end
0034 
0035 % calculate error if possible
0036 err=[]; % initialize err
0037 if exist('analytic')==1;
0038     if isa(analytic,'function_handle')==0;
0039         disp('analytic must be a function handle')
0040         return;
0041     else
0042         true=analytic(t);
0043         err=abs(true-y);
0044     end
0045 end
0046 
0047 end
0048

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