Implements backward 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 ] = backward_euler( f,t0,y0,h,T,analytic ) 0002 % Implements backward 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, 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 backward Euler scheme 0033 % find next iterate as zero of a non-linear function 0034 g = @(x) (x - y(k-1) - h * 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