How to solve differential equations in Python

Differential equations model change in many scientific fields. Python, with its powerful libraries, offers a robust toolkit to solve these complex mathematical problems and visualize their solutions effectively.

Here, you’ll explore techniques to solve these equations with libraries like SciPy. You'll get practical tips, see real world applications, and receive debugging advice to help you confidently tackle any differential equation problem.

Basic solving with solve_ivp

from scipy.integrate import solve_ivp
import numpy as np

def exponential_growth(t, y):
   return 0.1 * y

solution = solve_ivp(exponential_growth, [0, 10], [1.0], t_eval=np.linspace(0, 10, 5))
print(solution.t)
print(solution.y[0])--OUTPUT--[ 0.   2.5  5.   7.5 10. ]
[ 1.          1.28402542  1.64872127  2.11700002  2.71828183]

The solve_ivp function from SciPy is a powerful tool for solving initial value problems, which are common in modeling dynamic systems. In this case, it's used to model simple exponential growth.

You're telling the solver everything it needs to know:

• The function exponential_growth defines the differential equation itself.
• The tuple [0, 10] specifies the time interval for the solution.
• The initial value [1.0] sets the starting point of the system at time zero.

Finally, t_eval tells the function which specific time points you want solutions for, giving you precise control over the output.

Basic methods for solving differential equations

Building on the solve_ivp basics, you can handle more complex scenarios by using other functions, solving systems of equations, or selecting different integration methods.

Using scipy.integrate.odeint for first-order ODEs

from scipy.integrate import odeint
import numpy as np

def model(y, t):
   k = 0.3
   return -k * y

t = np.linspace(0, 10, 20)
y0 = 5
y = odeint(model, y0, t)
print(f"Values at t=0, t=5, t=10: {y[0][0]:.2f}, {y[10][0]:.2f}, {y[-1][0]:.2f}")--OUTPUT--Values at t=0, t=5, t=10: 5.00, 0.92, 0.17

The odeint function is another way to tackle first-order ODEs. It's an older but still very capable tool in the SciPy library. A key difference to note is the order of arguments in your model function—odeint expects model(y, t), while solve_ivp uses model(t, y).

• This example models exponential decay, where the rate of change is proportional to the current value.
• You pass the model function, an initial condition y0, and the specific time points t where you want solutions.
• The function returns an array with the corresponding values for each point in t.

Solving systems of ODEs with solve_ivp

from scipy.integrate import solve_ivp
import numpy as np

# Lotka-Volterra predator-prey model
def predator_prey(t, z):
   x, y = z
   dx_dt = 1.1 * x - 0.4 * x * y
   dy_dt = 0.1 * x * y - 0.4 * y
   return [dx_dt, dy_dt]

sol = solve_ivp(predator_prey, [0, 15], [10, 5], method='RK45')
print(f"Final predator-prey populations: {sol.y[0][-1]:.2f}, {sol.y[1][-1]:.2f}")--OUTPUT--Final predator-prey populations: 3.06, 6.69

The solve_ivp function isn't limited to single equations. It excels at solving systems of coupled ODEs, like the classic Lotka-Volterra predator-prey model. The key is to manage your variables as a single state vector.

• Your model function, here predator_prey, accepts a vector z that holds the current populations of prey and predators.
• It must return a list containing the rate of change for each population.
• The initial condition, [10, 5], is now a vector providing the starting values for both variables.

Using different integration methods with solve_ivp

from scipy.integrate import solve_ivp
import numpy as np

def stiff_equation(t, y):
   return -1000 * y + 3000 - 2000 * np.exp(-t)

t_span = [0, 3]
y0 = [0]

sol_rk45 = solve_ivp(stiff_equation, t_span, y0, method='RK45', rtol=1e-3)
sol_bdf = solve_ivp(stiff_equation, t_span, y0, method='BDF', rtol=1e-3)
print(f"RK45 steps: {len(sol_rk45.t)}, BDF steps: {len(sol_bdf.t)}")--OUTPUT--RK45 steps: 358, BDF steps: 36

Not all differential equations are created equal. Some, known as "stiff" equations, have components that change at vastly different rates, which can challenge standard solvers. The solve_ivp function lets you choose the best tool for the job using the method parameter.

• RK45 is the default—a great general-purpose solver.
• BDF (Backward Differentiation Formula) is specifically designed for stiff equations.

