KenCarp47 is giving wrong result no matter how tight the tolerance is. Is there a way to resolve this?

New to Julia
21

Changed every single number in the code to BigFloat. The result is unchanged.

using DifferentialEquations
using OrdinaryDiffEq
using ODEInterfaceDiffEq
using Plots
using Sundials
using LinearSolve
using BenchmarkTools
using Base.Threads
const θ = BigFloat(1.0e-6)
const Energy=BigFloat(50.0e0)
const Δmsq = BigFloat(2.5e-3*(1.0e6/197.3269631))
const Hvac=(Δmsq/(4.0*Energy))*[[-cos(2.0*θ)+0im, sin(2.0*θ)+0im] [sin(2.0*θ)+0im, cos(2.0*θ)+0im]]
const c = BigFloat(3e5)
const mup = BigFloat(1.0e5)
@time function eoms!(du,u,p,t)
    println(t)
    AngleBins=p[1]
    for i = 1:AngleBins*8
	    du[i]=BigFloat(0.0e0)
    end
	Hself0 = zeros(Complex{BigFloat},2,2)
    Hself1=zeros(Complex{BigFloat},2,2)
    dcosth=BigFloat(2.0e0/AngleBins)
    costharr=(range(BigFloat(-1.0+dcosth/2.0e0),BigFloat(1.0-dcosth/2.0e0),AngleBins))
	#splock = SpinLock()
    #Threads.@threads
	for i = 0:AngleBins-1
		#lock(splock)
        rho=[[u[i*8+1]+0im, u[i*8+3]-1im*u[i*8+4]] [u[i*8+3]+1im*u[i*8+4],u[i*8+2]+0im]]
        rhobar=[[u[i*8+5]+0im, u[i*8+7]-1im*u[i*8+8]] [u[i*8+7]+1im*u[i*8+8],u[i*8+6]+0im]]
		Hself0=Hself0+(rho-rhobar)*dcosth
	    Hself1=Hself1+(rho-rhobar)*costharr[i+1]*dcosth
		#unlock(splock)
	end
	#Threads.@threads
	for i = 0:AngleBins-1
		rho=[[u[i*8+1]+0im, u[i*8+3]-1im*u[i*8+4]] [u[i*8+3]+1im*u[i*8+4],u[i*8+2]+0im]]
        rhobar=[[u[i*8+5]+0im, u[i*8+7]-1im*u[i*8+8]] [u[i*8+7]+1im*u[i*8+8],u[i*8+6]+0im]]
		Hself=mup*(Hself0-Hself1*costharr[i+1])
    	dudt=-c*1.0e0im*((Hvac+Hself)*rho-rho*(Hvac+Hself))
        dudtbar=-c*1.0e0im*((-BigFloat(1.0e0)*Hvac+Hself)*rhobar-rhobar*(BigFloat(-1.0e0)*Hvac+Hself))
        du[i*8+1]=real(dudt[1,1])
        du[i*8+2]=real(dudt[2,2])
        du[i*8+3]=real(dudt[1,2])
        du[i*8+4]=imag(dudt[1,2])
    	du[i*8+5]=real(dudtbar[1,1])
        du[i*8+6]=real(dudtbar[2,2])
        du[i*8+7]=real(dudtbar[1,2])
        du[i*8+8]=imag(dudtbar[1,2])
    end
end
AngleBins=100
size=8*AngleBins
dcosth=BigFloat(2.0e0/AngleBins)
costharr=range(BigFloat(-1.0+dcosth/2.0e0),BigFloat(1.0-dcosth/2.0e0),AngleBins)
u0=zeros(BigFloat,size)
for i=0:AngleBins-1
    u0[8*i+1]=BigFloat(0.5e0)
    u0[8*i+5]=BigFloat(0.47e0)+BigFloat(0.05)*exp(-(costharr[i+1]-BigFloat(1.0e0))*(costharr[i+1]-BigFloat(1.0e0)))
end
p=[AngleBins]

