Scientific Computing

# Scientific Computing

</div>
</div>

# Part.3
# Solving Lineaer Systems of Equations
# 线性方程组的求解

</div>
</div>

## 线性方程组

**定义：** 设 $A$ 是 $n\times n$ 矩阵，$\mathbf{b}$ 是 $n$ 维列向量，线性方程组是指形如
$$A\mathbf{x}=\mathbf{b}
$$的方程组，其中 $\mathbf{x}$ 是未知的 $n$ 维列向量。

**例如：** 设给定一个线性方程组
$$
\begin{cases}
2x_1 + x_2 - x_3 = 8 \\\\
-3x_1 - x_2 + 2x_3 = -11 \\\\
-2x_1 + x_2 + 2x_3 = -3
\end{cases}
$$

则可表示为矩阵形式 $A\mathbf{x}=\mathbf{b}$，其中

$$A=\begin{pmatrix}2 & 1 & -1 \\\\
-3 & -1 & 2 \\\\
-2 & 1 & 2 \end{pmatrix},\quad \mathbf{b}=\begin{pmatrix}8 \\\\ -11 \\\\ -3 \end{pmatrix}.$$

## 线性方程组与其可解性

何时有解？解是否唯一？
- **唯一解**：当 $A$ 是非奇异矩阵（即可逆矩阵）时，线性方程组有唯一解，解为 $\mathbf{x}=A^{-1}\mathbf{b}$ ：
  - $A$ 是可逆的；
  - $\det(A) \neq 0$；
  - $A$ 的秩 $\text{rank}(A) = n$；
  - $A$ 的列（行）向量组线性无关；
  - $A$ 的特征值均非零；
  - $A$ 的奇异值 (singular values) 均不为零；

- **无解或无穷多解**：当 $A$ 是奇异矩阵时，线性方程组可能无解或有无穷多解。
  - **无解**：当 $A$ 的秩小于增广矩阵 $[A|\mathbf{b}]$ 的秩时，线性方程组无解；

## 求解：形如三角形的方程组

在攻克一般方程组之前，我们先解决一种特殊形式：三角形系统（Triangular Systems）。这是高斯消元法的最后一步，也是算法的基础。

**定义：**上三角矩阵是指主对角线下方全为 0 的矩阵。例如：
$$
A=\begin{pmatrix}2 & -1 & 0 \\\\
0 & 3 & 4 \\\\
0 & 0 & 5 \end{pmatrix}.
$$

求解 $Ax=b$ （其中 $A$ 为上三角）的方法是回代法（Back Substitution）。
- 从最后一个方程开始，逐步向上求解每个变量；

## 回代法计算量分析：

设有 $n$ 个未知数的上三角方程组：
$$\begin{cases}
a_{11}x_1 &+ a_{12}x_2 &+ \cdots &+ a_{1n}x_n &= b_1 \\\\
&a_{22}x_2 &+ \cdots &+ a_{2n}x_n &= b_2 \\\\
& & \vdots & & \\\\
 & & & a_{nn}x_n &= b_n
\end{cases}$$

回代法的计算量为？ $\color{red}{O(n^2)}$

</div>
<div class="col">

分析一下：

> **第$1$步：**
> *解出 $x_n$ : 一次除法；*
> 
> **第$k$步：**
> *$1$ 次除法，$k-1$ 次乘法和 $k-1$ 次加法，解出 $x_{n-k+1}$；*
> 
> **总计：**$n$ 次除法，$\displaystyle\frac{(n-1)\cdot n}{2}$ 次乘法和加法。

</div>
</div>

## 下三角方程组的求解

<p style="color: blue;">
其实与上三角方程组类似，只不过是从第一个方程开始，逐步向下求解每个变量。在此不再赘述。
</p>

## Guauss 消元法简介

高斯消元法（Gaussian Elimination）是求解线性方程组的基本方法。其核心思想是通过初等行变换将系数矩阵化为上三角矩阵，从而利用回代法求解。

**步骤：**

> 1. 对于每一列，选择一个主元（pivot），通常选择该列中绝对值最大的元素，以减少数值误差；
> 2. 使用主元对下面的行进行消元操作，使得该列下方的元素变为 0；
> 3. 重复上述步骤，直到矩阵变为上三角形式；
> 4. 最后，使用回代法求解上三角方程组。

这是本章的核心算法。对于一般的稠密矩阵 (dense matrix)，高斯消元法通常是最佳选择 。

## 高斯消元法的计算量分析

### 高斯消元的算法步骤：

```python
for k in range(n-1):
 for i in range(k+1, n):
 factor = A[i][k] / A[k][k]
  for j in range(k, n):
   A[i][j] -= factor * A[k][j]
   b[i] -= factor * b[k]
```

</div>
<div class="col">

![Sketch of Gaussian Elimination](./image.png)
</div>
</div>

## 高斯消元法的计算量分析

