《微分方程与动力系统（Differential Equations and Dynamical Systems）》听课笔记1

为了深入理解fMRI分析的原理，还是要学习一些微分方程与动力系统的东东，这是我的学习笔记，内容还没有很好的整理。

视频在这里1-Differential Equations and Dynamical Systems Overview_哔哩哔哩_bilibili

课程的原始地址在这里：ME 564 - Mechanical Engineering Analysis (washington.edu)

1. Overview

线性代数是微分方程的一部分，如果可以的话，大学课程应该将这两部分一起传授。这门课程除了基本的模型介绍以外，还会结合Matlab/Python的代码来进行具体的实现。

这门课程会介绍以下几个主题：

微积分：简单回顾微积分的计算，重点是泰勒展开
简单常微分方程（Simple Ordinary Differential Equations, ODEs）：即 $\dot{x}=\lambda x$ ，可以用于描述很多系统（如最常见的兔子吃草问题）
常微分方程组（System of ODEs）：即 $\dot X = AX$ ，常见的问题如在兔子吃草过程中引入其他动物（比如狼）
特征值与特征向量
非线性系统与混沌：即 $\dot X = f(X)$ ，常见的问题如行星运动
数值与计算：用Python或者Matlab

2. 微积分回顾：导数

导数反映的是一个函数 $f(x)$ 对于独立变量 $x$ 的变化率。其定义如下：

$\begin{equation} \frac{d f}{dx}=\lim\limits_{\Delta x \to 0}\frac{f(x+\Delta x)-f(x)}{\Delta x} \end{equation}$

链式法则（Chain Rule)：

$\frac{d}{d x}f(g(x))=\frac{df}{dx}(g(x))\frac{dg(x)}{dx}=f'(g(x))g'(x)$

例子： $f(x)=sin(x)$ 和 $g(x)=x^3$ ，则：

$f(g(x))=sin(x^3)$

$\frac{d}{dx}f(g(x))=3cos(x^3)x^2$

3. 向量和矩阵建模

【一个天气模型】在这个模型中，天气有三种状态，分别是多云(cloudy)、下雨(rainy)和晴天(nice),已知：
$\begin{array}{ll} P(cloudy|cloudy)=0.5 \\ P(rainy|cloudy)=0.25 \\ P(nice|cloudy)=0.25 \\ P(cloudy|rainy)=0.25 \\ P(rainy|rainy)=0.5 \\ P(nice|rainy)=0.25 \\ P(cloudy|nice)=0.5 \\ P(rainy|nice)=0.5 \\ P(nice|nice)=0 \end{array}$
上面公式表示一个概率，例如第一条表示如果前一天的天气为多云时，则第二天的天气为多云的概率是0.5，后面以此类推。此时，我们可以用一个向量来表示当前的天气情况：
$X_{today}=\begin{bmatrix} pr(R)\\ pr(N)\\ pr(C)\\ \end{bmatrix}=\begin{bmatrix} 0 \\ 1 \\ 0 \\ \end{bmatrix}$
上式中表示当前的天气为晴天。则第二天的前期状态可以根据之前的概率得到：
$X_{tomorrow}=AX_{today}=\begin{bmatrix} 0.5 \\ 0 \\ 0.5 \\ \end{bmatrix}$
不难看出， $A$ 为：
$A=\begin{bmatrix} 0.5&0.5&0.25 \\ 0.25&0&0.25 \\ 0.25&0.5&0.5\\ \end{bmatrix}$
此时，可以换一种形式来表示：
$\begin{align} X_{k+1}=&AX_k \\ =&A^{k}X_1 \end{align}$

clear all, close all, clc
A = [0.5 0.5 0.25;
0.25 0.0 0.25;
0.25 0.5 0.5];
xtoday = [1; 0; 0];
for k=1:50
    k
    xtomorrow = A * xtoday
    xtoday = xtomorrow;
end

运行上述代码可以看到，这个模型中天气会收敛到一个稳定的状态， $X=\begin{bmatrix}0.4\\0.2\\0.4\end{bmatrix}$ 。

有意思的事情是，这个稳态值和矩阵特征向量与特征值的关系。可以用eig函数来得到矩阵的特征值与特征向量：

[T,D]=eig(A)
T =

   -0.6667   -0.7071    0.4082
   -0.3333    0.0000   -0.8165
   -0.6667    0.7071    0.4082


D =

    1.0000         0         0
         0    0.2500         0
         0         0   -0.2500

此时如果对第一个特征值对应的特征向量进行归一化处理，可以得到：

T(:,1)/sum(T(:,1))

