Inaccurate results of integration using scipy solve_ivp

Question

I am trying to use solve_ivp to solve the following 1st order ODE:

$$ \frac{d \rho}{d z} = \frac{m \theta}{(1+\theta z)} \, \rho, $$

subject to $\rho(z=0)=1$, where $m$ and $\theta$ are constants. This ODE actually has a simple analytical solution: $$ \rho = (1+\theta z)^m. $$

When I try to solve this ODE numerically using solve_ivp, the resulting solution does not seem to accurately match the analytical solution. Please see the minimal working example below and the plot of the difference between the numerical and analytical solution.

Any help figuring out what is going on in here would be greatly appreciated. Thanks.

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

theta = 2.5
m = 1.4

def drho_func(z, rho):
    drho = rho*(m*theta)/(1.0+theta*z)
    return drho

zgrid = np.linspace(0,1,501)
y0 = [1]

sol = solve_ivp(drho_func , [zgrid[0], zgrid[-1]], y0, t_eval=zgrid)

rho_num = sol.y[0]
rho_true = (1+theta*zgrid)**m

plt.plot(rho_num-rho_true,zgrid)

Answer 1

The default tolerances for solve_ivp are rtol=1e-3 and atol=1e-6. The analytical value of $\rho$ varies between 1 and 6, so you can expect your numerical error to be somewhere on the order of $\epsilon \approx 0.001$. The fact that you actually see $\epsilon \approx 0.0002$ is not unexpected.

If you want to increase the accuracy of solve_ivp, simply specify a tighter numerical tolerance. For example, here I set rtol=1e-8 and atol=1e-8, and the difference between analytical and numerical solutions is on the order of 1e-8.

Answer 2

As the initial value problem was given and very well specified, the analytical solution can be easily calculated, where $y(t) = (1+2.5t)^{1.4}$ solution of the PVI problem, by the way, $\rho = y$ and $z = t$. With the initial value problem defined, and with the analytical solution in hand, we can solve the ODE numerically as an integrator scipy.integrate.odeint and plot the errors committed $\epsilon_\text{Absolute}(n) = |\rho_\text{Analytic}-\rho_\text{Numeric}|$ and $\epsilon_\text{Relative}(n) = |\frac{\rho_\text{Analytic}-\rho_\text{Numeric}}{\epsilon_\text{Analytic}}|$, where $n$ is the number of iterations.

Observation: I gave preference to odeint, because in several examples like the one specified above, odeint is chosen over solve_ivp.

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint, solve_ivp

# NUMERICAL SOLUTION APPLYING SCIPY'S ODEIN
# ORDINARY DIFFERENTIAL EQUATION dy/dt = f(t,y), t0 = t' and y(t0) = y0
def f(y,t, m, theta):
  # rho => y
  # z => t
  # m = 1.4
  # theta = 2.5
  return m*theta*y/(1+theta*t)

ti = 0.0
tf = 50.0
N  = 1000
t  = np.linspace(ti,tf,N)
y0 = 1.0

sol = odeint(f, y0, t, args=(1.4,2.5), rtol=1E-8, atol=1E-8)

sol_numeric = sol[:,0]

# ANALYTICAL SOLUTION OF THE INITIAL VALUE PROBLEM
def f_solution(t):
  return 1.0*(1+2.5*t)**(1.4)

sol_analytic = f_solution(t) # Observation, t  = np.linspace(ti,tf,N)

# ERROR
erro_absolut = np.abs(sol_analytic-sol_numeric)
erro_relativ = np.abs((sol_analytic-sol_numeric)/sol_analytic)
quantidade_iteracoes = np.linspace(0,N,N)

# PLOTTING OF RESULTS
plt.style.use('dark_background')
plt.figure(figsize = (21,6))
plt.plot(t,sol,'b.',t,f_solution(t),'r-')
plt.title('Comparison between analytical and numerical solution (Odeint)')
plt.xlabel('time')
plt.ylabel('y(t)')
plt.grid(lw = 1.0,color = 'y',linestyle = '-')
plt.show()

plt.figure(figsize = (21,6))
plt.style.use('dark_background')
plt.subplot(1,2,1)
#plt.plot(quantidade_iteracoes,erro_absolut,color = 'red', lw = 3.0)
plt.plot(quantidade_iteracoes,erro_absolut,'r.', lw = 3.0)
plt.title('Absolute error', fontsize=18)
plt.xlabel('num of iteration', fontsize=18)
plt.ylabel(r'$\epsilon_{Absolute}$', fontsize=18)
plt.grid(lw = 0.5,color = 'y',linestyle = '-')

plt.subplot(1,2,2)
#plt.plot(quantidade_iteracoes,erro_relativ,color = 'blue', lw = 3.0)
plt.plot(quantidade_iteracoes,erro_relativ,'b.', lw = 3.0)
plt.title('Relative error', fontsize=18)
plt.xlabel('num of iteration', fontsize=18)
plt.ylabel(r'$\epsilon_{relative}$', fontsize=18)
plt.grid(lw = 0.5,color = 'y',linestyle = '-')

# https://danielmuellerkomorowska.com/2021/02/16/differential-equations-with-scipy-odeint-or-solve_ivp/#:~:text=The%20primary%20advantage%20is%20that,UPDATE%2007.02.