欧拉法与龙格库塔法解常微分方程（附Matlab代码）

此资源是我自己以前写的一篇随笔（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和数值计算的初学者来说尤其有价值。