solve ode system with ode45

Hello everyone, I would like to solve a system of differential equations using ode45, but I don't know how to proceed :

• d^2 (x)/dt^2 = a * (d(x)/dt - d(y)/dt) + b * x^3
• d^2 (y)/dt^2 = c * (d(y)/dt - d(x)/dt) + b * y^3

of course my system is much more complicated. Thanks a lot of your help.

Answer 1

[Massimo Zanetti]

So, basically the solution can be implemented as follows

%parameters
A=1;
B=1;
C=1;
%time interval and initial conditions
t_interval = [0,10];
init_cond = [0,0,0,0]';
%solution
[t,y] = ode45(@(t,Y) odefcn(t,Y,A,B,C) , t_interval , init_cond);
%plot
plot(t,y(:,1),'b',t,y(:,2),'r');

where your system of equations transformed into ode by sobstitution x'=z, y'=w, and the order of variables is Y=[x,y,z,w]:

function dYdt = odefcn(t,Y,A,B,C)

dYdt = [ Y(3);
         Y(4);
         A*(Y(3)-Y(4)) + B*Y(1)^3;
         C*(Y(4)-Y(3)) + B*Y(2)^3];

end

Of course, you will notice that starting by null conditions (init_cond=[0,0,0,0]) the solution you get is x(t)=y(t)=0 for each t. Because they solve the equation. Change initial conditions according to your original problem.

[Bechara Saade]

Thank you! It's working, but if we have null initial conditions how we can solve it?

[Sofiya Vyshnya]

This may not be the perfect solution, but could you set the initial conditions to some really small value, such as 1E-10? This value is approximately equal to 0, but it won't make your x(t) and y(t) get zeroed out.

[Abbas Abouei]

Sir thanks for the comment, I am trying to solve a system of coupled equation only. i used your way. i can get the output but it seems that it is not right, the matlab is busy for long time and no output.it seems cpu also dose not occupied by matlab. coul you please help me through it?

the code is as follow:

clear all; close all; clc;

g=9.8;

gamma=0.0613;

M=8.7*10^-8;

v=5;

landa=2.5;

eta=2.5;

G=0.01;

b=1000;

%time interval and initial conditions

t_interval = [0,1];

init_cond = [0.2,1,1,0]';

%solution

[t,y] = ode45(@(t,Y) odefcn(t,Y,g,gamma,M,v,eta,G,b,landa) , t_interval , init_cond);

%plot

plot(t,y(:,1),'b');

function dYdt = odefcn(t,Y,g,gamma,M,v,eta,G,b,landa)

dYdt = [ Y(4);

2*Y(4)*Y(2)/Y(1)-(Y(2)-1)/(landa*(1-(Y(2)+2*Y(3))/b));

-1*Y(4)*Y(3)/Y(1)-(Y(3)-1)/(landa*(1-(Y(2)+2*Y(3))/b));

g-(gamma/M)*(v*3.14/Y(1))^0.5-3*eta*v*Y(4)/(M*Y(1)^2)-(v*G/M*Y(1))*(Y(2)-Y(1))/(1-(Y(2)+2*Y(3))/b);];

end

the original equuation is like

the initial condition for example is L(0)=0.1 L'(0)=0, azz(0)=1 arr(0)=1

I simplify the equation like L'=W and Y=[L,azz,arr,W] so (L(0)=0.1,azz(0)=1 arr(0)=1,W(0)=0)

anyhelp is appreciated

Answer 2

[Steven Lord]

See this video by Cleve or the "Solve Nonstiff Equation" example on the documentation page for the ode45 function for techniques you can use.

Answer 3

[piranha007]

Thanks guys for your help.