ans =

    0.4000
    0.2000
    0.4000

可以看到其最终结果就是系统的稳态分布，因此可以用这种方法来计算得到稳态分布，而不是使用迭代的方法。

上面过程的Python代码如下：

import numpy as np
import matplotlib.pyplot as plt
A = np.array([[0.5, 0.5, 0.25],[0.25, 0, 0.25],[0.25, 0.5, 0.5]])
xtoday = np.array([[1], [0], [0]])
the_weather = np.zeros((3,50))
for k in range(50):
    the_weather[:,k] = xtoday[:,0]
    xtomorrow = A@xtoday
    xtoday = xtomorrow
    print(k)
    print(xtomorrow)
    
plt.plot(the_weather.transpose())
plt.grid(True)

天气收敛

4. 简单常微分方程和它的指数解

简单常微分方程的形式： $\dot x = \lambda x$

【例子】假设在一个草原上有 $x$ 只兔子，兔子的增长速度与兔子的数量成正比。这个问题就可以描述为上面的公式：
$\begin{array}{ll} \dot x = \frac{dx}{dt} = \lambda x \\ x = x(t) \\ \end{array}$
其中， $\lambda$ 表示增长率， $x$ 是时间的函数。后面我们还会延续这个模型，引入一些其他的动物。

4.1 微分方程求解

要想计算这个微分方程，有很多种方法：

【method 1】
$\begin{align}{ll} \because \frac{dx}{dt} =& \lambda x \\ \therefore dx =& \lambda x dt \\ \Rightarrow \int dx =& \int \lambda x dt \\ \Rightarrow ln(x(t)) =& \lambda t + C \\ \Rightarrow x(t) =& e^{\lambda t +C}= Ke^{\lambda t}\\ \end{align}$
根据上式，不难得到以下结论：
$\begin{array}{ll} x(0)=K \\ \therefore x(t) = e^{\lambda t} x(0) \end{array}$

易知上述微分方程的差分形式为： $\Delta x=rx(t)\Delta t$

5. 使用幂级数求解微分方程

对于常微分方程 $\dot x=ax\Rightarrow x(t)=e^{at}x_0$ ，假设函数 $x(t)$ 可以写成幂级数的形式，可知：
$\begin{align} x(t) = c_0+c_1x+c_2x^2+...\\ \end{align}$
则
$\dot x(t)=c_1+2c_2x+3c_3x^2+...$
同时：
$\begin{align} ax(t) = ac_0+ac_1x+ac_2x^2+...\\ \end{align}$
$t$ 的对应系数应该相等，因此：
$\begin{align} c_1=&ac_0=ax_0 \\ 2c_2=&ac_1\Rightarrow c_2=\frac{a}{2}c_1=\frac{a^2}{2}x_0 \\ 3c_3 =&ac_2 \Rightarrow c_3=\frac{a}{3}c_2=\frac{a^3}{3!}x_0 \\ ...\\ (n+1)c_{n+1}=&ac_n \Rightarrow c_{n+1}=\frac{a}{(n+1)!}x_0 \end{align}$
所以：
$\begin{align} x(t)=&x_0+x_0at+x_0\frac{a^2t^2}{2!}+x_0\frac{a^3t^3}{3!}+...\\ =&x_0[1+at+\frac{a^2t^2}{2!}+\frac{a^3t^3}{3!}+...] \\ =&x_0e^{at} \end{align}$

6. 泰勒级数和幂级数

一个函数 $f(x)$ 可以在一个点 $x$ （ $f(x)$ 在 $x$ 处光滑）被泰勒展开：
$f(x+\Delta x)=f(x)+\frac{df(x)}{dx}\Delta x+\frac{d^2f(x)}{dx^2}\frac{\Delta x^2}{2!}+\frac{d^3f(x)}{dx^3}\frac{\Delta x^3}{3!}+...+\frac{d^nf(x)}{dx^n}\frac{\Delta x^n}{n!}+...$
如果在特定的点 $a$ （即 $x=a+\Delta x$ )展开，上述公式也可以写为如下形式：
$f(x)=f(a)+\frac{df}{dx}(x-a)+\frac{d^2f}{dx^2}\frac{(x-a)^2}{2!}+\frac{d^3f}{dx^3}\frac{(x-a)^3}{3!}+...$
【例1】当 $f(x)=sin(x)$ 时，其泰勒展开为：
$\begin{align} sin(x)=&sin(0)+cos(0)x-\frac{sin(0)}{2!}x^2-\frac{cos(0)}{3!}x^3+... \\ =&x-\frac{x^3}{3!}+\frac{x^5}{5!}-\frac{x^7}{7!}+\frac{x^9}{9!}+... \\ \end{align}$
【例2】当 $f(x)=cos(x)$ 时，其泰勒展开为：
$cos(x)=1-\frac{x^2}{2!}+\frac{x^4}{4!}-\frac{x^6}{6!}+\frac{x^8}{8!}+...$
Matlab代码如下：

