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?

Answer 1

[Torsten]

x^12.94 can give complex values for the derivatives if x < 0.

[Dimitris]

It could be, though all numbers are positive for t=0 and they remain positive for a lot of seconds. I get a complex solution from the beginning. Other than that, the 12.94 exponent isn't just a random number that I put there. It has an important physical meaning and I cannot really change it.

[Torsten]

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]);

ans = 1

ans = 1.0000

ans = 1.0000

ans = 1.0000

ans = 1.0000

ans = 1.0000

ans = 1.0000

ans = 1.0000

ans = 1.0000

ans = 1.0000

ans = 1.0000

ans = 1.0000

ans = 1.0000

ans = 1.0000

ans = 1.0000

ans = 1.0000

ans = 1.0000

ans = 1.0000

ans = 1.0000

ans = 1.0000

ans = 1.0000

ans = 1.0000

ans = 1.0000

ans = 1.0000

ans = 1.0000

ans = 1.0000

ans = 1.0000

ans = 1.0000

ans = 1.0000

ans = 1.0000

ans = 1.0000

ans = 1.0000

ans = 1.0000

ans = 1.0000

ans = 1.0000

ans = 1.0000

ans = 1.0000

ans = 1.0000

ans = 1.0000

ans = 1.0000

ans = 1.0000

ans = 1.0000

ans = 1.0000

ans = 0.9999

ans = 0.9999

ans = 0.9998

ans = 0.9998

ans = 0.9998

ans = 0.9998

ans = 0.9997

ans = 0.9996

ans = 0.9992

ans = 0.9991

ans = 0.9990

ans = 0.9990

ans = 0.9983

ans = 0.9979

ans = 0.9960

ans = 0.9957

ans = 0.9952

ans = 0.9952

ans = 0.9914

ans = 0.9895

ans = 0.9800

ans = 0.9784

ans = 0.9762

ans = 0.9763

ans = 0.9573

ans = 0.9479

ans = 0.9005

ans = 0.8919

ans = 0.8815

ans = 0.8818

ans = 0.8228

ans = 0.7934

ans = 0.6453

ans = 0.6172

ans = 0.5844

ans = 0.5873

ans = 0.5179

ans = 0.4832

ans = 0.3097

ans = 0.2789

ans = 0.2403

ans = 0.2403

ans = 0.1710

ans = 0.1363

ans = -0.0372

ans = -0.0680

ans = -0.1066 + 0.0000i

ans = -0.1066 - 0.0000i

ans = -0.1760 + 0.0000i

ans = -0.2107 - 0.0000i

ans = -0.3842 + 0.0000i

ans = -0.4150 + 0.0000i

ans = -0.4535 + 0.0000i

ans = -0.4535 + 0.0000i

ans = -0.5229 + 0.0000i

ans = -0.5576 + 0.0000i

ans = -0.7311 + 0.0000i

ans = -0.7619 + 0.0000i

ans = -0.8005 + 0.0000i

ans = -0.8005 + 0.0000i

ans = -0.8699 + 0.0000i

ans = -0.9047 + 0.0000i

ans = -1.0804 + 0.0005i

ans = -1.1123 + 0.0007i

ans = -1.1546 + 0.0014i

ans = -1.1544 + 0.0014i

ans = -1.2296 + 0.0029i

ans = -1.2743 + 0.0053i

ans = -1.7082 + 0.0871i

ans = -1.6454 - 0.0398i

ans = 4.9333 +76.1697i

ans = -2.9668e+05 - 9.1796e+04i

ans = -1.1619 + 0.0016i

ans = -1.1657 + 0.0017i

ans = -1.1851 + 0.0022i

ans = -1.1887 + 0.0023i

ans = -1.1931 + 0.0024i

ans = -1.1929 + 0.0024i

ans = -1.2008 + 0.0026i

ans = -1.2049 + 0.0028i

ans = -1.2257 + 0.0037i

ans = -1.2297 + 0.0040i

ans = -1.2345 + 0.0042i

ans = -1.2342 + 0.0041i

ans = -1.2706 + 0.0062i

ans = -1.2926 + 0.0085i

ans = -1.4538 + 0.0481i

ans = -1.4924 + 0.0410i

ans = -1.5142 + 0.5189i

ans = -1.6149 + 3.4362i

ans = -1.2378 + 0.0043i

ans = -1.2397 + 0.0045i

ans = -1.2490 + 0.0051i

ans = -1.2508 + 0.0052i

ans = -1.2529 + 0.0053i

ans = -1.2528 + 0.0053i

ans = -1.2567 + 0.0056i

ans = -1.2586 + 0.0058i

ans = -1.2686 + 0.0066i

ans = -1.2704 + 0.0068i

ans = -1.2727 + 0.0070i

ans = -1.2726 + 0.0070i

ans = -1.2931 + 0.0090i

ans = -1.3051 + 0.0108i

ans = -1.3782 + 0.0302i

ans = -1.3990 + 0.0354i

ans = -1.4155 + 0.0598i

ans = -1.4010 + 0.0652i

ans = -1.2793 + 0.0076i

ans = -1.2829 + 0.0080i

ans = -1.3012 + 0.0105i

ans = -1.3049 + 0.0112i

ans = -1.3091 + 0.0119i

ans = -1.3087 + 0.0116i

ans = -1.3166 + 0.0129i

ans = -1.3208 + 0.0138i

ans = -1.3431 + 0.0197i

ans = -1.3479 + 0.0216i

ans = -1.3532 + 0.0233i

