4 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 :)

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;

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

Start Hunting!