clear all, close all
x = -10:0.01:10;
y = sin(x);
plot(x, y, 'k', 'LineWidth',2)
axis([-10 10 -10 10])
grid on, hold on

%% 一阶泰勒展开
P = [1 0];
yT1 = polyval(P, x);
plot(x, yT1, 'b--', 'LineWidth', 1.2)

%% 三阶泰勒展开
P = [-1/factorial(3) 0 1 0];
yT3 = polyval(P, x);
plot(x, yT3, 'r--', 'LineWidth', 1.2)

%% 五阶泰勒展开
P = [1/factorial(5) 0 -1/factorial(3) 0 1 0];
yT5 = polyval(P, x);
plot(x, yT5, 'g--', 'LineWidth', 1.2)

%% 七阶泰勒展开
P = [-1/factorial(7) 0 1/factorial(5) 0 -1/factorial(3) 0 1 0];
yT7 = polyval(P, x);
plot(x, yT7, 'm--', 'LineWidth', 1.2)

%% 九阶泰勒展开
P = [1/factorial(9) 0 -1/factorial(7) 0 1/factorial(5) 0 -1/factorial(3) 0 1 0];
yT9 = polyval(P, x);
plot(x, yT9, 'c--', 'LineWidth', 1.2)

talyor.png

Python的泰勒展开结果如下：

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
from scipy.special import factorial

x = np.linspace(-10, 10, 2000)
y = np.sin(x)
plt.plot(x, y, 'k', linewidth=2)
plt.xlim(-10, 10)
plt.ylim(-10, 10)
plt.grid(True)

# 一阶泰勒展开
P = [1, 0]
yT1 = np.polyval(P, x)
plt.plot(x, yT1, 'b--', linewidth=1.2)

# 三阶泰勒展开
P = [-1/factorial(3), 0, 1, 0]
yT3 = np.polyval(P, x)
plt.plot(x, yT3, 'r--', linewidth=1.2)

# 五阶泰勒展开
P = [1/factorial(5), 0, -1/factorial(3), 0, 1, 0]
yT5 = np.polyval(P, x)
plt.plot(x, yT5, 'g--', linewidth=1.2)

# 七阶泰勒展开
P = [-1/factorial(7), 0, 1/factorial(5), 0, -1/factorial(3), 0, 1, 0]
yT7 = np.polyval(P, x)
plt.plot(x, yT7, 'm--', linewidth=1.2)

# 九阶泰勒展开
P = [1/factorial(9), 0, -1/factorial(7), 0, 1/factorial(5), 0, -1/factorial(3), 0, 1, 0]
yT9 = np.polyval(P, x)
plt.plot(x, yT9, 'c--', linewidth=1.2)

image.png

7.指数函数的泰勒级数和欧拉公式

本节课利用泰勒级数推导出欧拉公式，过程如下：
$\begin{align} e^x =& 1 + x+\frac{x^2}{2!}+\frac{x^3}{3!}+\frac{x^4}{4!}+\frac{x^5}{5!}+... \\ =& \sum_{k=0}^{\infty}\frac{x^k}{k!}\\ \Rightarrow e^{ix}=& 1 + ix+\frac{(ix)^2}{2!}+\frac{(ix)^3}{3!}+\frac{(ix)^4}{4!}+\frac{(ix)^5}{5!}+...\\ =&(1-\frac{x^2}{2!}+\frac{x^4}{4!}-\frac{x^6}{6!}+...) + i(x-\frac{x^3}{3!}+\frac{x^5}{5!}+...) \\ =&cos(x)+isin(x) \end{align}$
得到上述公式后，后面的计算会变得更加容易

8. 用二阶常微分方程描述谐振子

此处，我们引入了一个物理学中的常见概念——谐振子，并使用微分方程来对其状态进行建模。

问题描述如下：在一个光滑平面上，有一个质量为 $m$ 的物体，其一端用弹簧与墙壁相连，弹簧的系数为k，物体偏移静止状态（即除支撑力外，不受任何弹簧弹力）的位移为 $x$ ，则可知：

