Appendix¶
Editor: Weipeng Xu
Last modified: 24/06/2025
Reference: Sundnes, Joakim. Solving Ordinary Differential Equations in Python. Springer Nature, 2024.
Programming of difference equations¶
A sequence is simply a collection of numbers: $$ x_0, x_1, x_2,\cdots, x_n,\cdots $$ For certain sequence, we can write the complete sequence in a compact form, for example: $$ \begin{aligned} &(x_n)^{\infty}_{n=0},\quad x_n = 2n+1\\ &(x_n)^{\infty}_{n=0},\quad x_n = n!\\ &(x_n)^{\infty}_{n=0},\quad x_n = \sum_{j=0}^n\frac{x^j}{j!}\\ \end{aligned} $$ In real-life applications, sequences are usually finite, i.e., $(x_n)_{n=0}^N$. We usually difference equations to define the sequence, which defines $x_n$ by a relation involving $x_{n-1}$ and possibly earlier terms.
More examples of difference equations¶
Fibonacci numbers (2-nd order difference equation) $$ x_n = x_{n-1} + x_{n-2},\quad x_0=x_1=1 $$
Logistic growth (nonlinear difference equation)
$$ x_n = x_{n-1} + \frac{\rho}{100}x_{n-1}(1-\frac{x_{n-1}}{R}) $$
Factorial $$ x_n = nx_{n-1},\quad x_0=1 $$
Newton's method $$ x_n = x_{n-1} - \frac{f(x_{n-1})}{f^{\prime}(x_{n-1})} $$
def Newton(f, dfdx, x, epsilon=1e-7, max_n=100):
n = 0
while abs(f(x)) > epsilon and n <= max_n:
x = x - f(x) /dfdx(x)
n += 1
return x, n, f(x)
Systems of difference equations¶
The growth of two species, a predator and a prey, in the same ecosystem can be described by:
- Lotke-Volterra model $$ \begin{aligned} &x_n = x_{n-1} + ax_{n-1} - bx_{n-1}y_{n-1}\\ &y_n = y_{n-1} + dbx_{n-1}y_{n-1}-cy_{n-1} \end{aligned} $$ where $a$ is the natural growth rate of the prey in the absence of predators, $b$ is the death rate of prey per encounter of prey and predator, $c$ is the natural death rate of predators in the absence of food (prey), and $d$ is the efficiency of turning predated prey into predators.
Taylor Series and Approximations¶
Polynomial approximations have been used for centuries to compute exponentials, trigonometric functions and others. The most famous one is the Taylor series: $$ \begin{aligned} f(x) &= \sum_{k=0}^{\infty}\frac{1}{k!}(\frac{d^kf(0)}{dx^k})x^k\\ &=f(0) + f^{\prime}(0)x + \frac{1}{2}f^{\prime\prime}(0)x^2 + \frac{1}{6}f^{\prime\prime\prime}(0)x^3+\cdots \end{aligned} $$ For example: $$ \begin{aligned} &e^x = \sum_{k=0}^{\infty}\frac{x^k}{k!}\\ &\sin x = \sum_{k=0}^{\infty}(-1)^k\frac{x^{2k+1}}{(2k+1)!} \end{aligned} $$ The approximation accuracy will improve as $N$ is increased. We usually truncate the series to $N$, and we can also shift the variables: $$ f(x)\approx\sum_{k=0}^N\frac{1}{k!}(\frac{d^kf(a)}{dx^k})(x-a)^k $$ The most popular choice is $N=1$, which provides a reasonable approximation close to $x=0$.
We can reformulate the taylor series for $e^x$ around $0$ with $n$ terms $$ e_n = \sum_{k=0}^{n-1}\frac{x^k}{k!}=\sum_{k=0}^{n-2}\frac{x^k}{k!} + \frac{x^{n-1}}{(n-1)!} $$ as a difference equation $$ e_n = e_{n-1} + \frac{x^{n-1}}{(n-1)!},\quad e_0=0 $$ To avoid the repeated multiplications, we can design the following difference system: $$ \begin{aligned} &e_n = e_{n-1} + a_{n-1},\quad e_0=0\\ &a_n = \frac{x}{n}a_{n-1},\quad a_0=1 \end{aligned} $$
import numpy as onp
x = 0.5
N = 5
index_set = range(N + 1)
a = onp.zeros(len(index_set))
e = onp.zeros(len(index_set))
a[0] = 1
print(f'Exact: exp({x}) = {onp.exp(x):4.5f}')
for n in index_set[1:]:
e[n] = e[n-1] + a[n-1]
a[n] = x / n * a[n-1]
print(f'n = {n}, appr. = {e[n]:4.5f}, e = {onp.abs(e[n]-onp.exp(x)):4.5f}')
Exact: exp(0.5) = 1.64872 n = 1, appr. = 1.00000, e = 0.64872 n = 2, appr. = 1.50000, e = 0.14872 n = 3, appr. = 1.62500, e = 0.02372 n = 4, appr. = 1.64583, e = 0.00289 n = 5, appr. = 1.64844, e = 0.00028