Implements an explicit Runge-Kutta 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 b, column vector holding the Runge-Kutta weights c, column vector holding the Runge-Kutta nodes A, matrix holding the coefficients a(i,j) which prescribe the method. 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 ] = explicit_RK( f,t0,y0,h,T,b,c,A,analytic ) 0002 % Implements an explicit Runge-Kutta 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 % b, column vector holding the Runge-Kutta weights 0010 % c, column vector holding the Runge-Kutta nodes 0011 % A, matrix holding the coefficients a(i,j) which prescribe the method. 0012 % analytic (optional), function handle to the analytic solution; 0013 % analytic=@(t)() 0014 % Output arguments: 0015 % y, vector of solution at times defined by t 0016 % t, vector of time steps 0017 % err, pointwise error. If the analytic solution is not provided, err will be empty 0018 0019 % first check user inputs 0020 [Arows,Acols]=size(A); 0021 [brows,bcols]=size(b); 0022 [crows,ccols]=size(c); 0023 if Arows~=Acols; 0024 error('matrix must be a square matrix'); 0025 elseif ~istril(A) 0026 error(['matrix must be strictly lower triangular, else the prescribed '... 0027 'method is not explicit']); 0028 elseif norm(diag(A))~=0; 0029 error(['matrix must be strictly lower triangular, else the prescribed '... 0030 'method is not explicit']); 0031 elseif bcols~=1 || ccols~=1; 0032 error('input must be column vectors'); 0033 elseif brows~=Arows || crows~=Arows; 0034 error('dimensions do not match'); 0035 elseif ~isequal(c, sum(A,2)); 0036 error('the method is not consistent'); 0037 elseif isa(f,'function_handle')==0 0038 error('input must be a function handle'); 0039 elseif h<=0; 0040 error('step size must be positive'); 0041 elseif T<=t0; 0042 error('invalid interval'); 0043 end 0044 0045 n=ceil((T-t0)/h); % calculate the ceiling of (T-t0)/h. 0046 % This gives the number of steps 0047 t=linspace(t0,t0+h*n,n+1); t=t'; % initialize vector of time steps 0048 y=zeros(n+1,1); % initialize vector to hold y values 0049 y(1)=y0; % set initial y value 0050 0051 k=zeros(Arows,1); % initialize vector k 0052 0053 for i=1:n; 0054 k(1) = feval(f,t(i),y(i)); % because the method is explicit and consistent 0055 for j=2:Arows; 0056 s=0; % initialize sum 0057 for l=1:j-1; 0058 s=s+A(j,l)*k(l); % calculate sum 0059 end 0060 k(j)=feval(f,t(i)+c(j)*h,y(i)+h*s); % calculate k(j) 0061 end 0062 y(i+1)=y(i)+h*(b'*k); % calculate y(i+1) 0063 end 0064 0065 if exist('analytic')==1; 0066 if isa(analytic,'function_handle')==0; 0067 disp('analytic must be a function handle') 0068 return; 0069 else 0070 true=analytic(t); 0071 err=abs(true-y); 0072 end 0073 end 0074 0075 end