Solve a system of 3 second order ODE
10 views (last 30 days)
Show older comments
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 :
clc;
clear;
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;
end
end
f.m :
function grad = f(t,y)
grad = (0.1).^t - 3*y;
end
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 :)
0 Comments
Answers (1)
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];
end
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;
2 Comments
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!