微分方程數值解作業20250307-优快云博客

本文链接：https://blog.youkuaiyun.com/weixin_46616813/article/details/146123669

1. 求解初值问题

$\left\{\begin{array}{l} \frac{d u}{dt}=u-\frac{2 t}{u}, 0<t<1, \\ u(0)=1 \end{array}\right.$

该问题的精确解是 $u=\sqrt{1 + 2t}$ 。请分别用显式Euler法、梯形公式和预估 - 校正改进Euler公式求解上述初值问题，时间步长取 $\Delta t = 0.1$ ，列表给出精确解和各算法的数值解与计算误差。(计算结果保留 $6$ 位有效数字，梯形公式的收敛参数取为 $\varepsilon = 10^{-6}$ )。

- 显式Euler法：

对于初值问题 $\frac{du}{dt}=f(t, u)$ ， $u(t_0)=u_0$ ，显式Euler法的迭代公式为 $u_{n + 1}=u_n+h f(t_n, u_n)$ ，这里 $f(t, u)=u-\frac{2t}{u}$ ， $h = 0.1$ ， $t_0 = 0$ ， $u_0 = 1$ 。

- 梯形公式：

迭代公式为 $u_{n+1}=u_n+\frac{h}{2}[f(t_n, u_n)+f(t_{n + 1}, u_{n+1})]$ ，这是一个隐式公式，需要迭代求解。给定收敛参数 $\varepsilon = 10^{-6}$ ，一般采用迭代法，如先猜测 $u_{n+1}^0=u_n+h f(t_n, u_n)$ ，然后迭代 $u_{n+1}^{k + 1}=u_n+\frac{h}{2}[f(t_n, u_n)+f(t_{n + 1}, u_{n+1}^k)]$ ，直到 $|u_{n+1}^{k + 1}-u_{n+1}^k|<\varepsilon$ 。

- 预估 - 校正改进Euler公式：

先预估 $\overline{u}_{n+1}=u_n+h f(t_n, u_n)$ ，再校正 $u_{n+1}=u_n+\frac{h}{2}[f(t_n, u_n)+f(t_{n + 1},\overline{u}_{n+1})]$ 。

Matlab代碼：

% 显式Euler法
t = 0:0.1:1;
u_euler = zeros(size(t));
u_euler(1) = 1;
for n = 1:length(t)-1
u_euler(n + 1)=u_euler(n)+0.1*(u_euler(n)-2*t(n)/u_euler(n));
end

% 梯形公式
u_trapezoid = zeros(size(t));
u_trapezoid(1) = 1;
epsilon = 1e-6;
for n = 1:length(t)-1
u0 = u_trapezoid(n)+0.1*(u_trapezoid(n)-2*t(n)/u_trapezoid(n));
u = u0;
while true
u_new = u_trapezoid(n)+0.1/2*((u_trapezoid(n)-2*t(n)/u_trapezoid(n))+(u-2*t(n + 1)/u));
if abs(u_new - u)<epsilon
u_trapezoid(n + 1)=u_new;
break;
end
u = u_new;
end
end

% 预估 - 校正改进Euler公式
u_improved_euler = zeros(size(t));
u_improved_euler(1) = 1;
for n = 1:length(t)-1
u_bar = u_improved_euler(n)+0.1*(u_improved_euler(n)-2*t(n)/u_improved_euler(n));
u_improved_euler(n + 1)=u_improved_euler(n)+0.1/2*((u_improved_euler(n)-2*t(n)/u_improved_euler(n))+(u_bar-2*t(n + 1)/u_bar));
end

% 精确解
u_exact = sqrt(1 + 2*t);

% 计算误差
error_euler = abs(u_exact - u_euler);
error_trapezoid = abs(u_exact - u_trapezoid);
error_improved_euler = abs(u_exact - u_improved_euler);

% 输出结果
fprintf('t\t精确解\t显式Euler法数值解\t显式Euler法误差\t梯形公式数值解\t梯形公式误差\t预估 - 校正改进Euler公式数值解\t预估 - 校正改进Euler公式误差\n');
for i = 1:length(t)
fprintf('%0.1f\t%0.6f\t%0.6f\t%0.6e\t%0.6f\t%0.6e\t%0.6f\t%0.6e\n',t(i),u_exact(i),u_euler(i),error_euler(i),u_trapezoid(i),error_trapezoid(i),u_improved_euler(i),error_improved_euler(i));
end

t	精确解	显式Euler法数值解	显式Euler法误差	梯形公式数值解	梯形公式误差	预估 - 校正改进Euler公式数值解	预估 - 校正改进Euler公式误差
0.0	1.000000	1.000000	0.000000e+00	1.000000	0.000000e+00	1.000000	0.000000e+00
0.1	1.095445	1.100000	4.554885e-03	1.095656	2.107736e-04	1.095909	4.639759e-04
0.2	1.183216	1.191818	8.602225e-03	1.183594	3.778269e-04	1.184097	8.806126e-04
0.3	1.264911	1.277438	1.252677e-02	1.265441	5.296563e-04	1.266201	1.290297e-03
0.4	1.341641	1.358213	1.657181e-02	1.342323	6.819121e-04	1.343360	1.719365e-03
0.5	1.414214	1.435133	2.091936e-02	1.415058	8.449289e-04	1.416402	2.188366e-03
0.6	1.483240	1.508966	2.572656e-02	1.484267	1.026866e-03	1.485956	2.715905e-03
0.7	1.549193	1.580338	3.114490e-02	1.550429	1.235219e-03	1.552514	3.320753e-03
0.8	1.612452	1.649783	3.733188e-02	1.613929	1.477669e-03	1.616475	4.023233e-03
0.9	1.673320	1.717779	4.445929e-02	1.675083	1.762649e-03	1.678166	4.846311e-03
1.0	1.732051	1.784771	5.272002e-02	1.734151	2.099794e-03	1.737867	5.816593e-03