ans = -1.3523 + 0.0225i

ans = -1.3653 + 0.0273i

ans = -1.3726 + 0.0313i

ans = -1.4074 + 0.0673i

ans = -1.4199 + 0.0837i

ans = -1.4200 + 0.0818i

ans = -1.4142 + 0.0798i

ans = -1.3561 + 0.0239i

ans = -1.3580 + 0.0247i

ans = -1.3679 + 0.0294i

ans = -1.3696 + 0.0305i

ans = -1.3719 + 0.0317i

ans = -1.3719 + 0.0314i

ans = -1.3759 + 0.0336i

ans = -1.3779 + 0.0350i

ans = -1.3876 + 0.0426i

ans = -1.3890 + 0.0446i

ans = -1.3912 + 0.0464i

ans = -1.3915 + 0.0458i

ans = -1.3973 + 0.0515i

ans = -1.3997 + 0.0551i

ans = -1.4064 + 0.0727i

ans = -1.4021 + 0.0769i

ans = -1.4042 + 0.0821i

ans = -1.4113 + 0.0797i

ans = -1.3960 + 0.0501i

ans = -1.3979 + 0.0528i

ans = -1.4050 + 0.0658i

ans = -1.4040 + 0.0685i

ans = -1.4056 + 0.0719i

ans = -1.4085 + 0.0713i

ans = -1.4107 + 0.0768i

ans = -1.4112 + 0.0794i

ans = -1.4139 + 0.0907i

ans = -1.4146 + 0.0900i

ans = -1.4144 + 0.0925i

ans = -1.4141 + 0.0965i

ans = -1.4145 + 0.1020i

ans = -1.4144 + 0.1044i

ans = -1.4153 + 0.1159i

ans = -1.4179 + 0.1177i

ans = -1.4178 + 0.1197i

ans = -1.4140 + 0.1199i

ans = -1.4137 + 0.1245i

ans = -1.4134 + 0.1264i

ans = -1.4129 + 0.1364i

ans = -1.4136 + 0.1391i

ans = -1.4135 + 0.1410i

ans = -1.4122 + 0.1394i

ans = -1.4117 + 0.1442i

ans = -1.4115 + 0.1462i

ans = -1.4107 + 0.1568i

ans = -1.4109 + 0.1600i

ans = -1.4107 + 0.1620i

ans = -1.4103 + 0.1597i

ans = -1.4099 + 0.1646i

ans = -1.4098 + 0.1666i

ans = -1.4094 + 0.1772i

ans = -1.4092 + 0.1804i

ans = -1.4092 + 0.1824i

ans = -1.4094 + 0.1802i

ans = -1.4094 + 0.1851i

ans = -1.4095 + 0.1871i

ans = -1.4099 + 0.1978i

ans = -1.4096 + 0.2008i

ans = -1.4097 + 0.2028i

ans = -1.4104 + 0.2007i

ans = -1.4109 + 0.2056i

ans = -1.4112 + 0.2075i

ans = -1.4127 + 0.2180i

ans = -1.4125 + 0.2209i

ans = -1.4129 + 0.2229i

ans = -1.4138 + 0.2209i

ans = -1.4148 + 0.2255i

ans = -1.4155 + 0.2273i

ans = -1.4182 + 0.2371i

ans = -1.4182 + 0.2399i

ans = -1.4189 + 0.2417i

ans = -1.4197 + 0.2398i

ans = -1.4213 + 0.2439i

ans = -1.4222 + 0.2456i

ans = -1.4261 + 0.2543i

ans = -1.4263 + 0.2569i

ans = -1.4272 + 0.2585i

ans = -1.4279 + 0.2565i

ans = -1.4299 + 0.2601i

ans = -1.4309 + 0.2614i

ans = -1.4356 + 0.2688i

ans = -1.4361 + 0.2712i

ans = -1.4371 + 0.2725i

ans = -1.4375 + 0.2706i

ans = -1.4398 + 0.2735i

ans = -1.4408 + 0.2746i

ans = -1.4459 + 0.2806i

ans = -1.4467 + 0.2826i

ans = -1.4476 + 0.2837i

ans = -1.4477 + 0.2819i

ans = -1.4500 + 0.2842i

ans = -1.4510 + 0.2850i

ans = -1.4561 + 0.2898i

ans = -1.4571 + 0.2914i

ans = -1.4581 + 0.2923i

ans = -1.4578 + 0.2908i

ans = -1.4586 + 0.2914i

ans = -1.4590 + 0.2917i

ans = -1.4608 + 0.2931i

ans = -1.4612 + 0.2934i

ans = -1.4616 + 0.2937i

ans = -1.4616 + 0.2936i

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]

Yes, I see. Is there though a way to fully get rid of complex solutions? Cause if I understood your code correctly, I wrote

[t,ya] = ode45(@f,[0 30],[12075 1740 0 80000]);

function dydt = f(t,y)

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

Which still ends up giving me a lot of complex numbers to be honest.

[Torsten]

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]

Yes, I understand. Though, for it to become negative, y(3) must be greater than about 103000. My initial conditions define y(3) as 80000, and then it should be decreasing. So there shouldn't arise any case that 1-0.0000096*y(3) ever becomes negative.

[Torsten]

Then you should check your equations for errors. ode45 only digests what you feed to it.

[Dimitris]

Maybe there is an error. 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.

[Torsten]

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]

Oh, I didn't know it was based in the R-K method. It could then be the way the solver in interpreting the equations. I'll see if I can somehow change their form.

[Torsten]

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]

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}')