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
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