** 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?
Thanks
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)
push!(xvars,x)
ny,dy,k,d2ydy2 = p
for i in 1:ny
d2ydy2[i] = -2y[i]/dy^2
end
d2ydy2[:] += (circshift(y,1) + circshift(y,-1))/dy^2
dydx[:] = - (im/(2.0*k))*d2ydy2[:]
end
# compute dydx using analytical expression for second derivative along y
function f2!(dydx,y,p,x)
push!(xvars,x)
ny,dy,k,d2ydy2 = p
d2ydy2[:] .= -ky^2*y[:]
dydx[:] .= - (im/(2.0*k))*d2ydy2[:]
end
# 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!"))