Ordinary Differential Equations in 1D

Initial Value Problems and Nonlinear Rate Laws

This first session builds the numerical workflow used throughout the whole course: define a right-hand side, choose an initial condition, integrate forward in time, and interpret the solution (Strogatz 2024; Virtanen et al. 2020).

From One Variable to Small Systems

The core mathematical object is the initial value problem

\[ \dot y = f(t, y), \qquad y(t_0) = y_0. \tag{1}\]

In a scalar problem, \(y(t)\) is one number. In a small system, \(y(t)\) becomes a vector. The numerical workflow stays the same.

That is why this session starts with the label “1D” but already includes the SIR model: the main skill is not the number of components, but how you turn a rate law into a computable trajectory.

Case Studies

Start with a classical epidemic model and see how a coupled ODE system fits into the same solve_ivp pipeline.

SIR Epidemic Model

Then study a saturating nonlinear rate law through a Michaelis–Menten inspired example.

Michaelis–Menten Kinetics

Finally, combine phase-line reasoning, numerical integration, and interactivity in the spruce budworm model.

Spruce Budworm Model

When the workflow is clear, move to the assignment.

Assignment

Learning Goals

  • Formulate an ODE or a small ODE system in Python.
  • Use scipy.integrate.solve_ivp() to integrate an initial value problem.
  • Interpret trajectories in time and connect them to parameters.
  • Recognize threshold, saturation, and bistability effects in simple models.
  • Reuse the same solver workflow across very different applications.

Flow of This Session

  • Write the right-hand side of the model.
  • Choose parameters, a time span, and an initial condition.
  • Integrate with solve_ivp().
  • Plot the solution and interpret the result.
  • Compare how the behavior changes when you vary parameters or initial data.

What do we need?

scipy.integrate.solve_ivp

This is the main solver for initial value problems in the course. It takes the right-hand side fun(t, y), the interval t_span, the initial state y0, and optional arguments such as t_eval, method, and args.

import numpy as np
from scipy.integrate import solve_ivp


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


sol = solve_ivp(
    exponential_decay,
    t_span=(0, 10),
    y0=[2.0],
    t_eval=np.linspace(0, 10, 200),
)

numpy

Use NumPy arrays for initial conditions, evaluation grids, and post-processing of trajectories.

matplotlib.pyplot

Use Matplotlib to visualize time series, parameter comparisons, and phase-line style diagnostics.

Why start with solve_ivp()?

Because once you understand the initial value problem pipeline, you can reuse it almost unchanged in later sessions on planar systems, coupled oscillators, and PDE-inspired reductions.

Tip

The default RK45 method is a good first choice for most examples in this session. Later sessions will make the method choice itself part of the modeling discussion.

References

Strogatz, Steven H. 2024. Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry, and Engineering. Chapman; Hall/CRC.
Virtanen, Pauli, Ralf Gommers, Travis E. Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau, Evgeni Burovski, et al. 2020. “SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python.” Nature Methods 17 (3): 261–72. https://doi.org/10.1038/s41592-019-0686-2.