solving 2nd order differential equation system in matlab

Oct 2015
4
0
iran
Hi
I have to solve a 2nd order differential equation system(9*9) in matlab and I have some problems with it:
eq1=f(x1,x2,.....,x9,x1'',x2'',....,x9'')=0
.
.
eq9

they are linear equations.
The initial conditions: x1(0)=x2(0)=....x9(0)=0, x1'(0)=.....x9'(0)=0
One of the problem is when I want to solve them by ode45 ,it needs to introduce x1'', x2'' ,.... in each equation.
what should I do? please help me.
 
Apr 2015
62
6
Sweden
Well ode45 only solves 1st order differential equations, so first off you have to rewrite every n'th order differential equation to a connected system of n first order differential equations. And then solve all these as your final system.
 
Oct 2015
4
0
iran
Thank you for your answer, I did that.
this is the basic program:
[T OUTPUT]=ode45(@odefunc,t,IC)

And I have to write an ode function like this:

function output=odefun (s,t)

x1=s(1)
y1=s(2) %y1=x1'
x2=s(3)
y2=s(4) %y2=x2'
.
.
.
x9=s(17)
y9=s(18) %y9=x9'


output=[y1;
f1(x1,x2,.....,x9,x1'',x2'',....,x9'');
y2;
f2;
.
.
.];
Is there any way to replace (x'')s by a variable that could be familiar to matlab? I can't replace them by y1' or diff(y1)!
 
Oct 2015
4
0
iran
How could I show to matlab that for example x1'' is the first order differential of y1?
thank you.
 
Apr 2015
62
6
Sweden
Well as an example given the following to differential equations:

\(\displaystyle \frac{d^{2} x}{d{t}^{2}} = -\frac{C_{d} \cdot A \cdot \rho}{2 \cdot m} \cdot \left( \left( \frac{dx}{dt} \right)^{2}+\left( \frac{dy}{dt}^{2} \right) \right) ^{\frac{1}{2}} \cdot \frac{dx}{dt} \)
\(\displaystyle \frac{d^{2} y}{d{t}^{2}} = -\frac{g}{m}-\frac{C_{d} \cdot A \cdot \rho}{2 \cdot m} \cdot \left( \left( \frac{dx}{dt} \right)^{2}+\left( \frac{dy}{dt}^{2} \right) \right) ^{\frac{1}{2}} \cdot \frac{dy}{dt}\)

Can be rewritten where you introduce following:

\(\displaystyle z_{1} = x\)
\(\displaystyle z_{2} =y\)
\(\displaystyle z_{3} = \frac{dx}{dt}\)
\(\displaystyle z_{4} = \frac{dy}{dt}\)

Then the system of first order differential equation from the above will be:

\(\displaystyle \frac{dz_{1}}{dt} = z_{3} \)
\(\displaystyle \frac{dz_{1}}{dt} = z_{4} \)
\(\displaystyle \frac{dz_{1}}{dt} = -\frac{C_{d} \cdot A \cdot \rho}{2 \cdot m} \cdot \left( \left( x_{3} \right)^{2}+\left( z_{4}^{2} \right) \right) ^{\frac{1}{2}} \cdot z_{3} \)
\(\displaystyle \frac{dz_{1}}{dt} = -\frac{g}{m}-\frac{C_{d} \cdot A \cdot \rho}{2 \cdot m} \cdot \left( \left( z_{3} \right)^{2}+\left( z_{4}^{2} \right) \right) ^{\frac{1}{2}} \cdot z_{4} \)


Then the correspondent matlab code for the above will look like this:

function [dzdt] = odefunB(t, z,Cd,A,p,m,g)

dzdt = [
z(3);
z(4);
- (Cd*A*p/(2*m)) * sqrt((z(3)^2 + z(4)^2))*z(3);
-g/m - (Cd*A*p/(2*m)) * sqrt((z(3)^2 + z(4)^2))*z(4)];
end

Where each line is one of the right hand sides of the first order differential equations above, and then all the conditions comes like this:

fineness = 1000;
startangel = 30;
initialspeed = 500;
m = 0.55;
Cd = 0.5;
d = 0.07;
A = (pi*d^2)/4;
p = 1.2041;
g = 9.82;
tspan = linspace(0,30,fineness);
z0 = [0 2 cos(startangel/180 * pi)*initialspeed sin(startangel /180 * pi)*initialspeed];

[t,z] = ode45(@(t,z) odefunB(t,z,Cd,A,p,m,g),tspan,z0);

Hope the example can help you :)
 
  • Like
Reactions: 1 person
Apr 2015
62
6
Sweden
ohh i should mention that m,Cd,d,A,p and g are all constants
 
Oct 2015
4
0
iran
I couldn't explain what exactly is the problem.

the functions that I have are like these:
f1=A11*x1''+A12*x2''+A13*x3''+......+A19*x9''+B11*x1+B12*x2+.....+B19*x9
f2=A21*x1''+A22*x2''+A23*x3''+......+A29*x9''+B21*x1+B22*x2+.....+B29*x9
.
.
.f9=A91*x1''+A92*x2''+A93*x3''+......+A99*x9''+B91*x1+B92*x2+.....+B99*x9

%A & B are constant

so in function as I said:
output=[s(2);
-(A12*x2''+A13*x3''+......+A19*x9''+B11*s(1)+B12*s(3)+.....+B19*s(17))/A11; % the problem is exactly these (x2'' , x3'' ..... x9'') in every EQ
s(4);
-(A21*x1''+A23*x3''+......+A29*x9''+B21*s(1)+B22*s(3)+.....+B29*s(17))/A22;
.
.
.];

according to ur exanple , in those equations there is only ONE 2nd order differential in each equation. so the solution and program is simple. But in this problem it's kind of complicated .
thank you again for ur time. If you know how should I solve it, please let me know.