此资源是我自己以前写的一篇随笔（word格式），对欧拉法与龙格库塔法进行了讲解，并利用matlab进行2~4阶龙格库塔法解常微分方程的仿真，附带详细注释，并输出不同解法下的对比结果，对学习龙格库塔法和matlab的新手会有帮助
《欧拉法与龙格库塔法解常微分方程——Matlab实现》
常微分方程在科学计算中扮演着至关重要的角色，它广泛应用于物理学、工程学、生物学等多个领域。解决这类问题的方法多种多样，其中欧拉法和龙格库塔法是最常见的数值解法。本文将详细介绍这两种方法，并结合Matlab代码进行实际应用。
欧拉法是一种简单的数值积分方法，用于求解初值问题的一阶常微分方程。基本思想是从已知的初始点出发，通过直线近似来逼近函数的真实曲线。欧拉法的公式如下：
\[ y_{n+1} = y_n + h f(x_n, y_n) \]
这里的 \( h \) 是步长，\( f(x, y) \) 是微分方程的右侧函数，\( (x_n, y_n) \) 是当前点的坐标。欧拉法具有一阶精度，但其局部阶段误差与步长 \( h \) 成正比，是二阶无穷小量。
改进欧拉法，也称为Heun法，通过引入预报值和校正值来提高精度。它的迭代公式如下：
\[ y_{n+1} = y_n + \frac{h}{2}(f(x_n, y_n) + f(x_n + h, y_n + hf(x_n, y_n))) \]
这种方法的局部阶段误差是 \( h^3 \) 的量级，因此具有二阶精度。
龙格库塔法则是欧拉法和改进欧拉法的推广，它通过在每个时间步内考虑更多的函数值来提高精度。以 \( m \) 个点为例，龙格库塔法的迭代公式可表示为：
\[ y_{n+1} = y_n + h \sum_{j=1}^{m} \lambda_j K_j \]
其中 \( K_j = f(\xi_j, y_n + h \sum_{i=1}^{j-1} \beta_{ij} K_i) \)，\( \xi_j \) 在 \( [x_n, x_{n+1}] \) 区间内，\( \alpha_j, \beta_{ij}, \lambda_j \) 是特定于算法阶数的系数。二阶、三阶和四阶的龙格库塔公式分别如下：
- 二阶龙格库塔：\( \lambda_1 = \lambda_2 = \frac{1}{2} \)，\( \alpha_1 = \alpha_2 = 1 \)，\( \beta_{11} = 0 \)，\( \beta_{21} = \frac{1}{2} \)
- 三阶龙格库塔：\( \lambda_1 = \frac{1}{6} \)，\( \lambda_2 = \frac{2}{3} \)，\( \lambda_3 = \frac{1}{6} \)，\( \alpha_1 = 1 \)，\( \alpha_2 = \frac{1}{2} \)，\( \alpha_3 = 1 \)，\( \beta_{11} = 0 \)，\( \beta_{21} = \beta_{22} = \frac{1}{2} \)，\( \beta_{31} = \beta_{32} = \frac{1}{2} \)，\( \beta_{33} = 1 \)
- 四阶龙格库塔：\( \lambda_1 = \lambda_3 = \lambda_4 = \frac{1}{6} \)，\( \lambda_2 = \frac{1}{3} \)，\( \alpha_1 = 1 \)，\( \alpha_2 = \alpha_3 = \frac{1}{2} \)，\( \alpha_4 = 1 \)，\( \beta_{11} = 0 \)，\( \beta_{21} = \beta_{22} = \frac{1}{2} \)，\( \beta_{31} = \beta_{32} = \beta_{33} = \frac{1}{2} \)，\( \beta_{41} = \beta_{42} = \beta_{43} = 1 \)
Matlab提供了强大的工具来实现这些数值解法。以下代码展示了如何使用Matlab求解微分方程 \( y' = y\cos(x) \)，\( y(0) = 1 \) 的2、3、4阶龙格库塔解：
```matlab
syms x y; % 定义符号变量
f(x,y) = y*cos(x); % 微分方程导数
x0 = 0; y0 = 1; % 初始条件
h = 0.2; % 步长
x = x0:h:10; % x轴范围
y = exp(sin(x)); % 真实解
% 龙格库塔法求解
y2 = zeros(size(x)); % 二阶
y3 = zeros(size(x)); % 三阶
y4 = zeros(size(x)); % 四阶
y2(1) = y0;
y3(1) = y0;
y4(1) = y0;
for ii = 2:length(x)
% 二阶
K1 = f(x(ii-1), y2(ii-1));
K2 = f(x(ii-1)+h/2, y2(ii-1)+h*K1/2);
y2(ii) = y2(ii-1) + h*K2;
% 三阶
K1 = f(x(ii-1), y3(ii-1));
K2 = f(x(ii-1)+h/2, y3(ii-1)+h*K1/2);
K3 = f(x(ii-1)+h, y3(ii-1)+h*(K2*2-K1));
y3(ii) = y3(ii-1) + h*(K1+4*K2+K3)/6;
% 四阶
K1 = f(x(ii-1), y4(ii-1));
K2 = f(x(ii-1)+h/2, y4(ii-1)+h*K1/2);
K3 = f(x(ii-1)+h/2, y4(ii-1)+h*K2/2);
K4 = f(x(ii-1)+h, y4(ii-1)+h*K3);
y4(ii) = y4(ii-1) + h*(K1+2*K2+2*K3+K4)/6;
end
% 绘制解和误差图
figure
plot(x,y, x,y2, x,y3, x,y4)
title('龙格库塔求解值')
xlabel('x')
ylabel('y')
legend('原值', '2阶值', '3阶值', '4阶值')
figure
plot(x,y2-y, x,y3-y, x,y4-y)
title('龙格库塔误差')
xlabel('x')
ylabel('误差')
legend('2阶误差', '3阶误差', '4阶误差')
```
通过这段代码，我们可以直观地看到不同阶数的龙格库塔解与真实解的接近程度，以及随着阶数增加误差的减小。这为理解和掌握数值解法提供了实践基础，对于学习Matlab和数值计算的初学者来说尤其有价值。
评论2