龙格库塔法

参考视频

# 公式

$$y(t + \Delta t) = y(t) + y'(t) \Delta t + \frac{1}{2} y''(\Delta t)^2 + \cdots + \frac{1}{n!} y^{(n+1)}(t)(\Delta t)^n + R_n$$

其中，余项$R_n$ 表示为：

$$R_n = \frac{1}{(n + 1)!} f^{(n + 1)}(t + \theta t)(\Delta t)^{n + 1} \rightarrow O((\Delta t)^{n + 1})$$

# 示例

若$y = e^t$ 在$y = 0$ 处展开得：

$$y(0 + \Delta t) = y(\Delta t) = 1 + \Delta t + \frac{1}{2} (\Delta t)^2 + \frac{1}{3!} (\Delta t)^3 + \cdots$$

也可以写成：

$$y(t) = 1 + t + \frac{1}{2} t^2 + \frac{1}{3!} t^3 + \cdots$$

# 一级近似

用$y = e^t$ 的一级近似$1 + t$ 和原函数做对比：

clear

t = -1:0.001:1;

y = exp(t);

plot(t, y)

hold on

plot(t, 1 + t)

几何意义：在展开点处，一级近似的斜率与原函数相同。

# 二级近似

用$y = e^2$ 的二级近似$1 + t + \frac{1}{2} t^2$ 和原函数做对比：

clear

t = -1:0.001:1;

y = exp(t);

plot(t, y)

hold on

plot(t, 1 + t + 1/2 * t.^2)

几何意义：在展开点处的二级近似的一阶 / 二阶导数与原函数相同。

之后也可画出三级 / 四级近似和原函数的对比图，发现拟合程度越来越高，代表误差越来越小

# 欧拉法

$$\frac{dy}{dx} = f(y, t)$$

$$y(t_n) = y_n$$

$$y_{n+1} = y_n + f(y_n, t_n) \Delta t$$

对$y(t_n)$ 在$t_n$ 处做泰勒展开：

$$y(t_n + 1) = y(t_n + \Delta t) = y(t_n) + y'(t_n) \Delta t + O((\Delta t)^2) = y(t_n) + f(y_n, t_n) \Delta t + O((\Delta t)^2)$$

对比可发现，欧拉法的误差量级为$O((\Delta t)^2)$，这是一步产生的误差量级。

比如，从$0$ 模拟到$1$，共有$\frac{1}{\Delta t}$ 步；每一步都产生新的误差，这些误差会累积，产生的误差量级为：

$$O((\Delta t)^2) \cdot \frac{1}{\Delta t} = O(\Delta t)$$

因此，这称为一阶误差。

如果一步的误差是$O((\Delta t)^3)$，则从$0$ 到$1$ 产生的误差量级为：

$$O((\Delta t)^3) \cdot \frac{1}{\Delta t} = O((\Delta t)^2)$$

因此，这称为二阶误差

# 利用泰勒展开

对$\frac{d^2y}{dt^2} = \frac{df(y,t)}{dt}$ 进行展开：

$$\frac{d^2y}{dt^2} = \frac{df(y, t)}{dt} = \frac{\partial f(y, t)}{\partial y} \cdot \frac{dy}{dt} + \frac{\partial f(y, t)}{\partial t}$$

可以进一步表示为：

$$\frac{d^2y}{dt^2} = \frac{\partial f(y, t)}{\partial y} \cdot f(y, t) + \frac{\partial f(y, t)}{\partial t}$$

结合欧拉法：

$$y_{n+1} = y_n + f(y_n, t_n) \Delta t + \frac{1}{2} \left( \frac{\partial f(y, t)}{\partial y} \cdot f(y, t) + \frac{\partial f(y, t)}{\partial t} \right) (\Delta t)^2$$

# 实际例子

$\frac{dy}{dt} = y, \; y(0) = y_0$

$$y(t) = y_0 e^t$$

通过二阶近似：

$$y_{n+1} = y_n + y_n \Delta t + \frac{1}{2} y_n (\Delta t)^2$$

clear

T = 1;

N = 4;

dt = T/N;

t = 0:dt:T;

t1 = 0:0.001:1;

y = zeros(1, length(t));

y(1) = 8;

for n = 1:length(t)-1  

    y(n+1) = y(n) + y(n) * dt;

end

plot(t, y, t1, 8 * exp(t1))

欧拉法最终形成的误差和$\Delta t$ 成正比，可以试着

把 N 改成 40

把 N 改成 400

下面试试二阶近似

$$y_{n+1} = y_n + y_n\Delta t + \frac{1}{2} y_n \Delta t^2$$

clear

T=1;

N=4;

dt=T/N;

t=0:dt:T;

t1=0:0.001:1;

y=zeros(1,length(t));

y(1)=8;

y1=zeros(1,length(t));

y1(1)=8;

for n=1:length(t)-1  

    y(n+1)=y(n)+y(n)*dt;

    y1(n+1)=y1(n)+y1(n)*dt+1/2*y1(n)*dt^2;

end

plot(t,y,t,y1,t1,8*exp(t1))

可以看出，在 N 一样时，二阶近似比欧拉法的精度高上不少

# 改进欧拉法

改进欧拉法得到：

$$y_{n+1} = y_n + \frac{1}{2} \left( f(y_n, t_n) + f(y_{n+1}, t_{n+1}) \right) \Delta t$$

其中$f(y_{n+1}, t_{n+1})$ 用欧拉法公式计算，得到：

$$y_{n+1} = y_n + \frac{1}{2} \left( f(y_n, t_n) + f\left(y_n + f(y_n, t_n) \Delta t, t_n + \Delta t \right) \right) \Delta t$$

现在我们要证明这个式子是二阶精度的。

# 泰勒展开

对$f(y + \Delta y, t + \Delta t)$ 做泰勒展开：

$$f(y + \Delta y, t + \Delta t) = f(y, t) + \frac{\partial f}{\partial y} \Delta y + \frac{\partial f}{\partial t} \Delta t + O((\Delta t)^2)$$

替换$\Delta y = f(y, t) \Delta t$，得到：

$$f(y + \Delta y, t + \Delta t) = f(y_n, t_n) + \frac{\partial f}{\partial y} (f(y_n, t) \Delta t) + \frac{\partial f}{\partial t} \Delta t + O((\Delta t)^2)$$

# 代入改进欧拉法公式

将展开结果代入改进欧拉法：

$$y_{n+1} = y_n + \frac{1}{2} \left( f(y_n, t_n) + f\left(y_n + f(y_n, t_n) \Delta t, t_n + \Delta t \right) \right) \Delta t$$

=

$$y_n + \frac{1}{2} \left( f(y_n, t_n) + f(y_n, t_n) + \frac{\partial f}{\partial y} (f(y_n, t) \Delta t) + \frac{\partial f}{\partial t} \Delta t + O((\Delta t)^2) \right) \Delta t$$

=

$$y_n + f(y_n, t_n) \Delta t + \frac{1}{2} \left( \frac{\partial f(y_n, t_n)}{\partial y} f(y_n, t) + \frac{\partial f(y_n, t)}{\partial t} \right) (\Delta t)^2 + O((\Delta t)^3)$$

# 误差分析

和$y_n$ 在$t_n$ 处的泰勒展开公式做对比，可以看出改进的欧拉法的误差量级也是$O((\Delta t)^3)$。

# 以$k$ 表示

定义：

$$k_1 = f(y_n, t_n)$$

$$k_2 = f\left(y_n + k_1 \Delta t, t_n + \Delta t\right)$$

最终得到：

$$y_{n+1} = y_n + \frac{1}{2} (k_1 + k_2) \Delta t$$