tspan = (BigFloat(0.0),BigFloat(1e-6))
prob = ODEProblem(eoms!,u0,tspan,p)
#reltolp=1e-10
#abstolp=1e-10
reltolp=1e-20
abstolp=1e-20
#@time sol1 = solve(prob,ARKODE(Sundials.Explicit(),etable = Sundials.DORMAND_PRINCE_7_4_5),reltol=reltolp,abstol=abstolp,saveat=1e-8)
#@time sol2 = solve(prob,CVODE_BDF(linear_solver=:GMRES,krylov_dim=30),reltol=reltolp,abstol=abstolp,saveat=1e-8)
#@time sol3 = solve(prob,CVODE_Adams(linear_solver=:GMRES,krylov_dim=30),reltol=reltolp,abstol=abstolp,saveat=1e-8)
#@time sol4 = solve(prob,KenCarp47(linsolve=KrylovJL_GMRES(),concrete_jac=true),reltol=reltolp,abstol=abstolp,saveat=1e-8)
#@time sol5 = solve(prob,TRBDF2(linsolve=KrylovJL_GMRES(),concrete_jac=true),reltol=reltolp,abstol=abstolp,saveat=1e-8)
@time sol6 = solve(prob,VCABM(),reltol=reltolp,abstol=abstolp,saveat=1e-8)
#@time sol7 = solve(prob,ODEInterfaceDiffEq.radau(),reltol=reltolp,abstol=abstolp,saveat=1e-8)
#@time sol8 = solve(prob,RadauIIA5(),reltol=reltolp,abstol=abstolp,saveat=1e-8)
#reltolp=1e-14
#abstolp=1e-14
#@time sol9 = solve(prob,ODEInterfaceDiffEq.radau(),reltol=reltolp,abstol=abstolp,saveat=1e-8)
#@time sol10 = solve(prob,RadauIIA5(),reltol=reltolp,abstol=abstolp,saveat=1e-8)
#reltolp=1e-6
#abstolp=1e-6
#@time sol8 = solve(prob,Rosenbrock23(autodiff=true),reltol=reltolp,abstol=abstolp,saveat=1e-8)
using Plots
#plot(sol1,vars=(2),yscale=:log10,ylims=(1e-15,1e-2),label="DORMAND_PRINCE_7_4_5 24 sec")
#plot!(sol2,vars=(2),yscale=:log10,ylims=(1e-15,1e-3),label="CVODE_BDF 17.69 sec")
#plot!(sol3,vars=(2),yscale=:log10,ylims=(1e-15,1e-3),label="CVODE_Adams 6.37 sec")
#plot!(sol4,vars=(2),yscale=:log10,ylims=(1e-15,1e-3),label="KenCarp47 24.14 sec")
#plot!(sol5,vars=(2),yscale=:log10,ylims=(1e-15,1e-3),label="TRBDF2 6.75 sec")
#plot(sol7,vars=(2),yscale=:log10,ylims=(1e-15,1e-3),label="ODEInterfaceDiffEq.radau() tol=1e-12")
#plot!(sol8,vars=(2),yscale=:log10,ylims=(1e-15,1e-3),label="RadauIIA5 tol=1e-12")
#plot!(sol9,vars=(2),yscale=:log10,ylims=(1e-15,1e-3),label="ODEInterfaceDiffEq.radau() tol=1e-14")
#plot!(sol10,vars=(2),yscale=:log10,ylims=(1e-15,1e-3),label="RadauIIA5 tol=1e-14")
plot(sol6,vars=(2),yscale=:log10,ylims=(1e-15,1e-3),label="VCABM tol=1e-20")

BigFloat4
BigFloat41242×830 55.3 KB

1 Like
22

And the L-stable methods should flip, what do they look like?

23

One sidenote on BigFloats - this probably won’t affect the integrator, but you should assign BigFloat literals using the big"" string macro, which will correctly preserve the full precision of the input. Passing a Float64 literal to BigFloat limits you to Float64 precision. See:

julia> big(1.0e-2)
0.01000000000000000020816681711721685132943093776702880859375

julia> big"1.0e-2"
0.009999999999999999999999999999999999999999999999999999999999999999999999999999995
24

RadauIIA5, TRBDF2 and KenCarp47 are not working at tol=1e-20. It takes more than 10 minutes to reach 1e-10. So it will take several hours or days for these methods to work with tolerances this tight. RadauIIA5 is giving the correct result without BigFloats and tol=1e-14. I do not expect that to change when tol is changed to 1e-20. As far as TRBDF2 and KenCarp47 they are too slow.

