Solve a system of 3 second order ODE

조회 수: 5(최근 30일)
Nicolas Caro 2021년 5월 17일
댓글: James Tursa 2021년 5월 18일
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 :
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 :)

댓글을 달려면 로그인하십시오.

답변(1개)

James Tursa 2021년 5월 17일
편집: James Tursa 2021년 5월 17일
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표시숨기기 이전 댓글 수: 1
James Tursa 2021년 5월 18일
No need to post it if you have it working and have no further questions.

댓글을 달려면 로그인하십시오.

Community Treasure Hunt

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

Start Hunting!

Translated by