Solving ODEs (solve_ivp) in SciPy

Introduction

We use solve_ivp to find solutions to ordinary differential equations (ODEs). It helps us understand how things change over time.

Modeling how a population grows or shrinks over time.

Predicting the speed and position of a moving object.

Simulating how a chemical reaction changes with time.

Understanding how temperature changes in a system.

Studying the spread of diseases over days or weeks.

Syntax

SciPy

from scipy.integrate import solve_ivp

solution = solve_ivp(fun, t_span, y0, method='RK45', t_eval=None)

# fun: function defining the ODE system
# t_span: tuple with start and end times (t0, tf)
# y0: initial values of the variables
# method: integration method (default 'RK45')
# t_eval: times at which to store the solution

The fun function must take time t and current values y and return their derivatives.

If you want results at specific times, use t_eval to list those times.

Examples

This solves the simple decay ODE dy/dt = -0.5y from time 0 to 10 with initial value 5.

SciPy

def fun(t, y):
    return -0.5 * y

solution = solve_ivp(fun, (0, 10), [5])

This solves a system representing a simple harmonic oscillator with initial position 0 and velocity 1.

SciPy

def fun(t, y):
    return [y[1], -y[0]]

solution = solve_ivp(fun, (0, 6.28), [0, 1], t_eval=[0, 1, 2, 3, 4, 5, 6])

Sample Program

This program solves the exponential decay ODE dy/dt = -0.5y from time 0 to 10 starting at 5. It prints values at intervals and shows a plot of the decay curve.

SciPy

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

def exponential_decay(t, y):
    return -0.5 * y

# Solve from t=0 to t=10 with initial value y=5
solution = solve_ivp(exponential_decay, (0, 10), [5], t_eval=np.linspace(0, 10, 100))

# Print some results
for t, y in zip(solution.t[::20], solution.y[0][::20]):
    print(f"Time {t:.1f}, Value {y:.3f}")

# Plot the solution
plt.plot(solution.t, solution.y[0])
plt.title('Exponential Decay')
plt.xlabel('Time')
plt.ylabel('Value')
plt.grid(True)
plt.show()

Important Notes

Always check that your fun function returns derivatives as a list or array.

Use t_eval to get solution values at specific times for easier analysis or plotting.

Different methods like 'RK45' or 'BDF' can be chosen depending on the problem type.

Summary

solve_ivp helps solve ordinary differential equations by computing values over time.

You define the change rules in a function and give starting conditions and time range.

It returns results you can print, analyze, or plot to understand how the system behaves.