I usually come to solve outwards ode like schrodinger equation under columb potential and I recently realize that the initial condition and time span affect the stability of the solution as posted as following:
using DifferentialEquations,Plots,Roots;
pyplot()
function columbode!(du,u,p,r)
v(r)=-1/r
du[1]=u[2]
du[2]=(v(r)-p)*u[1]
end
f(x)=let u0=[0.0,x],rspan=(1e-10,80.0)
prob(e)=ODEProblem(columbode!,u0,rspan,e)
solf(e)=solve(prob(e), Tsit5(), reltol=1e-8, abstol=1e-8)
e_eigen=find_zero(e->solf(e).u[end][1],(-1.0,-0.1))
plot(solf(e_eigen),vars=[1])
end
display(plot(f(1e-3),f(3e-3),f(6e-3),f(1e-2),layout=(2,2)))
savefig("discourse_rescaling_ode.png")
I tested four different initial derivative and found the final solutions perform differently. Only one of them behaves normally. I wonder if there’s way to deal with stability of outwards ode. Currently, I wonder if I can normalize the intermediate sol during integration using callback or integrator interface and I have no good ideas since these operations may corrupt the interpolations stuff. Thanks for any help!