2nd order numerical differential equation system solving

5 views (last 30 days)
ahkrit
ahkrit on 4 Oct 2017
Commented: Torsten on 11 Oct 2017
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);
  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!

Accepted Answer

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

Sign in to comment.

More Answers (1)

Josh Meyer
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)
  1 Comment
ahkrit
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.
Thanks in advance!

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!