2nd order numerical differential equation system solving
2 次查看(过去 30 天)
显示 更早的评论
Hi!
Could you guys please help me with the following 2nd order equation system?
- G=6.673*10^-11;
- m1=1; m2=2; m3=3;
- syms x1(t) x2(t) x3(t);
- syms y1(t) y2(t) y3(t);
- syms u1(t) u2(t) u3(t);
- syms v1(t) v2(t) v3(t);
- %Körper 1/Mass1
- ode1 = u1==diff(x1,t);
- ode2 = v1==diff(y1,t);
- 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);
- 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);
- %Körper 2/Mass2
- ode5 = u2==diff(x2,t);
- ode6 = v2==diff(y2,t);
- 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);
- 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);
- %Körper 3/Mass3
- ode9 = u3==diff(x3,t);
- ode10 = v3==diff(y3,t);
- 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);
- 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);
- cond1 = x1(0) == 0;
- cond2 = x2(0) == 1;
- cond3 = x3(0) == 2;
- cond4 = y1(0) == 5;
- cond5 = y2(0) == 4;
- cond6 = y3(0) == 3;
- cond7 = u1(0) == 1;
- cond8 = u2(0) == 1;
- cond9 = u3(0) == 1;
- cond10 = v1(0) == 1;
- cond11 = v2(0) == 1;
- cond12 = v3(0) == 1;
- conds = [cond1; cond2; cond3; cond4; cond5; cond6; cond7; cond8; cond9; cond10; cond11; cond12];
- 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!
0 个评论
采纳的回答
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 个评论
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
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)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Equation Solving 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!