Complex solutions in ODE solver
Mostrar comentarios más antiguos
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?
Respuestas (1)
x^12.94 can give complex values for the derivatives if x < 0.
11 comentarios
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
el 5 de Sept. de 2023
Torsten
el 5 de Sept. de 2023
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
el 6 de Sept. de 2023
Torsten
el 6 de Sept. de 2023
Then you should check your equations for errors. ode45 only digests what you feed to it.
Dimitris
el 7 de Sept. de 2023
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
el 10 de Sept. de 2023
Torsten
el 10 de Sept. de 2023
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}')
Categorías
Más información sobre Ordinary Differential Equations en Centro de ayuda y File Exchange.
el 4 de Sept. de 2023
el 10 de Sept. de 2023
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
