Unexpected behavior DifferentialEquations (parabolic ODE) -- large # of steps

** This question is as much about numerics as about use of DifferentialEquations.jl **
** This is a toy problem with an analytical solution – I’d like to understand the behavior of DifferentialEquations.jl, not interested in the solution **

I solve a simple parabolic ODE dy/dx = - im/(2 k) d^2(y)/dx^2 between x=0 and x=xspan s.t. periodic boundary conditions along y. The excitation (yin) is periodic and so is the solution. I solve the problem using

prob = ODEProblem(f!,yin,xspan,p)
sol = solve(prob,Vern7(),abstol=myabstol,reltol=myreltol,save_everystep=true)

for two versions of f:

  • f1 computes d^2(y)/dx^2 via a second order finite difference equation
  • f2 computes it analytically from knowledge of the periodicity along y

Both yield the expected solution (the solutions are in essence identical and only different because of the small difference between the approximate wavenumber used if f1 vs the exact in f2).

But f2 requires 6000 steps, whereas f1 only needs 45 (which is more what I expected). The code below plots all relevant data.

How can this behavior be explained, and how can I call the ODE solver when using f1 to avoid this many steps?


using DifferentialEquations

ny = 1000    #  of points along y
dy = 0.1     # discretization step along y
xspan = 90.0   # size of domain along x
nhar = 2;      # number of periods along y

k = 2pi  # wavenumber
ky = nhar*k/(ny*dy)  # wavenumber along x
kx = sqrt(k^2-ky^2)   # wavenumber along y
phi = atan(ky/kx)     # angle of propagation

yin = [exp(-im*k*j*dy*sin(phi)) for j in 1:ny]  # periodic excitation

myabstol = 1e-3
myreltol = 1e-3

# compute dydx using central second difference derivative along y
function f1!(dydx,y,p,x) 
    ny,dy,k,d2ydy2 = p
    for i in 1:ny
      d2ydy2[i] = -2y[i]/dy^2
    d2ydy2[:] += (circshift(y,1) + circshift(y,-1))/dy^2
    dydx[:] = - (im/(2.0*k))*d2ydy2[:]

# compute dydx using analytical expression for second derivative along y
function f2!(dydx,y,p,x)
    ny,dy,k,d2ydy2 = p
    d2ydy2[:] .= -ky^2*y[:]
    dydx[:] .= - (im/(2.0*k))*d2ydy2[:]

# solve problem using f1, plot x vs step and compare to analytical solution
d2ydy2 = similar(fin)
p = (ny,dy,k,d2ydy2)
xvars = [];
prob1 = ODEProblem(f1!,yin,xspan,p)
sol1 = solve(prob1,Vern7(),abstol=myabstol,reltol=myreltol,save_everystep=true)
display(Plots.plot(xvars, title = "method 1: 6000 steps!"))
# compute ky resulting from numerical discretization
kyapp = sqrt(-(exp(im*ky*dy)-2.0+exp(-im*ky*dy))/dy^2)
fouta = [exp(-im*ky*j*dy+im*kyapp^2/(2*k)*xeval) for j in 1:ny]*exp(im*k*xeval)
display(Plots.plot(abs.(sol1(xspan)-fouta),title = "method 1: accurate!"))

# solve problem using f2, plot x vs step and compare to analytical solution
xvars = [];
prob2 = ODEProblem(f2!,yin,xspan,p)
sol2 = solve(prob2,Vern7(),abstol=myabstol,reltol=myreltol,save_everystep=true)
display(Plots.plot(xvars, title = "45 steps!"))
fouta = [exp(-im*ky*j*dy+im*ky^2/(2*k)*xeval) for j in 1:ny]*exp(im*k*xeval)
display(Plots.plot(abs.(sol2(xspan)-fouta),title = "method 2: accurate!"))

The first way introduces a stability criteria and thus a stepsize limit.

Thank you.

I am having some trouble getting other solver options to work. For example in the above code):

sol1 = solve(prob1,AutoTsit5(Rosenbrock23()),abstol=myabstol,reltol=myreltol,save_everystep=false)

=> MethodError: no method matching ForwardDiff.DerivativeConfig(::DiffEqDiffTools.TimeGradientWrapper{ODEFunction{true,typeof(f1!),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Array{Complex{Float64},1},Tuple{Int64,Float64,Float64,Array{Complex{Float64},1}}}, ::Array{Complex{Float64},1}, ::Float64)
Closest candidates are:
ForwardDiff.DerivativeConfig(::F, !Matched::AbstractArray{Y<:Real,N} where N, ::X<:Real) where {F, …

Am I correctly interpreting that I should rewrite the f1 function to act on, and return, Real quantities?
Any suggestions for solvers that operate directly on Complex quantities?

Just use AutoTsit5(Rosenbrock23(autodiff=false)) because of the AD issues with complex numbers.