> **计算量分析：**针对第 $k$ 步消元：
> - $n-k$ 次除法；
> - $(n-k)*(n-k)$ 次乘法和减法；
> - 总计约 $\displaystyle\frac{2}{3}n^3$ 次浮点运算

因此，计算大概有 $\displaystyle\sum_{k=1}^{n-1} \left[(n-k) + 2(n-k)^2\right] = 2\sum_{j=1}^{n-1}j(2j+1) $ 步。

总计的计算量为 $\color{red}{O(n^3)}$。

<p style="color: blue;">
如果矩阵的规模扩大 2 倍（例如从 $100 \times 100$ 到 $200 \times 200$），计算时间将增加约 $2^3 = 8$ 倍。这是大规模计算中必须考虑的瓶颈。😨
</p>

## MATLAB 中的运用 Backslash 运算符

在 MATLAB 中，求解线性方程组 $A\mathbf{x}=\mathbf{b}$ 最简单的方法是使用反斜杠运算符 `\`：

```python
A = [2 -1 0; 0 3 4; 0 0 5];
b = [8; -11; -3];
x = A\b;
```

该运算符会根据矩阵 $A$ 的性质选择最合适的算法来求解方程组，通常在矩阵很密的情况下会使用高斯消元法。

很稀疏的矩阵则会使用专门的稀疏矩阵算法，以提高计算效率。

## TODO：稀疏的线性方程组求解方法

> 对于稀疏矩阵（大部分元素为零的矩阵），直接使用高斯消元法可能会导致“填充”（fill-in）现象，即在消元过程中产生大量非零元素，增加计算量和存储需求。
> 
> 因此，针对稀疏矩阵的线性方程组求解，通常采用如共轭梯度法（Conjugate Gradient Method）、GMRES（Generalized Minimal Residual Method）等，这些方法通过迭代逼近解，适用于很大的稀疏矩阵。
>

## LU 分解简介

LU 分解是高斯消元法的自然延伸，也是实际工程中最常用的形式。

LU 分解（LU Decomposition）是将矩阵 $A$ 分解为一个下三角矩阵 $L$ 和一个上三角矩阵 $U$ 的乘积，即 $A=LU$。这种分解对于求解线性方程组非常有用。

> 什么是 LU 分解？
> - $U$ （Upper）是一个上三角矩阵，也就是高斯消元后得到的矩阵形式。
> - $L$ (Lower) 是一个下三角矩阵，主对角线上的元素均为 1，其余元素为消元过程中使用的乘数；
> 
> $A=LU$ ：把A看作是对应线性变换的矩阵，那么 $L$ 代表了消元过程中的行变换，而 $U$ 则是最终的上三角形式。

## 为什么可以求出解

给定线性方程组 $A\mathbf{x}=\mathbf{b}$，如果我们有 $A=LU$，则可以将方程组重写为：
$$LU\mathbf{x}=\mathbf{b}.$$

我们可以引入一个中间变量 $\mathbf{y}$，使得：
$$L\mathbf{y}=\mathbf{b},$$
$$U\mathbf{x}=\mathbf{y}.$$
这样，我们可以通过两个步骤来求解 $\mathbf{x}$：
1. **前代**：先解 $L\mathbf{y}=\mathbf{b}$，由于 $L$ 是下三角矩阵，可以使用前代法高效地求解 $\mathbf{y}$；
2. **回代**：然后解 $U\mathbf{x}=\mathbf{y}$，由于 $U$ 是上三角矩阵，可以使用回代法求解 $\mathbf{x}$。

## 为什么要使用 LU 分解求解线性方程组？

为什么还要用 $LU$ 分解，而不是直接使用高斯消元法？记录 $L$ 是不是多此一举？

> LU 分解的主要优势在于其高效性，特别是在需要多次求解同一系数矩阵但不同右端项的线性方程组时。也就是在于处理多组 $\mathbf{b}$ 的情况。（Multiple RHS）
>

## LU 分解求解多组线性方程组的计算量分析

运用 LU 分解求解线性方程组 $A\mathbf{x}=\mathbf{b}$ 的过程分为两个步骤：**预处理**：对矩阵 $A$ 进行 LU 分解，得到 $L$ 和 $U$；**多次求解**：对于每一个新的右端项 $\mathbf{b}$，先解 $L\mathbf{y}=\mathbf{b}$（前代），再解 $U\mathbf{x}=\mathbf{y}$（回代）。

> 具体来说：
> 1. **预处理**：通过一次 $LU$ 分解，我们可以将矩阵 $A$ 分解为 $L$ 和 $U$，这一步的计算量为 $O(n^3)$。
> 2. **多次求解**：对于每一个新的右端项 $\mathbf{b}$，我们只需进行两次三角形方程组的求解（一次前代，一次回代），每次的计算量为 $O(n^2)$。
> 
> **总计算量**：因此，对于 $m$ 个不同的右端项 $\mathbf{b}$，总的计算量为 $O(n^3 + mn^2)$，相比于每次都使用高斯消元法的 $O(mn^3)$，节省了大量的计算时间。

# Exercises for Part.3
# Solving Lineaer Systems of Equations

</div>
</div>

## 练习题解答

参见 [此处](/BLOGS/OTS/SC_ANS_P3.html)

# Part.7
# Eigenvalues and Eigenvectors 
# 特征值与特征向量

</div>
</div>

## 基本定义

对于一个 $n \times n$ 的矩阵 $A$，若存在标量 $\lambda$ 与非零向量 $x$ 满足以下线性关系：

$$Ax = \lambda x, \quad x \neq 0$$

则称 $\lambda$ 为矩阵 $A$ 的**特征值**，$x$ 为对应于 $\lambda$ 的**特征向量**。

## 一些性质

> - **非零：**必须强调 $x$ 不能为零向量，否则任何标量均满足方程，定义将失去意义。
> - **特征向量非唯一：**特征向量在方向上是确定的，但在模长上具有任意性。若 $x$ 为特征向量，则 $\alpha x$（$\alpha \neq 0$）亦为特征向量。
> - 在数值计算中，有时对其进行**归一化（Normalize）**处理，使其范数为1。
> - **复数：**即便矩阵 $A$ 为实矩阵，其特征值与特征向量仍可能为复数。比如旋转矩阵 $$R = \begin{pmatrix}0 & -1 \\\\ 1 & 0\end{pmatrix}$$ 的特征值为 $\lambda = i$ 和 $\lambda = -i$，对应的特征向量为复数向量。
> <p> </p>

## 特征多项式

我们将特征值问题 $Ax = \lambda x$ 重写为 $(A - \lambda I)x = 0$，其中 $I$ 为单位矩阵。为了使该方程有非零解（使得特征向量不为零向量），矩阵 $(A - \lambda I)$ 必须是奇异的，即其行列式为零：
$$\det(A - \lambda I) = 0.$$
该行列式展开后得到一个关于 $\lambda$ 的多项式，称为矩阵 $A$ 的**特征多项式**。特征值即为该多项式的根。

> - **特征多项式的性质**：特征多项式是一个 $n$ 次多项式，具有 $n$ 个根（包括重根），这些根即为矩阵的特征值。
> - **计算复杂度**：直接计算特征多项式的根在数值上不稳定，尤其对于高维矩阵。因此，实际应用中通常采用数值方法来近似计算特征值和特征向量。
>

## 特征值与多项式求根

一个常见的特征值应用是寻找多项式的根。如果我们有一个多项式

$$p(x) = c_0 + c_1 x + c_2 x^2 + \cdots + c_{n-1} x^{n-1} + c_{n} x^n,$$

我们可以构造一个伴随矩阵（companion matrix） $C$，其中：
$$C = \begin{pmatrix}
0 & 0 & 0 & \cdots & -\frac{c_0}{c_n} \\\\
1 & 0 & 0 & \cdots & -\frac{c_1}{c_n} \\\\
0 & 1 & 0 & \cdots & -\frac{c_2}{c_n} \\\\
\vdots & \vdots & \vdots & \ddots & \vdots \\\\
0 & 0 & 0 & \cdots & -\frac{c_{n-1}}{c_n}
\end{pmatrix}.$$

矩阵 $C$ 的特征值正是多项式 $p(x)$ 的根。

在 MATLAB 和 Octave 中，可以使用 `roots` 函数计算多项式的根，其就是通过构造这样一个伴随矩阵 $C$ 并计算其特征值来获得给定多项式的根。

## 通过数值方法求解特征值问题

对于大规模矩阵，直接计算特征多项式并求根是不切实际的。相反，我们使用数值方法来近似计算特征值和特征向量。

> - **幂迭代法（Power Iteration）：**一种简单的迭代方法，用于计算矩阵的主特征值和对应的特征向量。通过反复乘以矩阵并归一化，可以逐渐逼近主特征向量。
>
> <span style="color:grey;">
>
> -  **QR 算法：**一种更复杂但更强大的方法，适用于计算所有特征值和特征向量。该方法通过 QR 分解不断迭代矩阵，逐渐将其对角化，从而获得特征值。
> - **Lanczos 方法：**一种专门用于大型稀疏矩阵的迭代方法，适合计算少数几个特征值和特征向量。通过正交相似变换将对称矩阵转化为对称三对角矩阵，从而简化计算。
> - ......
> 
> </span>
>
> <p> </p>

## 幂迭代法的原理

假设 $A$ 是一个 $n \times n$ 矩阵，

对应的特征向量为 $x_1, x_2, \ldots, x_n$。

对于任意非零向量 $x_0$，我们可以将其表示为特征向量的线性组合：
$$x_0 = c_1 x_1 + c_2 x_2 + \cdots + c_n x_n,$$
其中 $c_1, c_2, \ldots, c_n$ 是标量系数，且 $c_1 \neq 0$。
通过反复乘以矩阵 $A$，我们得到：
$$A^k x_0 = c_1 \lambda_1^k x_1 + c_2 \lambda_2^k x_2 + \cdots + c_n \lambda_n^k x_n.$$
随着迭代次数 $k$ 的增加，项 $c_1 \lambda_1^k x_1$ 将主导整个表达式，因为 $\left|\lambda_1\right|$ 是最大的特征值的绝对值。因此，归一化后，$A^k x_0$ 将逐渐逼近主特征向量 $x_1$。

## 使用幂迭代法计算主特征值和特征向量

幂迭代法是一种简单而有效的数值方法，用于计算矩阵的主特征值和对应的特征向量。其基本思想是通过反复乘以矩阵并归一化，逐渐逼近主特征向量。

**算法步骤：**
 
> 1. 选择一个初始向量 $x_0$，通常为随机向量；
> 2. 对于每次迭代 $k$：
>    - 计算 $y_k = Ax_{k-1}$；
>    - 归一化 $x_k = \frac{y_k}{\|\|y_k\|\|}$；（为了防止数值溢出或下溢）
> 3. 重复步骤 2，直到收敛，即 $\|x_k - x_{k-1}\|$ 足够小。
> 4. 最终 $x_k$ 将近似为主特征向量，$\lambda_k = \frac{x_{k}^{T}Ax_k}{x_k^T x_k}$ 将近似为主特征值。
>
>

在这个算法中，选择哪种范数并不重要，但你应该保持一致，因为不同的范数会以不同的方式缩放向量。

## 幂迭代法的收敛性分析

幂迭代法的收敛速度取决于主特征值 $\lambda_1$ 与次大特征值 $\lambda_2$ 之间的比率 $\frac{|\lambda_2|}{|\lambda_1|}$。

具体来说：

> - **收敛比较快：**当 $\frac{|\lambda_2|}{|\lambda_1|}$ 较小时，幂迭代法收敛较快，因为次大特征值的影响迅速减小。
> - **收敛比较慢：**当 $\frac{|\lambda_2|}{|\lambda_1|}$ 接近 1 时，收敛速度较慢，因为次大特征值的影响较大，导致迭代过程中的误差减小较慢。
> - **可能不收敛：**如果主特征值不是唯一的最大值（即存在多个特征值具有相同的绝对值），则幂迭代法可能无法收敛到一个确定的特征向量。
> - **符号震荡：**如果主特征值是负数，迭代向量 $x_k$ 的符号可能在每一步跳变，导致收敛过程中的震荡现象。但这不代表不收敛，应关注特征值的估算值是否稳定。
>

## 逆转的幂迭代法

幂迭代法主要用于计算主特征值和特征向量，但有时我们需要计算最小特征值。为此，我们可以使用逆转的幂迭代法（Inverse Power Iteration）。
逆转的幂迭代法的基本思想也是通过求解线性方程组来间接计算最小特征值。

给定一对特征值 $\lambda$ 和特征向量 $x$，我们可以将特征值问题 $Ax = \lambda x$ 。

我们左乘 $A^{-1}$，得到
$$A^{-1}Ax = A^{-1}\lambda x \implies x = {\lambda}A^{-1}x,$$

也就是$$A^{-1}x = \frac{1}{\lambda}x,$$

这表明 $A^{-1}$ 的特征值 $\frac{1}{\lambda}$ 与 $A$ 的特征值 $\lambda$ 互为倒数。所以求解最小特征值及其特征向量的问题可以转化为求解 $A^{-1}x = \mu x$，其中 $\mu = \frac{1}{\lambda}$ 是 $A^{-1}$ 的最大特征值。

## 逆转的幂迭代法的具体步骤

只需要经过一些修改，我们就可以知道求解步骤如下：

> 1. 选择一个初始向量 $x_0$，通常为随机向量；
> 2. 对于每次迭代 $k$：
>    - 计算 $y_k = A^{-1}x_{k-1}$；
>    - 归一化 $x_k = \frac{y_k}{\|\|y_k\|\|}$；（为了防止数值溢出或下溢）
> 3. 重复步骤 2，直到收敛，即 $\|x_k - x_{k-1}\|$ 足够小。
> 4. 最终 $x_k$ 将近似为最小特征向量，$\mu_k = \frac{x_{k}^{T}Ax_k}{x_k^T x_k}$ 将近似为最小特征值。
>

## 继续寻找其他特征值和特征向量

一旦我们找到了最大/最小的特征值 $\lambda$ 和对应的特征向量 $x$，我们可以通过幂迭代法或逆转的幂迭代法继续寻找其他特征值和特征向量。

**以下是一个大致的步骤：**

> 我们可以使得特征值进行偏移：
> - 对于寻找次小特征值，可以考虑矩阵 $A - \sigma I$，其中 $\sigma$ 是一个接近最小特征值的值。通过逆转的幂迭代法应用于 $A - \sigma I$，可以找到次小特征值。
> 
> 如果之前矩阵的特征有 $\lambda_1, \lambda_2, \ldots, \lambda_n$，那么矩阵 $A - \sigma I$ 的特征值将变为 $\lambda_1 - \sigma, \lambda_2 - \sigma, \ldots, \lambda_n - \sigma$。通过选择适当的 $\sigma$，我们可以使得次小特征值成为新的最小特征值，从而应用逆转的幂迭代法找到它。
>

## 具体应用

特征值在生物数学中用于分析具有年龄结构的种群稳定性。

模型构建假设一种群分为不同年龄段 $P_1, P_2, \dots$，各阶段具有特定的出生率 $b_i$ 和死亡率 $d_i$。种群随时间的演化可表示为矩阵乘法：$$P^{n+1} = M P^n$$其中矩阵 $M$（通常称为 Leslie 矩阵）包含了出生率与存活率信息。

</div>
<div class="col">

其中矩阵 $M$ 为

$$
M = \begin{pmatrix}
b_1 & b_2 & b_3 & \cdots & b_k \\\\
s_1 & 0 & 0 & \cdots & 0 \\\\
0 & s_2 & 0 & \cdots & 0 \\\\
\vdots & \vdots & \vdots & \ddots & \vdots \\\\
0 & 0 & 0 & \cdots & s_{k-1}
\end{pmatrix}
$$

其中 $s_i = 1 - d_i$ 表示第 $i$ 个年龄段的存活率。

</div>

## 具体应用（续）

长期行为分析种群逐年的演化过程 $P^{n+1} = M P^n$ 在数学形式上等同于幂法的迭代过程。因此，种群的长期命运完全由转移矩阵 $M$ 的最大特征值 $\lambda_{max}$ 决定：
- $\lambda_{max} > 1$：种群数量将趋于无穷（指数增长）。
- $\lambda_{max} < 1$：种群将最终灭绝。
- $\lambda_{max} = 1$：种群数量将保持稳定。

$$
\begin{pmatrix}
P_1^{n+1} \\\\
P_2^{n+1} \\\\
P_3^{n+1} \\\\
\vdots \\\\
P_k^{n+1}
\end{pmatrix}=\begin{pmatrix}
b_1 & b_2 & b_3 & \cdots & b_k \\\\
s_1 & 0 & 0 & \cdots & 0 \\\\
0 & s_2 & 0 & \cdots & 0 \\\\
\vdots & \vdots & \vdots & \ddots & \vdots \\\\
0 & 0 & 0 & \cdots & s_{k-1}
\end{pmatrix}\begin{pmatrix}
P_1^{n} \\\\
P_2^{n} \\\\
P_3^{n} \\\\
\vdots \\\\
P_k^{n}
\end{pmatrix}
$$

# Exercises for Part.7
# Eigenvalues and Eigenvectors

</div>
</div>

## 练习题解答

参见 [此处](/BLOGS/OTS/homework7.html)

# Part.10
# Initial Value ODEs
# 常微分方程初值问题

</div>
</div>

## 常微分方程初值问题简介

对于高阶的常微分方程初值问题，我们通常将其转化为一阶常微分方程组来处理。设有一个 $n$ 阶常微分方程：
$$y^{(n)} = f\left(t, y, y', y'', \ldots, y^{(n-1)}\right),$$
我们可以引入新的变量：
$$y_1 = y, \quad y_2 = y',  \ldots, \quad y_n = y^{(n-1)},$$

</div>
<div class="col">

从而将原始的高阶常微分方程转化为以下一阶常微分方程组：

$$\begin{cases}
y_1' = y_2 \\\\
y_2' = y_3 \\\\
\vdots \\\\
y_{n-1}' = y_n \\\\
y_n' = f(t, y_1, y_2, \ldots, y_n)
\end{cases}$$

</div>
</div>

> 经历以上的降阶后，我们将重心放在如何数值求解**一阶**常微分方程的**初值问题**之上。

## 常用方法及其与数值积分公式的关系

在数值求解常微分方程初值问题时，常用的方法包括向前的欧拉法、向后的欧拉法、Crank-Nicolson 方法和经典的四阶龙格-库塔法。这些方法与数值积分公式密切相关。

> - **向前的欧拉法**：对应于矩形积分法，使用当前点的斜率来估计下一个点的值。
> - **向后的欧拉法**：对应于矩形积分法，使用下一个点的斜率来估计下一个点的值。
> - **Crank-Nicolson 方法**：对应于梯形积分法的隐式形式，结合了向前和向后的欧拉法的优点。
> - **四阶龙格-库塔法**：对应于 Simpson 积分法，使用多个斜率的加权平均值来估计下一个点的值。
> 
> 这些方法通过不同的方式近似积分，从而实现对常微分方程初值问题的数值求解。

## 向前的欧拉法（Forward Euler）

向前的欧拉法是一种显式数值方法，用于求解常微分方程初值问题。其基本思想是使用当前点的斜率来估计下一个点的值。利用从 $t_i$ 积分到 $t_{i+1}$ 的初值 ODE 给出一个数值积分公式的例子：

$$y_{i+1}-y_{i}=\int_{t_{i}}^{t_{i+1}} f(t, y(t)) d t \approx \Delta t f\left(t_{i}, y_{i}\right)$$

于是，给定常微分方程初值问题：
$$y' = f(t, y), \quad y(t_0) = y_0,$$
我们选择一个时间步长 $\Delta t$，并在离散时间点 $t_n = t_0 + n\Delta t$ 处近似解 $y_n \approx y(t_n)$。

所以我们得到向前的欧拉法的迭代公式为：
$$y_{n+1} = y_n + \Delta t f(t_n, y_n).$$

## 向后的欧拉法（Backward Euler）

向后的欧拉法是一种隐式数值方法，用于求解常微分方程初值问题。其基本思想是使用下一个点的斜率来估计下一个点的值。利用从 $t_i$ 积分到 $t_{i+1}$ 的初值 ODE 给出一个数值积分公式的例子：
$$y_{i+1}-y_{i}=\int_{t_{i}}^{t_{i+1}} f(t, y(t)) d t \approx \Delta t f\left(t_{i+1}, y_{i+1}\right)$$
于是，给定常微分方程初值问题：$y' = f(t, y)$，$y(t_0) = y_0,$

我们选择一个时间步长 $\Delta t$，并在离散时间点 $t_n = t_0 + n\Delta t$ 处近似解 $y_n \approx y(t_n)$。
所以我们得到向后的欧拉法的迭代公式为：
$$y_{n+1} = y_n + \Delta t f(t_{n+1}, y_{n+1}).$$
由于 $y_{n+1}$ 出现在等式的右侧，这使得向后的欧拉法成为隐式方法，通常需要通过迭代或数值求解方法来计算 $y_{n+1}$。

## Crank-Nicolson 方法

Crank-Nicolson 方法是一种隐式数值方法，用于求解常微分方程初值问题。它结合了向前和向后的欧拉法的优点，通过使用当前点和下一个点的斜率的平均值来估计下一个点的值。利用从 $t_i$ 积分到 $t_{i+1}$ 的初值 ODE 给出一个数值积分公式的例子：

$$
\begin{aligned}
y_{i+1}-y_{i}=\int_{t_{i}}^{t_{i+1}}f(t, y(t))dt \approx \frac{\Delta t}{2}\left[f\left(t_{i}, y_{i}\right)+f\left(t_{i+1}, y_{i+1}\right)\right]\\\\
\approx \Delta t f(\frac{t_{i+1}+t_{i}}{2}, \frac{y_{i+1}+y_{i}}{2})
\end{aligned}
$$

于是，给定常微分方程初值问题：$y' = f(t, y), \quad y(t_0) = y_0,$

我们选择一个时间步长 $\Delta t$，并在离散时间点 $t_n = t_0 + n\Delta t$ 处近似解 $y_n \approx y(t_n)$。
所以我们得到该方法的迭代公式为：
$$y_{n+1} = y_n + \Delta t f\left(\frac{t_{n+1}+t_n}{2}, \frac{y_{n+1}+y_n}{2}\right).$$

## The Simson's Rule

Simson 积分法是一种数值积分方法，用于近似计算定积分。其基本思想是通过使用二次插值多项式来近似被积函数，从而提高积分的精度。

对于区间 $[a, b]$ 上的函数 $f(t)$，Simson 积分法的公式为：
$$\int_{a}^{b} f(t) d t \approx \frac{b-a}{6}\left[f(a) + 4f\left(\frac{a+b}{2}\right) + f(b)\right].$$

在 RK4 方法中，我们使用类似的思想，通过计算多个斜率的加权平均值来估计下一个点的值，从而实现对常微分方程初值问题的数值求解。

## 四阶龙格-库塔法（RK4）

四阶龙格-库塔法（RK4）是一种显式数值方法，用于求解常微分方程初值问题。其基本思想是通过计算多个斜率的加权平均值来估计下一个点的值。利用从 $t_i$ 积分到 $t_{i+1}$ 的初值 ODE 给出一个数值积分公式的例子：
$$y_{i+1}-y_{i}=\int_{t_{i}}^{t_{i+1}} f(t, y(t)) d t \approx \frac{\Delta t}{6}\left[K_1 + 2K_2 + 2K_3 + K_4\right]$$
其中，
$$\begin{aligned}
K_1 & = f(t_i, y_i), \\\\
K_2 & = f\left(t_i + \frac{\Delta t}{2}, y_i + \frac{\Delta t}{2} K_1\right), \\\\
K_3 & = f\left(t_i + \frac{\Delta t}{2}, y_i + \frac{\Delta t}{2} K_2\right), \\\\
K_4 & = f(t_i + \Delta t, y_i + \Delta t K_3). 
\end{aligned}$$

## 四阶龙格-库塔法（RK4）

于是，给定常微分方程初值问题：
$$y' = f(t, y), \quad y(t_0) = y_0,$$
我们选择一个时间步长 $\Delta t$，并在离散时间点 $t_n = t_0 + n\Delta t$ 处近似解 $y_n \approx y(t_n)$。
所以我们得到四阶龙格-库塔法的迭代公式为：
$$y_{n+1} = y_n + \frac{\Delta t}{6}\left[K_1 + 2K_2 + 2K_3 + K_4\right].$$

## 隐式求解器与显式求解器的速度比较

对于每一步：

- **向前的欧拉法**：计算 $\displaystyle y_{n+1}=y_n+\Delta{t}f(t_n, y_n)$，只需简单计算函数值，速度较快。
- **向后的欧拉法**：计算 $\displaystyle y_{n+1}=y_n+\Delta{t}f(t_{n+1}, y_{n+1})$，需要求解隐式方程，速度较慢。对于求解隐式方程，通常需要使用迭代方法（如牛顿法）或数值求解器（如 MATLAB 的  `fsolve`  函数）。
  - 比如：
    - STEP-1: 定义一个函数 $\displaystyle g(x) = x - y_n - \Delta t f(t_{n+1}, x)$；
    - STEP-2: 迭代变为 $\displaystyle x_{k+1} = x_k - \frac{g(x_k)}{g'(x_k)}=\frac{x_k-y_i-\Delta t f(t_{n+1}, x_k)}{1-\Delta t f_{x}(t_{n+1}, x_k)}$，直到收敛为止。
    
    （在高维的情况下，$g'(x_k)$ 变为雅可比矩阵）

## 隐式求解器与显式求解器的速度比较

对于每一步：

- **Crank-Nicolson 方法**：计算 $\displaystyle y_{n+1}=y_n+\Delta{t}{f}\left(\frac{t_{n+1}+t_n}{2},\frac{y_{n+1}+y_n}{2}\right)$，同样需要求解隐式方程，速度较慢。
- **四阶龙格-库塔法**：计算 $\displaystyle y_{n+1}=y_n+\frac{\Delta t}{6}\left[K_1+2K_2+2K_3+K_4\right]$，需要多次计算函数值，但不需要求解隐式方程，速度较快。

## 稳定性比较

<span style="color:blue;">求解器“稳定性”的概念是基于这样的原则：如果 ODE 的真解不表现出指数增长，那么数值方法也不应该表现出指数增长 。</span>

在文中，选用了线性测试方程 $y' = \lambda y$ 来分析数值方法的稳定性。这个方程的解析解为 $y(t) = y_0 e^{\lambda t}$。指数增长发生在 $\text{Re}(\lambda) > 0$。当 $\text{Re}(\lambda) < 0$ 时，解趋于零，表示系统是稳定的。

> 通过对数值解是否会在 $\text{Re}(\lambda) < 0$ 时表现出指数增长来评估数值方法的稳定性。

## 稳定性比较： 向前的欧拉法（Forward Euler）

对于向前的欧拉法，应用于测试方程 $y' = \lambda y$，其迭代公式为：

$$y_{n+1} = y_n + \Delta t \lambda y_n = (1 + \Delta t \lambda) y_n.$$

通过递归，我们得到：

$$y_n = (1 + \Delta t \lambda)^n y_0.$$

为了确保数值解不表现出指数增长，我们需要满足以下条件：

$$|1 + \Delta t \lambda| \leq 1.$$

这个条件定义了向前的欧拉法的稳定区域。也就是说，如果 $\lambda$ 越大，则需要更小的时间步长 $\Delta t$ 来保持稳定性。这里的 $\lambda$ 表示局部线性化后得到的特征值，它决定了系统解的瞬态行为或局部增长/衰减率。

## 稳定性比较： 向后的欧拉法（Backward Euler）

对于向后的欧拉法，应用于测试方程 $y' = \lambda y$，其迭代公式为：
$$y_{n+1} = y_n + \Delta t \lambda y_{n+1}.$$
通过重排，我们得到：
$$y_{n+1} = \frac{1}{1 - \Delta t \lambda} y_n.$$
通过递归，我们得到：
$$y_n = \left(\frac{1}{1 - \Delta t \lambda}\right)^n y_0.$$

</div>
<div class="col">

为了确保数值解不表现出指数增长，我们需要满足以下条件：
$$\left|\frac{1}{1 - \Delta t \lambda}\right| \leq 1.$$
这个条件定义了向后的欧拉法的稳定区域。实际上，向后的欧拉法对于所有 $\text{Re}(\lambda) < 0$ 都是稳定的，因此它被称为**A-稳定**方法。

</div>

## 稳定性比较： Crank-Nicolson 方法

对于 Crank-Nicolson 方法，应用于测试方程 $y' = \lambda y$，其迭代公式为：
$$y_{n+1} = y_n + \frac{\Delta t}{2} \lambda (y_n + y_{n+1}).$$
通过重排，我们得到：
$$y_{n+1} = \frac{1 + \frac{\Delta t}{2} \lambda}{1 - \frac{\Delta t}{2} \lambda} y_n.$$
通过递归，我们得到：
$$y_n = \left(\frac{1 + \frac{\Delta t}{2} \lambda}{1 - \frac{\Delta t}{2} \lambda}\right)^n y_0.$$

</div>
<div class="col">

为了确保数值解不表现出指数增长，我们需要满足以下条件：
$$\left|\frac{1 + \frac{\Delta t}{2} \lambda}{1 - \frac{\Delta t}{2} \lambda}\right| \leq 1.$$
这个条件定义了 Crank-Nicolson 方法的稳定区域。Crank-Nicolson 方法在 $\text{Re}(\lambda) < 0$ 时也是稳定的，但其稳定性区域较向后的欧拉法略小。

</div>

## 稳定性比较： 四阶龙格-库塔法（RK4）

从计算所用的公式中可以预料到，分析 RK4 的稳定性比分析之前的那些方法要复杂得多。因此，书中直接给出了结果：RK4 的稳定性要求：

$$\left|1 + \lambda \Delta t + \frac{1}{2}(\lambda \Delta t)^2 + \frac{1}{6}(\lambda \Delta t)^3 + \frac{1}{24}(\lambda \Delta t)^4\right| \leq 1.$$

这最终得出与向前欧拉法类似的**条件稳定性**结果。

## ODE 求解器的比较

书中编写了MATLAB代码来验证各种ODE求解器的精度。通过比较数值解与解析解之间的误差，观察不同时间步长下的误差变化情况。

我们也来实践一下。

依旧使用书中的测试方程，也就是$y' = -0.8y^{1.5} + 20000(1 - e^{-3t})$，初始条件 $y(0) = 2000$。

| Method / N         |    21    |    51    |   101    |   501    |
| :----------------- | :------: | :------: | :------: | :------: |
| **Forward Euler**  | 0.000337 | 0.000559 | 0.001105 | 0.005142 |
| **Backward Euler** | 0.018156 | 0.034840 | 0.068626 | 0.314986 |
| **Crank-Nicolson** | 0.015286 | 0.036128 | 0.071233 | 0.329998 |
| **RK4**            | 0.000780 | 0.001624 | 0.003293 | 0.015901 |

在下一页，我们可以看到这些结果的图形表示。

![FIG: The Results](./CODE/resultofCODE10.png)

</div>
</div>

## ODE 求解器的准确度（Accuracy）

| Method \ E     | 21       | 51       | 101      | 501      |
| -------------- | -------- | -------- | -------- | -------- |
| Forward Euler  | 0.120458 | 0.046101 | 0.022736 | 0.004498 |
| Backward Euler | 0.104925 | 0.043678 | 0.022129 | 0.004473 |
| Crank-Nicolson | 0.003172 | 0.000507 | 0.000127 | 0.000005 |
| RK4            | 0.000007 | 0.000000 | 0.000000 | 0.000000 |

我们可以看出：
- 向前欧拉法和向后欧拉法的误差随着时间步长的减小而线性减少，符合一阶方法的预期。
- Crank-Nicolson 方法的误差随着时间步长的减小而二次减少，符合二阶方法的预期。
- 四阶龙格-库塔法的误差随着时间步长的减小而极速减少，符合预期。

使用的代码参见 [此处](/slides-scr/Scientific_Computing/CODE/code10.m) 。

## 总结

我们可以写出一个总结表格比较各种ODE求解器的速度、稳定性和准确度：

| 方法                | 速度 | 稳定性   | 准确度    |
| ------------------- | ---- | -------- | --------- |
| 向前欧拉法          | 快   | 条件稳定 | 低（1阶） |
| 向后欧拉法          | 慢   | A-稳定   | 低（1阶） |
| Crank-Nicolson 方法 | 慢   | A-稳定   | 中（2阶） |
| 四阶龙格-库塔法     | 快   | 条件稳定 | 高（4阶） |

根据以上比较，可以选择适合特定问题的ODE求解器。

# 敬请老师和同学们批评指正！

<hr/>

参考资料：[Scientific Computing: For Scientists and Engineers](https://ieeexplore.ieee.org/servlet/opac?bknumber=10783660)

我的部分联系方式：[hsingyu.yen@icloud.com](mailto:hsingyu.yen@icloud.com)

</div>
</div>