# 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!