Complex solutions in ODE solver

7 ビュー (過去 30 日間)
Dimitris
Dimitris 2023 年 9 月 4 日
コメント済み: Sam Chak 2023 年 9 月 10 日
Hi guys. I am trying to solve a system of 4 ODEs to run a physics simulation. The code is:
f = @(t,y) [-0.0000017*(1-0.0000097*y(3))^12.94*(y(1)^2+y(2)^2)/sqrt((y(2)/y(1))^2+1);+9.59-0.0000017*(1-0.0000097*y(3))^12.94*y(2)*(y(1)^2+y(2)^2)/(y(1)*sqrt((y(2)/y(1))^2+1));y(1);-y(2)+y(4)/(sqrt(6371000^2+y(4)^2))*y(1);]; % y(1)=vx, y(2)=vy, y(3)=y
[t,ya] = ode45(f,[0 30],[12075 1740 0 80000]); % vx(0), vy(0), x(0), y(0)
Now, the solver works pretty well if the first two initial conditions are small. However, when I get into big numbers, like here, the solver shows complex solutions for some reason. This shouldn't be happening, as every root is in the form of sqrt(a^2+b^2), which is always positive, and every denominator is also non-zero. Is there a way to fix this?

サインインしてこの質問に回答する。

回答 (1 件)

Torsten
Torsten 2023 年 9 月 4 日
移動済み: Torsten 2023 年 9 月 4 日
x^12.94 can give complex values for the derivatives if x < 0.
  11 件のコメント
9 件の古いコメントを表示
Torsten
Torsten 2023 年 9 月 10 日
It doesn't matter what solution method for the system of differential equations you use. All of them should yield the same solution. There is no "interpreting of the equations" - you supply the time derivatives, the solver solves for the functions for which you supplied the time derivatives.
Sam Chak
Sam Chak 2023 年 9 月 10 日
Ran in:
Hi @Dimitris
Sometimes, studying the system may help to understand the behavior of the state trajectories.
From the first differential equation
if we assume that , and modify some constants, then two state equations can be analyzed as a nonlinear second-order differential equation.
From the stream plot, when is positive, then rapidly increases. In the actual scenario, as increases, then this term '' will become negative at one point, as shown by @Torsten.
[x, y] = meshgrid (-3:0.15:3, -3:0.15:3);
u = y;
v = - 0.17*(1 - 0.97*x).^(1).*(y.^2);
l = streamslice(x, y, u, v);
axis tight
xlabel('y_{3}'), ylabel('y_{1}')
You really need to check the equations again for errors. One thing to observe is that if I change , then the first three states show stable convergent trajectories. No more complex numbers.
f = @(t,y) [ - 0.0000017*(1 - 0.0000097*y(3))^12.94*(y(1)^2 + y(2)^2)/sqrt((y(2)/y(1))^2 + 1);
9.59 - 0.0000017*(1 - 0.0000097*y(3))^12.94*y(2)*(y(1)^2 + y(2)^2)/(y(1)*sqrt((y(2)/y(1))^2 + 1));
- y(1); % <-- change the sign here
- y(2) + y(4)/(sqrt(6371000^2 + y(4)^2))*y(1)];
tspan = [0 300];
y0 = [12075 1740 0 80000];
[t, ya] = ode45(f, tspan, y0);
plot(t, ya), grid on
xlabel('t')
legend('y_{1}', 'y_{2}', 'y_{3}', 'y_{4}')

サインインしてコメントする。

カテゴリ

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