Solving Matrix Differential Equations using 4th Order Runge-Kutta Method
15 views (last 30 days)
Show older comments
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
0 Comments
Answers (1)
Joel
on 29 Mar 2023
Edited: Joel
on 29 Mar 2023
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 !
0 Comments
See Also
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!