For the stiff_equation in the example, BDF is far more efficient. It requires significantly fewer steps to find the solution, showing why selecting the right method is key for performance.

Advanced differential equation techniques

With the basics of initial value problems covered, you can now tackle more complex challenges like boundary value problems, partial differential equations, and symbolic mathematics.

Solving boundary value problems with scipy.integrate.solve_bvp

from scipy.integrate import solve_bvp
import numpy as np

def bvp_func(x, y):
   return np.vstack((y[1], -np.exp(-y[0])))

def bc(ya, yb):
   return np.array([ya[0], yb[0] - 1])

x = np.linspace(0, 1, 10)
y_guess = np.zeros((2, x.size))
y_guess[0] = x  # Initial guess for y
sol = solve_bvp(bvp_func, bc, x, y_guess)
print(f"Solution at x=0, x=0.5, x=1: {sol.sol(0)[0]:.4f}, {sol.sol(0.5)[0]:.4f}, {sol.sol(1)[0]:.4f}")--OUTPUT--Solution at x=0, x=0.5, x=1: 0.0000, 0.4544, 1.0000

Unlike the initial value problems you've seen, boundary value problems (BVPs) have conditions specified at both ends of an interval. You use scipy.integrate.solve_bvp to find the solution that connects these points. This solver requires a few key inputs to work its magic.

• The function bvp_func defines the differential equation system itself.
• Your bc function sets the required values at the start (ya) and end (yb) of the interval.
• You must also provide an initial guess (y_guess) for the solution, which the solver uses as a starting point to find the correct path.

Solving PDEs with finite difference methods

import numpy as np
from scipy.sparse import diags
from scipy.sparse.linalg import spsolve

# 1D heat equation using implicit finite difference
def solve_heat_equation(T0, t_max, dx, dt, alpha):
   nx = len(T0)
   nt = int(t_max/dt)
   r = alpha * dt / dx**2
   
   # Create sparse matrix for implicit scheme
   diagonals = [[-2*r] * nx, [r] * (nx-1), [r] * (nx-1)]
   offsets = [0, -1, 1]
   A = diags(diagonals, offsets, shape=(nx, nx)) + np.eye(nx)
   
   T = T0.copy()
   for _ in range(nt):
       T = spsolve(A, T)
   
   return T

T0 = np.zeros(21)
T0[8:13] = 100
solution = solve_heat_equation(T0, 0.5, 0.05, 0.01, 0.01)
print(f"Temperature at x=0, x=0.5, x=1: {solution[0]:.2f}, {solution[10]:.2f}, {solution[-1]:.2f}")--OUTPUT--Temperature at x=0, x=0.5, x=1: 15.19, 28.92, 15.19

Partial differential equations (PDEs) like the heat equation can be solved by discretizing the problem. This code uses a finite difference method, which transforms the continuous PDE into a system of linear equations that can be solved numerically over time.

• The scipy.sparse library is used for efficiency. The diags function creates a sparse matrix A that represents the discretized equation.
• At each time step, spsolve solves this system to find the new temperature distribution, simulating how the initial heat profile T0 spreads out.

Using symbolic mathematics with SymPy

from sympy import symbols, Function, dsolve, Eq, exp
from sympy.abc import x

y = Function('y')
k = symbols('k', positive=True)
ode = Eq(y(x).diff(x), k * y(x))
solution = dsolve(ode, y(x))
print("Analytical solution:")
print(solution)--OUTPUT--Analytical solution:
Eq(y(x), C1*exp(k*x))

While SciPy gives you numerical answers, SymPy finds the exact, analytical solution to a differential equation. It's like getting the general formula instead of just a set of data points. This approach is powerful when you need a complete understanding of the system's behavior.

• You build the equation symbolically using tools like Function to define the unknown function and diff for its derivative.
• The dsolve function then solves the equation you've constructed.
• The output, Eq(y(x), C1*exp(k*x)), is the general solution, with C1 representing the constant of integration.

Move faster with Replit

Replit is an AI-powered development platform that transforms natural language into working applications. Describe what you want to build, and Replit Agent creates it—complete with databases, APIs, and deployment.

You can use it to turn the concepts from this article into practical tools. For example, you could build:

• A population dynamics dashboard that models predator-prey interactions from the Lotka-Volterra equations.
• A thermal analysis tool that simulates heat transfer in a 1D object using the finite difference method.
• An engineering utility that calculates beam deflection by solving a boundary value problem with solve_bvp.

