I'm attempting to solve a very simple IVP in Python in order to do error analysis:
import numpy as np
from scipy.integrate import solve_ivp
def dh_dt_zerodriver(t, h):
return -2 / h
t = 50
steps = 10
dt = t / steps
t_span = [0, t]
t_eval = np.arange(0, t + dt, dt)
h0 = 5
sol = solve_ivp(dh_dt_zerodriver, t_span=t_span, y0=[h0], t_eval=t_eval)
However, the solution will not compute as it runs for an indefinite amount of time. I have used solve_ivp to solve more complex IVPs without analytical solutions, but this one seems to have me stumped despite it seemingly being a fairly simply problem.
My paper is based on using the RK45 method, but if I must use some different numerical method then it is okay.
The problem is the domain you want your solution. Your ODE is h'(t) = -2/h(t) where h(0) = 5. The analytical solution is h(t) = sqrt(25 - 4t). For t < 6.25, the solution is real, but for t >= 6.25, the solution becomes complex. If you were to instead integrate over t = [0, 6], you would quickly get the correct result.