2nd order numerical differential equation system solving

1 次查看（过去 30 天）
ahkrit
ahkrit 2017-10-4
评论： Torsten 2017-10-11
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!

请先登录，再回答此问题。

采纳的回答

Torsten
Torsten 2017-10-6
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 个评论
ahkrit
ahkrit 2017-10-10
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
Torsten 2017-10-11
https://de.mathworks.com/help/symbolic/solve-differential-algebraic-equations.html#bvh12tx-2
might help.
Best wishes
Torsten.

请先登录，再进行评论。

更多回答（1 个）

Josh Meyer
Josh Meyer 2017-10-5
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 个评论
ahkrit
ahkrit 2017-10-5
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!

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by