Solving Matrix Differential Equations using 4th Order Runge-Kutta Method

Good day all,

I am trying to create a script to employ the 4th order Runge Kutta method to solve a matrix differential equation where: d{V}/dt = [F(V)], where V is a 2x1 vector and F is a 2x2 matrix.

Previously I have successfully used the code below to solve the differential equation dy/dt = y*t^2 - 1.1*y

How should I adapt this code so I can input vectors?

close all; clear all; clc;

% This programme employs the 4th Order RK method

% Function is defined in a separate file funct.m where:

% function f = funct(t,y)

% f = y*t^2 - 1.1*y;

% end

% Create a 1 x 6 matrix A to contain all values req for the RK method

% Row 1: values of t

% Row 2: numerical values of y

% Row 3 to 6: values of k1 to k4 respectively

A = zeros(1,6);

% Input value of h

h = input('Input value of h: ');

% Input initial value of y and insert into row 1 column 2 of A

y0 = input('Input initial value of y: ');

A(1,2) = y0;

% Input lower and upper limit of t

L = input('Input lower limit of t: ');

T = input('Input upper limit of t: ');

% Loop to insert values of t in column 1 of A in increments of h

A(1,1) = L;

n = 1;

while n < (T-L)/h + 1

A(n+1,1) = A(n,1) + h;

n = n+1;

end

% Loop to calc k1 to k4 in columns 3-6 of A and find y(i+1) in column 2

n = 1;

while n < (T-L)/h + 1

% k1 in column 3

t = A(n,1);

y = A(n,2);

A(n,3) = funct(t,y);

% k2 in column 4

t = A(n,1) + 0.5*h;

y = A(n,2) + 0.5*h*A(n,3);

A(n,4) = funct(t,y);

% k3 in column 5

y = A(n,2)+0.5*h*A(n,4);

A(n,5) = funct(t,y);

% k4 in column 6

t = A(n,1)+h;

y = A(n,2)+h*A(n,5);

A(n,6) = funct(t,y);

% Insert y(i+1) into column 2 of the next row

A(n+1,2) = A(n,2) + (A(n,3)+2*(A(n,4)+A(n,5))+A(n,6))*h/6;

n = n+1;

end

A % Display matrix A (if req) to check values

% Plot result with column 1 (t) as hor-axis and column 2 (y) as vert-axis

h_axis=A(:,1);

v_axis=A(:,2);

plot(h_axis,v_axis); grid on

Answer 1

[Joel]

Hi,

When you say Matrix Differential equations, I assume you mean a system of first order ODEs which can be represented in the Matrix format. The procedure largely remains the same as RK4 for a single ODE. In the case of 1 ODE, we usually calculate K1, K2, K3, K4 in one iteration. For a system of 2 ODEs, you will have to calculate (K1,L1), (K2,L2), (K3,L3) and (K4,L4) in one iteration. You should also specify two initial conditions say Y(xo) and Z(xo).

This is the general algorithm:

Say,

Y’ = f1(x,y,z) and Z’ = f2(x,y,z)

Y(xo) and Z(xo) are defined.

h is step size

K1 = h*f1(xn,yn,zn)
L1 = h*f2(xn,yn,zn)
K2 = h*f1(xn+h/2, yn+K1/2, zn+L1/2)
L2 = h*f2(xn+h/2, yn+K1/2, zn+L1/2)
K3 = h*f1(xn+h/2, yn+K2/2, zn+L2/2)
L3 = h*f2(xn+h/2, yn+K2/2, zn+L2/2)
K4 = h*f1(xn+h/2, yn+K3/2, zn+L3/2)
L4 = h*f2(xn+h/2, yn+K3/2, zn+L3/2)

%Update both Y and Z:
Yn+1 = Yn + (1/6)*(K1+ 2*K2 + 2*K3 + K4)
Zn+1 = Zn + (1/6)*(K1+ 2*K2 + 2*K3 + K4)

Hope this helps !