Complex solutions in ODE solver
古いコメントを表示
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 件)
x^12.94 can give complex values for the derivatives if x < 0.
11 件のコメント
If you define your time derivatives in a function instead of a function handle, you can control when and how your derivatives become complex:
[t,ya] = ode45(@f,[0 30],[12075 1740 0 80000]);
function dydt = f(t,y)
1-0.0000097*y(3)
dydt = [-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)];
end
Dimitris
2023 年 9 月 5 日
Torsten
2023 年 9 月 5 日
I just wanted to show you when the complex numbers arise, namely at the time when (1-0.0000097*y(3)) changes sign from positive to negative. Of course there are ways to avoid this - the question is how to define (1-0.0000097*y(3)) in this case in order to make sense for your simulated functions.
Dimitris
2023 年 9 月 6 日
Torsten
2023 年 9 月 6 日
Then you should check your equations for errors. ode45 only digests what you feed to it.
Dimitris
2023 年 9 月 7 日
Or maybe I could potentially use the Runge-Kutta method to solve them, if I can find a way to write a code for it.
So you think your solver will be better than the MATLAB solver (especially because ode45 is a Runge-Kutta solver) ? That's ... hubris.
Depending on the background of your equations, x^a for negative x is sometimes interpreted as -abs(x)^a. Summarized: x^a = sign(x)*abs(x)^a. Test it whether it makes sense for your case.
Dimitris
2023 年 9 月 10 日
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.
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.
, 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.
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.
, 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}')
カテゴリ
ヘルプ センター および 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 
