ODE solver giving wrong results

General Usage
1

Hey everyone!
I want to solve f(u,p,t) = An*u + b, where:

An = [0.0+0.0im 1.0+0.0im
5.0+0.0im -1.0+0.0im],
b = [1.0+0.0im, 0.0+0.0im],
x = [0.0+0.0im, 1.0+0.0im] (initial condition)

Which means that I have a system like this:
x’ = y+1, x(0) = 0
y’ = 5x-y, y(0) = 1

And I get the following ODE: x"(t)+x’(t)-5x(t)=1. The analytic solution gives me x(1) = 3.1259897764930082
However, when I solve it via Julia’s ODE solver, I get x(1) = 4.21229.

Am I missing something? Here’s my code:

An = [0.0+0.0im 1.0+0.0im
	5.0+0.0im -1.0+0.0im]
b = [1.0+0.0im, 0.0+0.0im]
x = [0.0+0.0im, 1.0+0.0im]
tspan = (0.0,1.0)
f(u,p,t) = An*u + b;
prob = ODEProblem(f, x, tspan)
sol = solve(prob, Vern7(), dt = 0.1, adaptive = false, reltol=1e-8)
2

Seems strange, when I simply copy and paste your code, i get

julia> sol = solve(prob, Vern7(), dt = 0.1, adaptive = false, reltol=1e-8);

julia> sol(1)
2-element Vector{ComplexF64}:
 3.125989776501887 + 0.000000000000000im
 5.058513100378909 + 0.000000000000000im

julia>

which corresponds well enough to your analytical solution

2 Likes
3

Another thing even though unrelated to your original problem: Setting adaptive = false, should basically invalidate your reltol = 1e-8 option, since the ODE solver is only allowed to go in steps of 0.1 without adjusting the stepsize due to error. It should always be safer to use the default, which is adaptive = true.

2 Likes
4

The analytic solution to your equation is

x(t) = \frac{1-\sqrt{21}}{10}\exp{\left(-\frac{1+\sqrt{21}}{2}t\right)} + \frac{1+\sqrt{21}}{10}\exp{\left( \frac{\sqrt{21}-1}{2}t\right)} - \frac{1}{5}

which gives x(0) = 0 and x(1) \approx 3.1259....

Two comments:

  • Why do you use complex numbers in An, b, and x (initial value)?
  • The solution represents an unstable system, where numerical errors easily can “explode”.
2 Likes
5

Thank you for replying.
Yes, I tried again later, and I got the same result this time. I guess that was a bug or something.

6

Thanks!
I’m using complex numbers only because the other parts of my code require complex inputs.
Also, yes - I’m mainly focusing on stable systems to get lower errors (which is not this case though).