|
1 | 1 | --- |
2 | 2 | title: "Ordinary Differential Equations in 1D" |
3 | | -subtitle: "Numerical Integration with SciPy" |
| 3 | +subtitle: "Initial Value Problems and Nonlinear Rate Laws" |
4 | 4 | format: html |
5 | 5 | --- |
6 | 6 |
|
7 | | -Numerical integration is fundamental for solving ordinary differential equations (ODEs) that don't have analytical solutions [@strogatz2024nonlinear; @virtanen2020scipy]. In this first session, you will learn how to: |
| 7 | +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 [@strogatz2024nonlinear; @virtanen2020scipy]. |
8 | 8 |
|
9 | | -- Formulate ODEs in Python. |
10 | | -- Use SciPy's `solve_ivp` to numerically integrate ODEs. |
11 | | -- Visualize solutions and explore parameter spaces. |
12 | | -- Apply these techniques to real-world models. |
| 9 | +## From One Variable to Small Systems |
13 | 10 |
|
| 11 | +The core mathematical object is the initial value problem |
14 | 12 |
|
15 | | -## Contents |
| 13 | +$$ |
| 14 | +\dot y = f(t, y), |
| 15 | +\qquad |
| 16 | +y(t_0) = y_0. |
| 17 | +$$ {#eq-ode-1d-index-1} |
16 | 18 |
|
17 | | -- [SIR Epidemic Model](sir.qmd) |
18 | | -- [Michaelis–Menten Enzyme Kinetics](michaelis-menten.qmd) |
19 | | -- [Spruce Budworm Population Model](spruce-budworm.qmd) |
20 | | -- [Assignment](assignment.qmd) |
| 19 | +In a scalar problem, $y(t)$ is one number. In a small system, $y(t)$ becomes a vector. The numerical workflow stays the same. |
21 | 20 |
|
22 | | -## The Initial Value Problem |
| 21 | +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. |
23 | 22 |
|
24 | | -An **initial value problem (IVP)** consists of: |
| 23 | +## Case Studies |
25 | 24 |
|
26 | | -$$\frac{dy}{dt} = f(t, y), \quad y(t_0) = y_0$$ {#eq-ode-1d-index-1} |
| 25 | +Start with a classical epidemic model and see how a coupled ODE system fits into the same `solve_ivp` pipeline. |
27 | 26 |
|
28 | | -Where: |
29 | | -- $f(t, y)$ is the rate of change function |
30 | | -- $y_0$ is the initial condition at time $t_0$ |
| 27 | +[SIR Epidemic Model](sir.qmd){.btn .btn-primary} |
31 | 28 |
|
32 | | -## SciPy's `solve_ivp` |
| 29 | +Then study a saturating nonlinear rate law through a Michaelis–Menten inspired example. |
33 | 30 |
|
34 | | -The `scipy.integrate.solve_ivp` function is the standard tool for solving ODEs in Python. Try the following code: |
| 31 | +[Michaelis–Menten Kinetics](michaelis-menten.qmd){.btn .btn-primary} |
35 | 32 |
|
36 | | -```python |
37 | | -import numpy as np |
38 | | -from scipy.integrate import solve_ivp |
| 33 | +Finally, combine phase-line reasoning, numerical integration, and interactivity in the spruce budworm model. |
39 | 34 |
|
40 | | -def exponential_decay(t, y): |
41 | | - return -0.5 * y |
| 35 | +[Spruce Budworm Model](spruce-budworm.qmd){.btn .btn-primary} |
42 | 36 |
|
43 | | -t_span = [0, 10] |
44 | | -y0 = [2, 4, 8] |
45 | | -sol = solve_ivp( |
46 | | - fun=exponential_decay, |
47 | | - t_span=t_span, |
48 | | - y0=y0) |
| 37 | +When the workflow is clear, move to the assignment. |
49 | 38 |
|
50 | | -print(sol.t) |
| 39 | +[Assignment](assignment.qmd){.btn .btn-secondary} |
51 | 40 |
|
52 | | -print(sol.y) |
53 | | -``` |
| 41 | +## Learning Goals |
54 | 42 |
|
55 | | -**How does `solve_ivp` work?** Let's understand its parameters (copied from the [documentation](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html)): |
| 43 | +- Formulate an ODE or a small ODE system in Python. |
| 44 | +- Use `scipy.integrate.solve_ivp()` to integrate an initial value problem. |
| 45 | +- Interpret trajectories in time and connect them to parameters. |
| 46 | +- Recognize threshold, saturation, and bistability effects in simple models. |
| 47 | +- Reuse the same solver workflow across very different applications. |
56 | 48 |
|
57 | | -- `fun`: Right-hand side of the system: the time derivative of the state $y$ at time $t$. The calling signature is `fun(t, y)`, where `t` is a scalar and `y` is an ndarray with `len(y) = len(y0)`. Additional arguments need to be passed if `args` is used (see documentation of `args` argument). `fun` must return an array of the same shape as `y`. |
58 | | - In our example, `exponential_decay` defines the ODE $\frac{dy}{dt} = -0.5y$, which models exponential decay $y(t) = y_0 e^{-0.5t}$. |
59 | | -- `t_span`: Interval of integration $(t_0, t_f)$. The solver starts with $t=t_0$ and integrates until it reaches $t=t_f$. Both $t_0$ and $t_f$ must be floats or values interpretable by the float conversion function. |
60 | | -- `y0`: Initial state. For problems in the complex domain, pass `y0` with a complex data type (even if the initial value is purely real). |
| 49 | +## Flow of This Session |
61 | 50 |
|
62 | | -You can also specify additional parameters: |
63 | | -- `method`: Integration method to use. Common choices include `'RK45'` (default), `'RK23'`, `'DOP853'`, `'Radau'`, `'BDF'`, and `'LSODA'`. |
64 | | -- `t_eval`: Times at which to store the computed solution, must be sorted and lie within `t_span`. If None (default), use points selected by the solver. |
65 | | -- `args`: Additional arguments to pass to the user-defined functions. If, for example, `fun` has the signature `fun(t, y, a, b, c)`, then `args=(a, b, c)`. |
| 51 | +- Write the right-hand side of the model. |
| 52 | +- Choose parameters, a time span, and an initial condition. |
| 53 | +- Integrate with `solve_ivp()`. |
| 54 | +- Plot the solution and interpret the result. |
| 55 | +- Compare how the behavior changes when you vary parameters or initial data. |
66 | 56 |
|
67 | | -## Systems of ODEs |
| 57 | +## What do we need? |
68 | 58 |
|
69 | | -For multiple coupled equations, return a list or array of derivatives: |
| 59 | +### `scipy.integrate.solve_ivp` |
| 60 | +
|
| 61 | +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`. |
70 | 62 |
|
71 | 63 | ```python |
72 | | -def sir_model(t, y, beta, gamma): |
73 | | - S, I, R = y |
74 | | - N = S + I + R |
75 | | - |
76 | | - dSdt = -beta * S * I / N |
77 | | - dIdt = beta * S * I / N - gamma * I |
78 | | - dRdt = gamma * I |
79 | | - |
80 | | - return [dSdt, dIdt, dRdt] |
| 64 | +import numpy as np |
| 65 | +from scipy.integrate import solve_ivp |
| 66 | +
|
| 67 | +
|
| 68 | +def exponential_decay(t, y): |
| 69 | + return -0.5 * y |
| 70 | +
|
| 71 | +
|
| 72 | +sol = solve_ivp( |
| 73 | + exponential_decay, |
| 74 | + t_span=(0, 10), |
| 75 | + y0=[2.0], |
| 76 | + t_eval=np.linspace(0, 10, 200), |
| 77 | +) |
81 | 78 | ``` |
82 | 79 |
|
83 | | -## Best Practices |
| 80 | +### `numpy` |
| 81 | +
|
| 82 | +Use NumPy arrays for initial conditions, evaluation grids, and post-processing of trajectories. |
84 | 83 |
|
85 | | -1. **Always check convergence**: Plot solutions at different tolerances |
86 | | -2. **Use appropriate methods**: `'RK45'` (default) works well for most problems |
87 | | -3. **Vectorize when possible**: Makes code faster and cleaner |
88 | | -4. **Document parameters**: Keep track of units and meanings |
89 | | -5. **Validate against known solutions**: Test your implementation |
| 84 | +### `matplotlib.pyplot` |
90 | 85 |
|
91 | | -## Next Steps |
| 86 | +Use Matplotlib to visualize time series, parameter comparisons, and phase-line style diagnostics. |
92 | 87 |
|
93 | | -Apply these techniques to the classical models in Session 1: |
94 | | -- SIR epidemiological model |
95 | | -- Spruce budworm population dynamics |
96 | | -- Michaelis–Menten enzyme kinetics |
| 88 | +### Why start with `solve_ivp()`? |
97 | 89 |
|
98 | | -## Resources |
| 90 | +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. |
99 | 91 |
|
100 | | -- [SciPy solve_ivp documentation](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html) |
101 | | -- [Python for ODEs tutorial](https://python-numerical-methods.berkeley.edu/notebooks/chapter22.00-ODE-Initial-Value-Problems.html) |
| 92 | +::: {.callout-tip} |
| 93 | +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. |
| 94 | +::: |
0 commit comments