$\begin{align} F=&-Kx=ma \\ v=&\frac{dx}{dt} \\ a=& \frac{dv}{dt}= \frac{d^2x}{dt^2} \\ \end{align}$
令 $m=1$ ， $K=1$ ，则：
$\begin{align} \frac{d^2x}{dt^2}=-x \Rightarrow \ddot x = -x \end{align}$
求解上述方程的方法有很多种：

【Method 1】瞪眼法

因为 $x(0)=0$ ， $\dot x(0)=v_0=0$ ，结合前面的泰勒展开结果，充分考虑这些初始条件，则可以猜测
$\begin{align} x(t)=cos(t)x_0 \end{align}$
则
$\begin{align} \dot x(t)=&-sin(t)x_0\\ \ddot x(t)=&-cos(t)x_0\\ \end{align}$
而上述结果满足之前的 $\ddot x = -x$ ，于是瞪眼成功。

对于一般的 $K$ 和 $m$ ，
$\begin{align} x(t)=cos(\sqrt\frac{K}{m}t)x(0) \end{align}$

【Method 2】泰勒级数展开

$\begin{align} x(t)=&c_0+c_1t+c_2t^2+c_3t^3+c_4t^4+...\\ \dot x(t)=&c_1+2c_2t+3c_3t^2+4c_4t^3+5c_5t^4... \\ \ddot x(t)=&2c_2+3\cdot2c_3t+4\cdot3c_4t^2+5\cdot4c_5t^3... \\ \end{align}$

因此
$\begin{align} c_0=&x_0 \\ c1=&v_0 \\ 2c_2=&-c_0 \Rightarrow c_2=-\frac{1}{2}x_0 \\ 3\cdot2c_3=&-c_1 \Rightarrow c_3=-\frac{1}{3!}v_0 \\ 4\cdot3c_4=&-c_2 \Rightarrow c_4 = \frac{1}{4!}x_0 \\ 5\cdot4c_5=&-c_3 \Rightarrow c_5 = \frac{1}{5!}v_0 \\ ... \end{align}$
将其带回，可得：
$\begin{align} x(t)=&x_0+v_0t-\frac{1}{2!}x_0t^2-\frac{1}{3!}v_0t^3+\frac{1}{4!}x_0t^4+\frac{1}{5!}v_0t^5+... \\ =&x_0(1-\frac{t^2}{2!}+\frac{t^4}{4!}-\frac{t^6}{6!}+...)+v_0(t-\frac{t^3}{3!}+\frac{t^5}{5!}+...) \\ =&cos(t)x_0+sin(t)v_0 \end{align}$

【Method 3】瞪眼法2

一个函数求导之后仍然是自身的函数，最容易想到的就是指数函数，因此不妨猜测 $x$ 是一种指数函数。
$\begin{align} x(t)=&e^{\lambda t}\\ \dot x=&\lambda e^{\lambda t} \\ \ddot x=&\lambda^2 e^{\lambda t} \end{align}$

将其代入微分方程，可得：
$\begin{align} \lambda^2 e^{\lambda t}=&-e^{\lambda t} \Rightarrow \lambda ^2=-1 \end{align}$
所以：
$\begin{align} x(t)=&c_1e^{it}+c_2e^{-it}\\ =&(c_1+c_2)cos(t)+i(c_1-c_2)sin(t) \end{align}$

【Method 4】Suspend variable

在该方法中，引入一个新的变量 $v=\frac{dx}{dt}$
$\left. \begin{matrix} \dot x=v \\ \dot v=-x \end{matrix} \right\}\Rightarrow \frac{d}{dt} \begin{bmatrix} x \\ v \end{bmatrix}= \begin{bmatrix} 0 & 1 \\ -1 & 0 \end{bmatrix} \begin{bmatrix} x \\ v \end{bmatrix}$
令：
$A= \begin{bmatrix} 0 & 1 \\ -1 & 0 \end{bmatrix} \\ \underline x= \begin{bmatrix} x \\ v \end{bmatrix}$
则：
$\dot {\underline x} = A \underline x \Rightarrow \underline x=e^{At}\underline x(0)$

9. 二阶ODE，阻尼震荡器

阻尼震荡器

新的问题变成了添加了一个阻尼的震荡子，如上图所示，此时模型变为：
$\begin{align} F=ma \Rightarrow& -Kx-d\dot x = m\ddot x \\ \Rightarrow& m \ddot x+d \dot x+ Kx = 0 \end{align}$
令
$\begin{align} \zeta =& \frac{d}{m} \\ \omega^2 =& \frac{K}{m} \end{align}$

所以：

