图

中点法解常微分方程（组）

预备知识　常微分方程（组）的数值解

　　 我们先来尝试用欧拉法解一阶微分方程

\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} \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 + h f(x_n, y_n, t_n)\\ y_{n+1/2} = y_n + h g(x_n, y_n, t_n) \end{cases} \end{equation}
然后有
\begin{equation} \begin{cases} x_{n+1} = x_n + h f \qtyRound{x_{n+1/2}, y_{n+1/2}, h/2}\\ y_{n+1} = y_n + h g \qtyRound{x_{n+1/2}, y_{n+1/2}, h/2} \end{cases} \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 元，我们几天内就能脱离亏损状态，并保证网站能在接下来的一整年里向所有读者继续免费提供优质内容。感谢您的支持。
—— 小时（项目创始人）

编辑词条（需要权限）　返回目录　返回主页　捐助项目　© 小时物理百科 保留一切权利