我们先来尝试用欧拉法解一阶微分方程
\begin{equation}
y'(t) = y
\end{equation}
令初始条件为 $y(0) = 1$.令步长为 $h = 0.5$, 步数为 $5$, 结果如
图 1 所示(代码见词条最后).
图1:欧拉法数值解(蓝)和解析解(红)
由“”我们知道该方程的解析解为 $y = \E^{t}$. 对比数值解和解析解, 不难分析出误差产生的原因: 我们仅用每段步长区间左端的导数预测整个区间的曲线增量. 如果我们能利用每个区间中点的导数计算整个区间的增量, 这个预测将会比欧拉法更精确.
考虑微分方程 $y'(t) = f(y, t)$ 在区间 $[t_n, t_{n+1}]$ 的曲线, 若我们已知区间左端的函数值为 $y_n$, 我们可以先用微分近似估计曲线中点的函数值为
\begin{equation}
y \qtyRound{t_n + \frac h2} = y_n + \frac h2 f(y_n, t_n)
\end{equation}
然后再求出这个近似中点的导数为
\begin{equation}
y' \qtyRound{t_n + \frac h2} = f \qtySquare{y_n + \frac h2 f(t_n, y_n), t_n + \frac h2}
\end{equation}
最后我们利用这个导数估算该区间的曲线增量
\begin{equation}
y_{n+1} = y_n + hy' \qtyRound{t_n + \frac h2} = y_n + h f \qtySquare{y_n + \frac h2 f(t_n, y_n), t_n + \frac h2}
\end{equation}
这就是解常微分方程的
中点法.
我们再来用中点法取同样的步长计算式 1 , 结果如图 2 所示.
图2:中点法数值解(蓝)和解析解(红)
可见虽然中点法每一步的计算过程比欧拉法稍微复杂一些, 但精度却大大地提高了.
中点法同样适用于微分方程组, 例如对于常微分方程组
\begin{equation}\leftgroup{
x'(t) &= f(x, y, t)\\
y'(t) &= g(x, y, t)
}\end{equation}
首先计算近似中点为
\begin{equation}\leftgroup{
x_{n+1/2} &= x_n + h f(x_n, y_n, t_n)\\
y_{n+1/2} &= y_n + h g(x_n, y_n, t_n)
}\end{equation}
然后有
\begin{equation}
\leftgroup{
x_{n+1} &= x_n + h f \qtyRound{x_{n+1/2}, y_{n+1/2}, \frac h2}\\
y_{n+1} &= y_n + h g \qtyRound{x_{n+1/2}, y_{n+1/2}, \frac h2}
}\end{equation}
odeMid.m
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31
|
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));
|
致读者: 小时物理百科一直以来坚持所有内容免费且不做广告,这导致我们处于日渐严重的亏损状态。长此以往很可能会最终导致我们不得不选择商业化,例如大量广告,内容付费,会员制,甚至被收购。因此,我们鼓起勇气在此请求广大读者
热心捐款,使网站得以健康发展。如果看到这条信息的每位读者能慷慨捐助 10 元,我们几天内就能脱离亏损状态,并保证网站能在接下来的一整年里向所有读者继续免费提供优质内容。感谢您的支持。
—— 小时(项目创始人)