# Inaccurate results of integration using scipy solve_ivp

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)


• Leaving t_eval out you get the discretization that is actually used internally. I suspect that it is rather sparse with the default tolerances, with nodes for $z>1$ corresponding to the extrema of the error function, either close to the extrema themselves or to the midpoints between them. Jul 26, 2022 at 14:25

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.

• Thanks for this answer! That fixes my problem. I had not realised that the default tolerances were this "large", and in retrospect should have checked that. Many thanks for the help. Jul 26, 2022 at 11:06

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)

# 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.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.ylabel(r'$$\epsilon_{relative}$$', fontsize=18)

• You only need to write out the optional parameters that you actually change from their default value, like here atol and rtol. This makes for an easier reading of that call. Dec 4, 2022 at 20:04