Solve a system of 3 second order ODE

5 views (last 30 days)
Nicolas Caro
Nicolas Caro on 17 May 2021
Commented: James Tursa on 18 May 2021
Hello !
I'd like to implement and solve this second order ODE (Ordinary Differencial Equation) in a matlab program.
What I have for now is a program which can solve a first order ODE. I know that I have to split my second order ODE into 2 first order but I don't know how to implement it in Matlab. Can you help me set up the code for this ?
This is my code for now
Main :
h = 0.01; %Time step
y0 = 1; %This was for an example. For my problem, correct Initial Conditions are at next line
%y0 = [0,0,0,1,0,0]; %Initial condition dx1/dt = 1 m.s all others at 0
tfinal = 20;
tarray = 0:h:tfinal;
ytRK4 = RK4(y0,h,tfinal);
RK4.m :
function yt = RK4(y0,h,tfinal)
yt = y0;
%here I'll need to do : yt = zeros(numel(y0),length(h:tfinal)); for Memory Allocation
i = 2;
for t = h : h : tfinal %RK4 loop (going to change)
k1 = f(t-h , yt(i-1));
k2 = f(t - (0.5*h) , yt(i-1)+0.5*k1*h);
k3 = f(t - (0.5*h) , yt(i-1)+0.5*k2*h);
k4 = f(t , yt(i-1)+(k3 * h));
yt(i) = yt(i-1) + (1/6 * (k1 + 2*k2 + 2*k3 + k4)* h);
i = i + 1;
f.m :
function grad = f(t,y)
grad = (0.1).^t - 3*y;
I tried to comment as much as possible so you know my thoughts.
Tell me if some things aren't clear.
Thanks in advance for the time you'll take to respond :)

Answers (1)

James Tursa
James Tursa on 17 May 2021
Edited: James Tursa on 17 May 2021
You will need to change your scalar equations to vector equations, and you will need to change your derivative function to a vector function. For the derivative function, use the states x1, x2, x3, x1dot, x2dot, and x3dot in a 6-element vector called y. Then the code would look something like this:
function dy = deriv(t,y,k1,k2,k3,m1,m2,m3)
x1 = y(1);
x2 = y(2);
x3 = y(3);
x1dot = y(4);
x2dot = y(5);
x3dot = y(6);
x1dotdot = ____; % you fill this in
x2dotdot = ____; % you fill this in
x3dotdot = ____; % you fill this in
dy = [x1dot;x2dot;x3dot;x1dotdot;x2dotdot;x3dotdot];
You could define your f function in the RK4 routine as (where k1, k2, k3, m1 ,m2, m3 are previously defined):
f = @(t,y) deriv(t,y,k1,k2,k3,m1,m2,m3)
Your RK4 code then needs to be vectorized for the 6-element states. E.g., replace yt(i) with yt(:,i), and replace yt(i-1) with yt(:,i-1). Also the initialization of yt could be this:
yt = zeros(6,ceil(tfinal/h));
yt(:,1) = y0;
James Tursa
James Tursa on 18 May 2021
No need to post it if you have it working and have no further questions.

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by