$\ddot x+\zeta \dot x+\omega^2x=0$
看到这个模型，不难继续要使用瞪眼法，令：

$\left. \begin{align} x(t)=&e^{\lambda t} \\ \dot x(t) =& \lambda e^{\lambda t}\\ \ddot x(t) =& \lambda^2 e^{\lambda t} \end{align} \right\} \Rightarrow \lambda^2e^{\lambda t}+\zeta\lambda e^{\lambda t}+\omega^2 e^{\lambda t}=0$

所以：

$\lambda^2+\zeta\lambda+\omega^2=0$
这个方程被称为特征方程。

其解为：
$\lambda=\frac{-\zeta \pm \sqrt{\zeta^2-4\omega^2}}{2}$
最终
$x(t)=c_1e^{\lambda_1t}+c_2e^{\lambda_2t}$
还可以将模型写为矩阵形式：
$\left\{ \begin{align} \dot x =& v \\ \dot v =& -\omega^2x-\zeta v \\ \end{align} \right. \Rightarrow \frac{d}{dt} \begin{bmatrix} x \\ v \\ \end{bmatrix} = \begin{bmatrix} 0 & 1\\ -\omega^2 & -\zeta \\ \end{bmatrix} \begin{bmatrix} x \\ v\\ \end{bmatrix}$
下面我们在Matlab中对这个过程进行模拟，需要注意的是：
$\begin{align} \because \frac{dx}{dt} =& Ax \\ \therefore \Delta x =& Ax\Delta t \\ x_{k+1} =& x_{k}+Ax\Delta t \\ t_k =& k\Delta t \end{align}$
看懂了上式就可以对上述过程进行编码了，具体代码如下

clear all; close all;

w = 2*pi; 
d = 0.25;

% dt = Ax，此处和文中给出的公式并不一致，相当于k=.5，但是只要知道这是在做一些系统初始化的设置即可
A = [0 1; -w^2 -2*d*w];

dt = .01;   % time step
T = 10;     % amount of time to integrate

% 初始位置在2处，初始速度为0
x0 = [2; 0];

xF(:,1) = x0;
tF(1) = 0;
for k=1:T/dt
    tF(k+1) = k*dt;
    xF(:,k+1) = (eye(2) + dt*A)*xF(:,k);
end

plot(tF, xF(1,:), 'k')

% 使用龙格库塔（Runge Kutta）得到的结果
[t, xGood] = ode45( @(t, x) A*x, 0:dt:T, x0);
hold on
plot(t, xGood(:,1), 'r')
xlabel('Time (s)')
ylabel('Position (m)')
legend('Forward Euler', 'ODE45(RK4)')

figure
plot(xGood(:,1),xGood(:,2),'k')
xlabel('Position [m]')
ylabel('Velocity [m/s]')

位置随时间变化图

状态图

import numpy as np
from matplotlib import pyplot as plt
from scipy.integrate import solve_ivp

# play around with different 'd' and 'dt' to see various behavior!

w = 2 * np.pi  # natural frequency
d = .25        # damping ratio

# spring-mass-damper system
A = np.array([[0, 1], [-w ** 2, -2 * d * w]])  # \dot{x} = Ax

dt = 0.01  # time step
T = 10     # amount of time to integrate
n = int(T / dt)
t = np.linspace(0, T, n)

x0 = [2, 0]  # initial condition (x=2, v=0)

# iterate forward Euler
xF = np.zeros((2, n))
xF[:, 0] = x0
for k in range(n - 1): 
    xF[:, k + 1] = (np.eye(2) + dt * A) @ xF[:, k]

# compute better integral using built-in python code
# 4th-order Runge Kutta 
def linear_ode(t, x):
    return A @ x  # @ symbol for matrix-vector product here

linear_ode_solution = solve_ivp(linear_ode, (0, T), x0, t_eval=t) 
xGood = linear_ode_solution.y

plt.figure(figsize=(20, 4))
plt.subplot(1, 2, 1)
plt.plot(t, xF[0, :], 'k')
plt.plot(t, xGood[0, :],'r')
plt.xlabel('Time [s]')
plt.ylabel('Position [m]')
plt.legend(['Forward Euler', 'ODE45 (RK4)'])
plt.grid(True)

plt.subplot(1, 2, 2)
plt.plot(xF[0, :], xF[1, :], 'k')
plt.plot(xGood[0, :], xGood[1, :], 'r')
plt.xlabel('Position [m]')
plt.ylabel('Velocity [m/s]')
plt.legend(['Forward Euler', 'ODE45 (RK4)'])
plt.grid(True)