Home > K25104 > ODEs > RK_order.m

RK_order

PURPOSE ^

Identifies the order up to 4 of the Runge-Kutta scheme given by A,b,c.

SYNOPSIS ^

function [ p ] = RK_order( b,c,A )

DESCRIPTION ^

 Identifies the order up to 4 of the Runge-Kutta scheme given by A,b,c.
 Input arguments:
   b, column vector holding the Runge-Kutta weights
   c, column vector holding the Runge-Kutta nodes
   A, matrix holding the Runge-Kutta coefficients
 Output arguments:
   p, the order of the scheme up to 4.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [ p ] = RK_order( b,c,A )
0002 % Identifies the order up to 4 of the Runge-Kutta scheme given by A,b,c.
0003 % Input arguments:
0004 %   b, column vector holding the Runge-Kutta weights
0005 %   c, column vector holding the Runge-Kutta nodes
0006 %   A, matrix holding the Runge-Kutta coefficients
0007 % Output arguments:
0008 %   p, the order of the scheme up to 4.
0009 
0010 % check user inputs
0011 [Arows,Acols]=size(A);
0012 [brows,bcols]=size(b);
0013 [crows,ccols]=size(c);
0014 if Arows~=Acols;    
0015     error('input must be a square matrix');
0016 elseif bcols~=1 || ccols~=1;
0017     error('input must be column vectors');
0018 elseif brows~=Arows || crows~=Arows;
0019     error('input dimensions do not agree');
0020 end
0021 
0022 E=1e-12;    % we use this as a small tolerance on either side of each
0023             % condition, in case of possible rounding errors
0024 p=0;           
0025 if (sum(b)>1-E) && (sum(b)<1+E);                    % condition 1
0026     p=1;                    % update p
0027     if (b'*c>1/2-E) && (b'*c<1/2+E);                % condition 2
0028         p=2;                % update p
0029         c2=c.^2;            % '.^' performs element-wise exponentiation
0030         if (b'*c2>1/3-E) && (b'*c2<1/3+E);          % condition 3
0031             if (b'*A*c>1/6-E) && (b'*A*c<1/6+E)     % condition 4
0032                 p=3;        % update p
0033             end
0034         else return; 
0035         end
0036     else return; 
0037     end
0038 else return; 
0039 end
0040 
0041 c3=c.^3;
0042 if (b'*c3>1/4-E) & (b'*c3<1/4+E);                       % condition 5
0043     bc=b.*c;                % '.*' performs element-wise multiplication
0044     if (bc'*A*c>1/8-E) && (bc'*A*c<1/8+E);              % condition 6
0045         if (b'*A*c2>1/12-E) && (b'*A*c2<1/12+E);        % condition 7
0046             if (b'*A*A*c>1/24-E) && (b'*A*A*c<1/24+E)   % condition 8
0047                 p=4;        % update p
0048             end
0049         end
0050     end
0051 end
0052 
0053 
0054 end
0055

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