Complex solutions in ODE solver

5 次查看（过去 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