I want to solve numerically a second-order ODE (see first equation below) which depends on relatively easy functions (g, g'/g and f/g are superpositions of cos,sine,cosh,sinh). The main problem is that these functions depend on a function alpha (e.g. cos(alpha*phi)), and only its derivative is known (see second equation below, sigma and d are fixed).
This is obviously difficult to solve, so I tried first something simpler, setting alpha=phi². My code returns the warning Warning: Automatic dt set the starting dt as NaN, causing instability. followed by Warning: NaN dt detected. Likely a NaN value in the state, parameters, or derivative value caused this outcome. The code runs, but plotting the solution gives an error ArgumentError: reducing over an empty collection is not allowed, which I assume is telling me my solution is empty.
I have found a post with a similar issue, in which the error occurred because the parameters were giving a divergent solution, but in my case all the functions in the ODE are well-defined on the interval of interest (-5,5), so I doubt this would be my case. I have three questions then:
- How can I know where is the NaN value exactly?
- How can I avoid this error?
- Is there a way to implement the full expression (d alpha/d phi) into an integrator?
Many thanks!
EDIT: Following the comment of @Lutz Lehmann, I tried to reduce the 2nd order ODE and solve the three equations simultaneously, but Julia tells me UndefVarError: α not defined, even though I explicitly define it as the third equation. What am I doing wrong?
EDIT2: The error disappeared after renaming the functions correctly. However, I receive the same warnings and error as before, with the addition of Warning: First function call produced NaNs. Exiting. in first instance.
EDIT3: I have now written the expression for dα is a simpler way, and corrected two mistakes in the code (sign pointed out by @Lutz Lehmann and a missing closing parenthesis). After giving a small non-zero initial value to α and A', I can see the solution blows up at time ϕ>2.12, so I guess the problem lies there.
EDIT4: I am now closer to the source of the problem, but I don't know how to proceed further. I have tried to see which terms are problematic in the equation for dα/dϕ, and it seems the instability arises for the sinh and cosh terms. For example, trying to solve dα/dϕ = 2d*sinh(s2*α*ϕ) upsets the solver immediately, which returns Warning: Instability detected. Aborting. But Julia should be able to solve this, so what am I missing?
using DifferentialEquations
using Plots
const global lp=1.
const global s2=81.
const global d=-9.0e-4
const global B=1. # Parameter in the gaussian f(\phi) defined below
const global β =1. # Parameter in the gaussian f(\phi) defined below
const global k=1. # Wave number
@variables α,ϕ
# Useful functions
#α(ϕ) = ϕ^2
g1(α,ϕ) = s2^2*α^2*sinh(s2*α*ϕ)*sin(2d*α)
g2(α,ϕ) = cos(2d*α)+cosh(s2*α*ϕ)
g3(α,ϕ) = -α*s2*sin(2d*α)+2d*cos(2d*α)+2d*cosh(s2*α*ϕ)
g(α,ϕ) = exp(-2α)*g3(α,ϕ)/g2(α,ϕ)/2/lp
dgg(α,ϕ) = g1(α,ϕ)/g2(α,ϕ)/g3(α,ϕ) # Derivative of g over g
f1(α,ϕ) = -4B*lp*exp(2α)*ϕ*exp(-ϕ^2/β^2)/β^2
f2(α,ϕ) = cos(2d*α)+cosh(s2*α*ϕ)
f3(α,ϕ) = -α*s2*sin(2d*α)+2d*cos(2d*α)+2d*cosh(s2*α*ϕ)
dfg(α,ϕ) = f1(α,ϕ)*f2(α,ϕ)/f3(α,ϕ) # Derivative of f over g
α1(α,ϕ) = ϕ*s2*sin(2d*α)+2d*sinh(s2*α*ϕ)
function eom(du,u,p,ϕ)
α = u[1]
A = u[2]
dA = u[3]
du[1] = α1(α,ϕ)/g3(α,ϕ)
du[2] = dA
du[3] = -2*dgg(α,ϕ)*dA-2*dfg(α,ϕ)*dA-k^2*A/g(α,ϕ)
end
A₀ = [0.1,0.,0.1] # initial state vector
tspan = (-5.0,5.0) # time interval
prob = ODEProblem(eom,A₀,tspan)
sol = solve(prob,Tsit5())
plot(sol)
# Full solution to implement in the future
# dαdϕ1(ϕ) = ϕ*s2*sin(2d*α)+2d*sinh(s2*α*ϕ)
# dαdϕ(ϕ) = dαdϕ1(ϕ)/g3(ϕ)
# Original problem
# function eom(ddu,du,u,p,ϕ)
# ddu[1] = -2dgg(ϕ)*du[1]-2dfg(ϕ)*du[1]-k^2*u[1]/g(ϕ)
# end
# u0 = [0.] # initial state vector
# du0 = [0.]
# tspan = (-5.0,5.0) # time interval
# prob = SecondOrderODEProblem(eom,du0,u0,tspan)
# sol = solve(prob,Tsit5())
# plot(sol)