用 python 解簡單的常微分方程式

# sudo apt install python3-pip
# python3 -m venv venv
# cd venv
# . bin/activate
# pip3 install numpy matplotlib
import numpy as np
import matplotlib.pyplot as plt

n = 100
t = np.linspace(0, 1.0, n)
y = np.zeros(n)
y[0] = 5     # initial boundry condition

# x(t) := a*exp(w*t) - b => b = a*exp(0) - x(0), dy/dt = w*a*exp(w*t) = w * [x(t) + b]
w = -5       # weight, - to decay, + to explode
a = 10       # amplitude
b = a - y[0] # bias
x = a * np.exp(w*t) - b # time evolution function, x(0) = a - b => b = a - x[0], 常微分 dx/dt = w*(x +b)
df_dt = lambda t, x: w * (x + b) # def df_dt(t, x): return w * (x + b)

z = y
dt = t[1] - t[0]
for i in range(n - 1): # to update [i + 1]
    # 1. 4th Order Runge-Kutta method to solve y for dy/dt = f(y, t)
    k1 = df_dt(t[i], y[i])
    k2 = df_dt(t[i] + dt / 2, y[i] + dt * k1 / 2)
    k3 = df_dt(t[i] + dt / 2, y[i] + dt * k2 / 2)
    k4 = df_dt(t[i] + dt    , y[i] + dt * k3    )
    y[i + 1] = y[i] + dt * (k1 + 2 * k2 + 2 * k3 + k4) / 6
    # 2. 1st Order Euler method to solve z for dz/dt = f(z, t)
    z[i + 1] = z[i] + dt * df_dt(t[i], z[i])
    
fig = plt.figure()

plt.subplot(311)
plt.plot(y)
plt.ylabel('y: RK4')

plt.subplot(312)
plt.plot(z)
plt.ylabel('z',  rotation=75)
plt.yticks(rotation=90)

plt.subplot(313)
plt.plot(x)
plt.ylabel('x:=a*exp(wt)-b')
plt.xlabel(f"t, a={a}, b={b}, w={w}")
plt.show()