# 2nd order numerical differential equation system solving

5 views (last 30 days)
ahkrit on 4 Oct 2017
Commented: Torsten on 11 Oct 2017
Hi!
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);
7. %Körper 1/Mass1
8. ode1 = u1==diff(x1,t);
9. ode2 = v1==diff(y1,t);
10. 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);
11. 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);
12. %Körper 2/Mass2
13. ode5 = u2==diff(x2,t);
14. ode6 = v2==diff(y2,t);
15. 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);
16. 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);
17. %Körper 3/Mass3
18. ode9 = u3==diff(x3,t);
19. ode10 = v3==diff(y3,t);
20. 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);
21. 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);
22. cond1 = x1(0) == 0;
23. cond2 = x2(0) == 1;
24. cond3 = x3(0) == 2;
25. cond4 = y1(0) == 5;
26. cond5 = y2(0) == 4;
27. cond6 = y3(0) == 3;
28. cond7 = u1(0) == 1;
29. cond8 = u2(0) == 1;
30. cond9 = u3(0) == 1;
31. cond10 = v1(0) == 1;
32. cond11 = v2(0) == 1;
33. cond12 = v3(0) == 1;
34. conds = [cond1; cond2; cond3; cond4; cond5; cond6; cond7; cond8; cond9; cond10; cond11; cond12];
35. 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!

Torsten on 6 Oct 2017
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.
##### 2 CommentsShowHide 1 older comment
Torsten on 11 Oct 2017
https://de.mathworks.com/help/symbolic/solve-differential-algebraic-equations.html#bvh12tx-2
might help.
Best wishes
Torsten.

Josh Meyer on 5 Oct 2017
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 on 5 Oct 2017
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.