Home > K25104 > ODEs > explicit_RK.m

explicit_RK

PURPOSE ^

Implements an explicit Runge-Kutta method to solve y'=f.

SYNOPSIS ^

function [ y,t,err ] = explicit_RK( f,t0,y0,h,T,b,c,A,analytic )

DESCRIPTION ^

 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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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

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