贡献者: addis
我们先来尝试用欧拉法解一阶微分方程
\begin{equation}
y'(t) = y~.
\end{equation}
令初始条件为 $y(0) = 1$。令步长为 $h = 0.5$,步数为 $5$,结果如
图 1 所示(代码见文章最后)。
图 1:欧拉法数值解(蓝)和解析解(红)
我们知道该方程的解析解为 $y = \mathrm{e} ^{t}$。对比数值解和解析解,不难分析出误差产生的原因:我们仅用每段步长区间左端的导数预测整个区间的曲线增量。如果我们能利用每个区间中点的导数计算整个区间的增量,这个预测将会比欧拉法更精确。
考虑微分方程 $y'(t) = f(y, t)$ 在区间 $[t_n, t_{n+1}]$ 的曲线,若我们已知区间左端的函数值为 $y_n$,我们可以先用微分近似估计曲线中点的函数值为
\begin{equation}
y \left(t_n + \frac h2 \right) = y_n + \frac h2 f(y_n, t_n)~.
\end{equation}
然后再求出这个近似中点的导数为
\begin{equation}
y' \left(t_n + \frac h2 \right) = f \left[y_n + \frac h2 f(t_n, y_n), t_n + \frac h2 \right] ~.
\end{equation}
最后我们利用这个导数估算该区间的曲线增量
\begin{equation}
y_{n+1} = y_n + hy' \left(t_n + \frac h2 \right) = y_n + h f \left[y_n + \frac h2 f(t_n, y_n), t_n + \frac h2 \right] ~.
\end{equation}
这就是解常微分方程的
中点法。
我们再来用中点法取同样的步长计算式 1 ,结果如图 2 所示。
图 2:中点法数值解(蓝)和解析解(红)
可见虽然中点法每一步的计算过程比欧拉法稍微复杂一些,但精度却大大地提高了。
中点法同样适用于微分方程组,例如对于常微分方程组
\begin{equation}
\begin{cases}
x'(t) = f(x, y, t)\\
y'(t) = g(x, y, t)
\end{cases}~.
\end{equation}
首先计算近似中点为
\begin{equation}
\begin{cases}
x_{n+1/2} = x_n + \frac{h}{2} f(x_n, y_n, t_n)\\
y_{n+1/2} = y_n + \frac{h}{2} g(x_n, y_n, t_n)
\end{cases}~.
\end{equation}
然后有
\begin{equation}
\begin{cases}
x_{n+1} = x_n + h f \left(x_{n+1/2}, y_{n+1/2}, t_n + h/2 \right) \\
y_{n+1} = y_n + h g \left(x_{n+1/2}, y_{n+1/2}, t_n + h/2 \right)
\end{cases}~.
\end{equation}
代码 1:odeMid.m
% 设置参数
N = 6;
h = 0.5;
t = linspace(0, (N-1)*h, N); % 自变量
t0 = linspace(0, (N-1)*h, 100); % 用于画图
% 欧拉法
y = zeros(1,N); % 预赋值
y(1) = 1; % 初值
for ii = 1:N-1
y(ii+1) = y(ii) + h*y(ii);
end
% 画图
figure;
plot(x,y,'+-');
hold on;
plot(t0, exp(t0));
%中点法
y = zeros(1,N);
y(1) = 1;
for ii = 1:N-1
y(ii+1) = y(ii) + h*(y(ii) + 0.5*h*y(ii));
end
% 画图
figure;
plot(x,y,'+-');
hold on;
plot(x0, exp(x0));