25

They will take awhile, that’s for sure, but it would be good to confirm convergence instead of just assuming it. Try progressively decreasing it, so told=1e-16 then 17 etc. Since all of the methods are convergent, there should be a point at which one of them flips to the behavior of the other, and then there would be reasonable confidence that L-stability is overdamping. If you don’t see that flip though then it’s just speculation. I’m having an internet outage though and can’t do run it from my phone.

4 Likes
26

Don’t worry. I will do it. I will just submit it on HPC instead of running it on my laptop.

2 Likes
27

Here is the result for KenCarp47 and TRBDF2 with BigFloat and tolerance of 1e-20. Both show correct results. I also ran the same code with a tolerance of 1e-18 and it gave wrong results. I am a bit puzzled. I was of the opinion that L-stability is something desirable in an algorithm. What we are seeing here is that L-stability is something that is causing a problem.

KC47andTRBDF

3 Likes
28

Well, it dampens at infinity in order to stabilize the solving. I’ve always known theoretically that it could dampen solutions which have true high frequency oscillations, but this is the first time I’ve seen it in practice.

You can see references to this phenomena all around. For example, In MATLAB’s ODE docs it says to use the equivalent to Trapezoid instead of the equivalent to TRBDF2 in that case:

Use ode23t if the problem is only moderately stiff and you need a solution without numerical damping.

We state it similarly in our docs: https://diffeq.sciml.ai/stable/solvers/ode_solve/

  • Trapezoid - A second order A-stable symmetric ESDIRK method. “Almost symplectic” without numerical dampening. Also known as Crank-Nicolson when applied to PDEs. Adaptive timestepping via divided differences on the memory. Good for highly stiff equations which are non-oscillatory.
  • TRBDF2 - A second order A-B-L-S-stable one-step ESDIRK method. Includes stiffness-robust error estimates for accurate adaptive timestepping, smoothed derivatives for highly stiff and oscillatory problems.

Though before today I couldn’t point to a clear example of where that actually comes into play. Now I can! Cheers and thanks for taking the time to dig into this. This is a fun example of an analysis that really requires the special concoction of features that we’ve put together, so I find this really fun! :sweat_smile:

4 Likes
29

Is there a way to detect whether dampening due to L-stability is happening within a problem in Julia? It can be seen that the solution is unstable by looking at the eigenvalues of the differential equations. Are there techniques that calculate the eigenvalues of the equations to determine whether the solution has real instability or not?

I think I found the winner in the process. For tol of 1e-12 VCABM has the best performance. I will just avoid L-stable methods for the time being.

1 Like
30

Well I will say that this is a very very very rare phenomena. I help people with at least 5-10 different ODE systems a day, for about 6 years now, and this is the first time I’ve seen a reason to reject the L-stable result. So again, theoretically this has been known, but in practice it’s a bit wild to finally see it (the only reason why I could immediately point out the difference was that I was giddy inside that this could finally be the case).

I think the reason why it’s so rare is that people naturally tend to be simulating stable systems, and L-stability helps the simulation of a stable system stay stable. It’s fairly rare for someone to be simulating something that is just going unstable like that, in which case L-stability will dampen errors in such a way that it will make it stable. You can see this in the linear stability analysis for how implicit methods will make u' = a*u stable for some positive a’s, which is “incorrect”, though normally with nonlinear models that “incorrectness” helps stabilize what could be numerical divergences/oscillations in intermediate steps. But of course, that makes it obvious that if your model is truly supposed to diverge then in theory some implicit methods could make them not, which seems to be the case here.

The more common issue with L-stability is probably the effect on high frequency oscillations, where the damping effect can reduce their size. This is seen sometimes in models of circuits, which can bias the spectrum of the result which can matter for Fourier analysis of circuitry. But that’s usually a smaller numerical effect and not as dramatic as this.

Every method has trade-offs, and the trade-off of highly stable methods is that they try really really hard to be stable, which might be too hard in some cases. For your case, indeed I would say VCABM seems to be a good choice.

13 Likes