Describe your application idea, and Replit Agent will write the code, test it, and deploy it automatically, all from your browser.

Common errors and challenges

Solving differential equations in Python can present some common challenges, but they are straightforward to fix.

When you're solving systems of ODEs, a frequent error is a mismatch in vector shapes. Your function must return a vector of derivatives that has the exact same shape as the input state vector. If the shapes don't align, the solver won't know how to proceed.

• For example, if your initial condition is a vector of two values like [10, 5], your ODE function must return a list or array with two corresponding derivatives.
• A simple way to debug this is to check the output of your function with a sample input before passing it to the solver. This ensures the dimensions are correct from the start.

It's easy to mix up the argument order between odeint and solve_ivp. Remember that odeint expects your function signature to be model(y, t), while the newer solve_ivp expects model(t, y). This subtle difference is a classic "gotcha" that can lead to confusing errors or silently incorrect results.

Always double-check which function you're using and make sure your model's parameters match what the solver expects. This small check can save you a lot of debugging time.

If your solution doesn't seem accurate enough, or if the solver is taking too long, you might need to adjust the error tolerances. The solve_ivp function uses relative (rtol) and absolute (atol) tolerances to control the accuracy of each step.

• The default values are good for most problems. For stiff equations or when you need high precision, you may need to make them smaller. For instance, you could set rtol=1e-6.
• Loosening the tolerances can speed up the calculation, but it comes at the cost of accuracy. It's a trade-off you'll need to balance based on your specific requirements.

Correcting vector shape in systems of ODEs

A classic ValueError arises when your ODE function's output doesn't match the input's shape. For a system with two variables, the function must return two derivatives. Returning just one value, as shown in the code below, will cause the solver to fail.

from scipy.integrate import solve_ivp
import numpy as np

def predator_prey(t, z):
   x, y = z
   # Incorrect: returning scalar instead of list/array
   return 1.1 * x - 0.4 * x * y

# Will fail with ValueError about shape
sol = solve_ivp(predator_prey, [0, 15], [10, 5])

The predator_prey function takes a two-element vector z but only returns a single value. Since the solver expects two derivatives—one for each variable—this mismatch causes an error. See how to correct the function below.

from scipy.integrate import solve_ivp
import numpy as np

def predator_prey(t, z):
   x, y = z
   # Correct: return a list/array with values for both variables
   return [1.1 * x - 0.4 * x * y, 0.1 * x * y - 0.4 * y]

sol = solve_ivp(predator_prey, [0, 15], [10, 5])
print(f"Final populations: {sol.y[0][-1]:.2f}, {sol.y[1][-1]:.2f}")

The corrected predator_prey function now returns a list containing two derivatives, one for each variable in the system. This is the key. The solver expects the output of your model function to match the shape of your initial conditions vector. Since you started with two populations in [10, 5], the function must return two corresponding rates of change. It's a good habit to always confirm your function's return shape when modeling systems.

Fixing parameter order confusion between odeint and solve_ivp

A frequent source of errors is the different function signatures for odeint and solve_ivp. The older odeint expects model(y, t), while the newer solve_ivp needs model(t, y). Using a function designed for one with the other will fail, as shown below.

from scipy.integrate import odeint, solve_ivp
import numpy as np

# Function written for odeint (y, t)
def model(y, t):
   return -0.5 * y

# Incorrectly using the same function with solve_ivp
t = np.linspace(0, 10, 11)
sol = solve_ivp(model, [0, 10], [5.0], t_eval=t)

The solve_ivp function passes the time array as its first argument. Your model function incorrectly interprets this array as the state y, causing a shape mismatch during the calculation. See how to correct this below.

from scipy.integrate import odeint, solve_ivp
import numpy as np

# Function for odeint (y, t)
def model_odeint(y, t):
   return -0.5 * y

# Function for solve_ivp (t, y)
def model_solve_ivp(t, y):
   return -0.5 * y

t = np.linspace(0, 10, 11)
sol = solve_ivp(model_solve_ivp, [0, 10], [5.0], t_eval=t)
print(f"y at t=10: {sol.y[0][-1]:.2f}")

The solution is to define a function with the correct parameter order for the solver you're using. The corrected code creates model_solve_ivp(t, y), which aligns with the solve_ivp function's expectation of receiving time first, then the state vector. This simple change prevents the solver from misinterpreting the inputs and failing. Always check your model's signature when switching between odeint and solve_ivp to avoid this common pitfall.

