Why does Lorenz system trajectory look janky using Julia?

Question

I'm messing around with solving the Lorenz system over a long timespan (0 to 100,000) in both MATLAB and Julia in an attempt to make a crude runtime comparison (I'm using ode45 in MATLAB, and DP5 in Julia, with both abstol and reltol set to 1e-6, hoping that what they compute is equivalent).

I got that MATLAB took about 18 seconds, and produced this trajectory plot:

Julia timed in at 14.517 ms but I got the following warning:

Warning: Interrupted. Larger maxiters is needed. If you are using an integrator for non-stiff ODEs or an automatic switching algorithm (the default), you may want to consider using a method for stiff equations.

And as well as that, the trajectory plot looks... well... janky!

Why is this happening? Anyone got any ideas? Are ode45 and DP5 even comparable? Should I be coding up my own differential equation solver to make sure I can compare the two fairly?

I tried making sure that I used equivalent solvers, used the same relative and absolute tolerances so that the step sizes the ODE solvers took would be similar, but got very different looking results.

Edit: Thanks so much for the tips so far! As per suggestion, here is my code for the Lorenz trajectory in Julia (wouldn't be surprised if it had a bug that I haven't noticed yet...)

# Lorenz system ODE function
function lorenz!(du, u, p, t)
    sigma, rho, beta = p
    du[1] = sigma * (u[2] - u[1])
    du[2] = u[1] * (rho - u[3]) - u[2]
    du[3] = u[1] * u[2] - beta * u[3]
end

function test()
    # Lorenz system parameters as StaticArrays
    params = @SVector [10.0, 28.0, 8/3.0]

    # Initial conditions as StaticArrays
    u0 = [1.0, 0.0, 0.0]

    # Time span
    tspan = (0.0, 100000.0)

    # Solve the ODE using DP5
    quick_prob = ODEProblem(lorenz!, u0, 10, params)
    # function is compiled when we solve quick_prob
    quick_sol = solve(quick_prob, DP5(), abstol=1e-6, reltol=1e-6; dt = 0.1)

    # solve actual long problem and record execution time, excluding compilation time
    prob = ODEProblem(lorenz!, u0, tspan, params)
    @time sol = solve(prob, DP5(), abstol=1e-6, reltol=1e-6)

    # Plot the trajectory
    plot(sol, vars=(1, 2, 3), label="Lorenz System", xlabel="X", ylabel="Y", zlabel="Z", legend=:bottomright)
end

test()

Answer 1

This happens if the step size that is internally computed gets too large and the integrator only returns the internal nodes. If you sample a curve too sparsely and then plot the resulting polyline, you will get pointy angles at the nodes like you can see at the borders in the Julia plot.

This means on the one hand that the Julia integrator is so good that it can use rather large step sizes to meet the error target.

On the other hand, the Matlab solvers have a Refine parameter that defaults to 4, which means that for every internal step 3 points in the step interval are interpolated and added to the output, with the explicit aim to produce smooth plots fast, that is, using not very small error tolerances. One could say that this simulates the node density for a nodes-only output at error tolerances that are a factor 4^5=1024 smaller than the given ones.

_{From the Doc: The default value of Refine for most solvers is 1, but ode45 uses a default value of 4, while ode78 and ode89 use a default value of 8. These solvers use a larger default value to compensate for their tendency to take large steps.}

This is a useful middle position between just having the "raw" internal points and interpolating at externally given equidistant times regardless of how fast or slow the solution moves.

The Julia solvers have a similar parameter, but only in the event handling via continuous callbacks. There the parameter interp_points defaults to 10 points or segments per internal step to detect multiple sign changes inside the step interval.

Of course one could externally construct such an adaptive time sequence using interpolation. If there are N+1 times from N internal steps, make an equidistant list K4 with 4*N+1 points at distance 0.25 and construct T4 = interpolate(K4, K1, T), where (K1,T) is the table of integer indices and their times in the node list. Then use the dense output to get the states at the K4 times.

---

You could do your own solver experiments to see why you can get pointy plots. Take the circle DE x'=-y, y'=x, select your method and search for the step size where the solution plot still looks like a closed circle. With the low-order methods this requires still enough points that the result looks round. With RK4 it starts to get pointy with only about 15 points.

Answer 2

ode45 and DP5 are comparable. You don't need to write your own ode solver. most likely your julia code for the integrator has a bug. The easiest way to check this is to call your integrand on a random input in both Matlab and Julia to make sure you actually have the same equation.