2nd order numerical differential equation system solving

Hi!

Could you guys please help me with the following 2nd order equation system?

1. G=6.673*10^-11;
2. m1=1; m2=2; m3=3;
3. syms x1(t) x2(t) x3(t);
4. syms y1(t) y2(t) y3(t);
5. syms u1(t) u2(t) u3(t);
6. syms v1(t) v2(t) v3(t);
8. %Körper 1/Mass1
9. ode1 = u1==diff(x1,t);
10. ode2 = v1==diff(y1,t);
11. ode3 = diff(u1,t)*m1==((G*m1*m2)/((x2-x1)^2+(y2-y1)^2)^(3/2))*(x2-x1)+((G*m1*m3)/((x3-x1)^2+(y3-y1)^2)^(3/2))*(x3-x1);
12. ode4 = diff(v1,t)*m1==((G*m1*m2)/((x2-x1)^2+(y2-y1)^2)^(3/2))*(y2-y1)+((G*m1*m3)/((x3-x1)^2+(y3-y1)^2)^(3/2))*(y3-y1);
14. %Körper 2/Mass2
15. ode5 = u2==diff(x2,t);
16. ode6 = v2==diff(y2,t);
17. ode7 = diff(u2,t)*m2==((G*m2*m3)/((x3-x2)^2+(y3-y2)^2)^(3/2))*(x3-x2)+((G*m1*m2)/((x1-x2)^2+(y1-y2)^2)^(3/2))*(x1-x2);
18. ode8 = diff(v2,t)*m2==((G*m2*m3)/((x3-x2)^2+(y3-y2)^2)^(3/2))*(y3-y2)+((G*m1*m2)/((x1-x2)^2+(y1-y2)^2)^(3/2))*(y1-y2);
20. %Körper 3/Mass3
21. ode9 = u3==diff(x3,t);
22. ode10 = v3==diff(y3,t);
23. ode11 = diff(u3,t)*m3==((G*m3*m1)/((x1-x3)^2+(y1-y3)^2)^(3/2))*(x1-x3)+((G*m3*m2)/((x2-x3)^2+(y2-y3)^2)^(3/2))*(x2-x3);
24. ode12 = diff(v3,t)*m3==((G*m3*m1)/((x1-x3)^2+(y1-y3)^2)^(3/2))*(y1-y3)+((G*m3*m2)/((x2-x3)^2+(y2-y3)^2)^(3/2))*(y2-y3);
26. cond1 = x1(0) == 0;
27. cond2 = x2(0) == 1;
28. cond3 = x3(0) == 2;
30. cond4 = y1(0) == 5;
31. cond5 = y2(0) == 4;
32. cond6 = y3(0) == 3;
34. cond7 = u1(0) == 1;
35. cond8 = u2(0) == 1;
36. cond9 = u3(0) == 1;
38. cond10 = v1(0) == 1;
39. cond11 = v2(0) == 1;
40. cond12 = v3(0) == 1;
42. conds = [cond1; cond2; cond3; cond4; cond5; cond6; cond7; cond8; cond9; cond10; cond11; cond12];
43. odes = [ode1; ode2; ode3; ode4; ode5; ode6; ode7; ode8; ode9; ode10; ode11; ode12];

I tried to solve it with dsolve. How could it be solved with ode45? Thanks in advance!

Answer 1

[Torsten]

function main
y0=[0; 5; 1; 1; 1; 4; 1; 1; 2; 3; 1; 1];
t0=0;
tfinal=10;
[T Y] = ode45(@odesNew,[t0 tfinal],y0)

function dy = odesNew(t,y)
G=6.673*10^-11;
m1=1; m2=2; m3=3;
dy=zeros(12,1);
x1=y(1);
x2=y(2);
x3=y(3);
y1=y(4);
y2=y(5);
y3=y(6);
u1=y(7);
u2=y(8);
u3=y(9);
v1=y(10);
v2=y(11);
v3=y(12);
%Körper 1/Mass1
dy(1)=u1;
dy(4)=v1;
dy(7)=(((G*m1*m2)/((x2-x1)^2+(y2-y1)^2)^(3/2))*(x2-x1)+((G*m1*m3)/((x3-x1)^2+(y3-y1)^2)^(3/2))*(x3-x1))/m1;
dy(10)=(((G*m1*m2)/((x2-x1)^2+(y2-y1)^2)^(3/2))*(y2-y1)+((G*m1*m3)/((x3-x1)^2+(y3-y1)^2)^(3/2))*(y3-y1))/m1;
%Körper 2/Mass2
dy(2)=u2;
dy(5)=v2;
dy(8)=(((G*m2*m3)/((x3-x2)^2+(y3-y2)^2)^(3/2))*(x3-x2)+((G*m1*m2)/((x1-x2)^2+(y1-y2)^2)^(3/2))*(x1-x2))/m2;
dy(11)=(((G*m2*m3)/((x3-x2)^2+(y3-y2)^2)^(3/2))*(y3-y2)+((G*m1*m2)/((x1-x2)^2+(y1-y2)^2)^(3/2))*(y1-y2))/m2;
%Körper 3/Mass3
dy(3)=u3;
dy(6)=v3;
dy(9)=(((G*m3*m1)/((x1-x3)^2+(y1-y3)^2)^(3/2))*(x1-x3)+((G*m3*m2)/((x2-x3)^2+(y2-y3)^2)^(3/2))*(x2-x3))/m3;
dy(12)=(((G*m3*m1)/((x1-x3)^2+(y1-y3)^2)^(3/2))*(y1-y3)+((G*m3*m2)/((x2-x3)^2+(y2-y3)^2)^(3/2))*(y2-y3))/m3;

Best wishes

Torsten.

[ahkrit]

Thanks a lot! I would like to just know, is there another way of solving this with the original form of the equations? I mean with more functions as i wrote it first.

Thanks in advance!

[Torsten]

https://de.mathworks.com/help/symbolic/solve-differential-algebraic-equations.html#bvh12tx-2

might help.

Best wishes

Torsten.

Answer 2

[Josh Meyer]

Use the Symbolic Math Toolbox function odeFunction to convert the odes variable into a function handle. Once you have that, you just need to construct a numeric vector of initial conditions y0 (similar to conds) and decide what time span to solve over. The syntax will be

[t,y] = ode45(@odesNew,[t0 tfinal],y0)

[ahkrit]

Hi! Could you write down exactly how it should look like for the first 4 equations, i just cannot figure it out how it should work.

Thanks in advance!