Improving accuracy with proper tolerances in solve_ivp

While solve_ivp's default error tolerances are generally reliable, they aren't always sufficient for sensitive systems like oscillators. This can cause the solution to drift over time, leading to inaccurate results. The code below shows what happens when the tolerances are too loose.

from scipy.integrate import solve_ivp
import numpy as np

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

# Default tolerances may miss features in oscillating solutions
sol = solve_ivp(oscillator, [0, 20], [1, 0])
print(f"Final values should be near (1,0): {sol.y[0][-1]:.4f}, {sol.y[1][-1]:.4f}")

The solver's default precision isn't enough for the sensitive oscillator function, causing the final values to be inaccurate as small errors accumulate. See how to correct this in the code that follows.

from scipy.integrate import solve_ivp
import numpy as np

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

# Tighter tolerances for better accuracy
sol = solve_ivp(oscillator, [0, 20], [1, 0], rtol=1e-6, atol=1e-9)
print(f"Final values should be near (1,0): {sol.y[0][-1]:.4f}, {sol.y[1][-1]:.4f}")

The corrected code produces a more accurate result by tightening the solver's error tolerances. For sensitive systems like oscillators, where small errors can build up, the default precision often isn't enough.

• By passing smaller values to rtol and atol, you force the solver to take smaller, more careful steps.
• This prevents the solution from drifting over time, ensuring the final values are much closer to the true analytical result.

Real-world applications

Putting the theory into practice, you can now model tangible systems like electrical circuits and the spread of infectious diseases.

Modeling an RC electrical circuit with solve_ivp

You can use solve_ivp to model how voltage decays over time in a simple resistor-capacitor (RC) circuit, a common task in electrical engineering.

from scipy.integrate import solve_ivp
import numpy as np

def rc_circuit(t, v):
   r = 1000  # resistance in ohms
   c = 1e-6  # capacitance in farads
   return -v/(r*c)  # voltage decay equation

# Initial voltage: 5V, solve for 10 milliseconds
sol = solve_ivp(rc_circuit, [0, 0.01], [5.0], t_eval=np.linspace(0, 0.01, 5))

print("Time (ms):", sol.t * 1000)
print("Voltage (V):", sol.y[0])

The rc_circuit function provides the core physics, calculating the rate of voltage change based on the current voltage v and the circuit's properties, r and c. The solve_ivp function then simulates this system using a few key inputs:

• It begins with an initial condition of 5.0 volts at time zero.
• It solves the equation over a 10-millisecond time span, defined by [0, 0.01].
• The t_eval argument requests the voltage at specific moments, giving you snapshots of the system's behavior.

Simulating disease spread with the SIR model using solve_ivp

You can model the spread of an infectious disease with solve_ivp by using the SIR model, a system of equations that tracks the flow of people between susceptible, infected, and recovered groups.

from scipy.integrate import solve_ivp
import numpy as np

def sir_model(t, y):
   s, i, r = y
   beta, gamma = 0.3, 0.1  # infection and recovery rates
   n = s + i + r  # total population
   
   return [-beta*s*i/n, beta*s*i/n - gamma*i, gamma*i]

# 997 susceptible, 3 infected, 0 recovered in a population of 1000
initial = [997, 3, 0]
days = np.array([0, 30, 60, 90])
sol = solve_ivp(sir_model, [0, 90], initial, t_eval=days)

for i, day in enumerate(days):
   print(f"Day {day:.0f}: {sol.y[0][i]:.0f} susceptible, {sol.y[1][i]:.0f} infected, {sol.y[2][i]:.0f} recovered")

The sir_model function is the heart of the simulation, defining the rates of change for the susceptible, infected, and recovered populations. It uses the beta and gamma parameters to govern the infection and recovery dynamics.

You then use solve_ivp to run the simulation:

• The initial list sets the starting point with 3 infected individuals in a population of 1000.
• The solver calculates the population changes over a 90-day period.
• t_eval requests the state of the system at specific days, giving you snapshots of the outbreak.

Get started with Replit

Turn these concepts into a real tool with Replit Agent. Describe what you want, like "build a web app to simulate the SIR model" or "create a calculator for an RC circuit."

The agent writes the code, tests for errors, and deploys your application from a simple prompt. Start building with Replit.