四阶龙格库塔法(RK4)与欧拉法在MATLAB中的实现与对比分析
四阶龙格库塔法(RK4)与欧拉法在MATLAB中的实现与对比分析
四阶龙格库塔法(RK4)与欧拉法在MATLAB中的实现与对比分析
1. 问题描述
我们考虑以下一阶线性微分方程: d y d x = − 2 y + 2 \frac{dy}{dx} = -2y + 2 dxdy=−2y+2
初始条件为 y ( 0 ) = 1 y(0) = 1 y(0)=1。
我们的微分方程是:
d
y
d
x
+
2
y
=
2
\frac{dy}{dx} + 2y = 2
dxdy+2y=2
这是一个线性微分方程,使用积分因子法求解:
- 积分因子: μ ( x ) = e ∫ 2 d x = e 2 x \mu(x) = e^{\int 2 dx} = e^{2x} μ(x)=e∫2dx=e2x
- 乘以积分因子:
e 2 x d y d x + 2 e 2 x y = 2 e 2 x e^{2x}\frac{dy}{dx} + 2e^{2x}y = 2e^{2x} e2xdxdy+2e2xy=2e2x - 左边为导数:
d d x ( e 2 x y ) = 2 e 2 x \frac{d}{dx}\left(e^{2x}y\right) = 2e^{2x} dxd(e2xy)=2e2x - 积分:
e 2 x y = ∫ 2 e 2 x d x = e 2 x + C e^{2x}y = \int 2e^{2x} dx = e^{2x} + C e2xy=∫2e2xdx=e2x+C - 求解 ( y ):
y ( x ) = 1 + C e − 2 x y(x) = 1 + Ce^{-2x} y(x)=1+Ce−2x
根据初始条件 y ( 0 ) = 1 y(0) = 1 y(0)=1:
1
=
1
+
C
⟹
C
=
0
1 = 1 + C \implies C = 0
1=1+C⟹C=0
所以解析解为:
y
(
x
)
=
1
y(x) = 1
y(x)=1
注意:由于初始条件
y
(
0
)
=
1
y(0) = 1
y(0)=1 是方程的稳态解,因此无论使用哪种数值方法,理论上数值解都应保持
y
=
1
y = 1
y=1。为了更好地展示数值方法的性能,我们将稍微修改初始条件为
y
(
0
)
=
2
y(0) = 2
y(0)=2,使得解析解为:
y
(
x
)
=
1
+
(
2
−
1
)
e
−
2
x
=
1
+
e
−
2
x
y(x) = 1 + (2 - 1)e^{-2x} = 1 + e^{-2x}
y(x)=1+(2−1)e−2x=1+e−2x
2. MATLAB代码实现
以下是MATLAB代码,包括欧拉法和四阶龙格库塔法的实现,以及两者与解析解的对比。
%% 清空环境
clear;
clc;
close all;
%% 定义微分方程
% dy/dx = -2y + 2
f = @(x, y) -2*y + 2;
%% 解析解
% y(x) = 1 + (y0 - 1)*exp(-2x)
y_exact = @(x, y0) 1 + (y0 - 1)*exp(-2*x);
%% 初始条件
x0 = 0; % 初始x值
y0 = 2; % 初始y值
xf = 5; % 结束x值
h = 0.1; % 步长
n = floor((xf - x0)/h); % 步数
%% 初始化结果存储
x_euler = zeros(1, n+1);
y_euler = zeros(1, n+1);
x_rk4 = zeros(1, n+1);
y_rk4 = zeros(1, n+1);
x_exact = linspace(x0, xf, n+1);
y_exact_vals = y_exact(x_exact, y0);
% 设置初始值
x_euler(1) = x0;
y_euler(1) = y0;
x_rk4(1) = x0;
y_rk4(1) = y0;
%% 欧拉法迭代过程
for i = 1:n
x_euler(i+1) = x_euler(i) + h;
y_euler(i+1) = y_euler(i) + h * f(x_euler(i), y_euler(i));
end
%% 四阶龙格库塔法(RK4)迭代过程
for i = 1:n
k1 = f(x_rk4(i), y_rk4(i));
k2 = f(x_rk4(i) + h/2, y_rk4(i) + h*k1/2);
k3 = f(x_rk4(i) + h/2, y_rk4(i) + h*k2/2);
k4 = f(x_rk4(i) + h, y_rk4(i) + h*k3);
y_rk4(i+1) = y_rk4(i) + (h/6)*(k1 + 2*k2 + 2*k3 + k4);
x_rk4(i+1) = x_rk4(i) + h;
end
%% 可视化结果
figure;
hold on;
plot(x_exact, y_exact_vals, 'k-', 'LineWidth', 2, 'DisplayName', '解析解');
plot(x_euler, y_euler, 'b--o', 'LineWidth', 1.5, 'MarkerSize', 4, 'DisplayName', '欧拉法');
plot(x_rk4, y_rk4, 'r-+', 'LineWidth', 1.5, 'MarkerSize', 4, 'DisplayName', '四阶RK4法');
title('欧拉法与四阶龙格库塔法(RK4)对比');
xlabel('x');
ylabel('y');
legend('Location', 'best');
grid on;
hold off;
%% 误差分析
% 计算欧拉法和RK4与解析解的误差
error_euler = abs(y_euler - y_exact_vals);
error_rk4 = abs(y_rk4 - y_exact_vals);
% 绘制误差曲线
figure;
hold on;
plot(x_exact, error_euler, 'b--', 'LineWidth', 1.5, 'DisplayName', '欧拉法误差');
plot(x_exact, error_rk4, 'r-', 'LineWidth', 1.5, 'DisplayName', '四阶RK4法误差');
title('欧拉法与四阶龙格库塔法(RK4)误差对比');
xlabel('x');
ylabel('误差');
legend('Location', 'best');
grid on;
hold off;
3. 代码详解
3.1 清空环境
clear;
clc;
close all;
clear
:清除工作区中的所有变量。clc
:清除命令窗口。close all
:关闭所有打开的图形窗口。
3.2 定义微分方程
f = @(x, y) -2*y + 2;
- 定义微分方程的右侧函数 ( f(x, y) = -2y + 2 )。
3.3 解析解
y_exact = @(x, y0) 1 + (y0 - 1)*exp(-2*x);
- 定义解析解函数,便于后续误差分析和可视化。
3.4 初始条件
x0 = 0; % 初始x值
y0 = 2; % 初始y值
xf = 5; % 结束x值
h = 0.1; % 步长
n = floor((xf - x0)/h); % 步数
- 设置初始和结束点、步长以及迭代次数。
3.5 初始化结果存储
x_euler = zeros(1, n+1);
y_euler = zeros(1, n+1);
x_rk4 = zeros(1, n+1);
y_rk4 = zeros(1, n+1);
x_exact = linspace(x0, xf, n+1);
y_exact_vals = y_exact(x_exact, y0);
% 设置初始值
x_euler(1) = x0;
y_euler(1) = y0;
x_rk4(1) = x0;
y_rk4(1) = y0;
- 预先分配数组以提高效率。
x_exact
和y_exact_vals
用于存储解析解的值。
3.6 欧拉法迭代过程
for i = 1:n
x_euler(i+1) = x_euler(i) + h;
y_euler(i+1) = y_euler(i) + h * f(x_euler(i), y_euler(i));
end
- 简单的欧拉法迭代公式:
y n + 1 = y n + h ⋅ f ( x n , y n ) y_{n+1} = y_n + h \cdot f(x_n, y_n) yn+1=yn+h⋅f(xn,yn)
3.7 四阶龙格库塔法(RK4)迭代过程
for i = 1:n
k1 = f(x_rk4(i), y_rk4(i));
k2 = f(x_rk4(i) + h/2, y_rk4(i) + h*k1/2);
k3 = f(x_rk4(i) + h/2, y_rk4(i) + h*k2/2);
k4 = f(x_rk4(i) + h, y_rk4(i) + h*k3);
y_rk4(i+1) = y_rk4(i) + (h/6)*(k1 + 2*k2 + 2*k3 + k4);
x_rk4(i+1) = x_rk4(i) + h;
end
- RK4的经典迭代公式:
k 1 = f ( x n , y n ) k 2 = f ( x n + h 2 , y n + h 2 k 1 ) k 3 = f ( x n + h 2 , y n + h 2 k 2 ) k 4 = f ( x n + h , y n + h k 3 ) y n + 1 = y n + h 6 ( k 1 + 2 k 2 + 2 k 3 + k 4 ) \begin{aligned} k_1 &= f(x_n, y_n) \\ k_2 &= f\left(x_n + \frac{h}{2}, y_n + \frac{h}{2}k_1\right) \\ k_3 &= f\left(x_n + \frac{h}{2}, y_n + \frac{h}{2}k_2\right) \\ k_4 &= f(x_n + h, y_n + h k_3) \\ y_{n+1} &= y_n + \frac{h}{6}(k_1 + 2k_2 + 2k_3 + k_4) \end{aligned} k1k2k3k4yn+1=f(xn,yn)=f(xn+2h,yn+2hk1)=f(xn+2h,yn+2hk2)=f(xn+h,yn+hk3)=yn+6h(k1+2k2+2k3+k4)
3.8 可视化结果
figure;
hold on;
plot(x_exact, y_exact_vals, 'k-', 'LineWidth', 2, 'DisplayName', '解析解');
plot(x_euler, y_euler, 'b--o', 'LineWidth', 1.5, 'MarkerSize', 4, 'DisplayName', '欧拉法');
plot(x_rk4, y_rk4, 'r-+', 'LineWidth', 1.5, 'MarkerSize', 4, 'DisplayName', '四阶RK4法');
title('欧拉法与四阶龙格库塔法(RK4)对比');
xlabel('x');
ylabel('y');
legend('Location', 'best');
grid on;
hold off;
- 绘制解析解、欧拉法解和RK4法解的曲线,便于直观对比。
3.9 误差分析
% 计算欧拉法和RK4与解析解的误差
error_euler = abs(y_euler - y_exact_vals);
error_rk4 = abs(y_rk4 - y_exact_vals);
% 绘制误差曲线
figure;
hold on;
plot(x_exact, error_euler, 'b--', 'LineWidth', 1.5, 'DisplayName', '欧拉法误差');
plot(x_exact, error_rk4, 'r-', 'LineWidth', 1.5, 'DisplayName', '四阶RK4法误差');
title('欧拉法与四阶龙格库塔法(RK4)误差对比');
xlabel('x');
ylabel('误差');
legend('Location', 'best');
grid on;
hold off;
- 计算并绘制欧拉法和RK4法的误差曲线,展示两种方法的精度差异。
4. 运行结果与分析
4.1 解的对比
- 解析解:蓝色实线。
- 欧拉法解:蓝色虚线带圆圈标记。
- RK4法解:红色实线带加号标记。
从图中可以看到:
- 解析解保持 y = 1 + e − 2 x y = 1 + e^{-2x} y=1+e−2x 的形状。
- 欧拉法在步长较大时,误差逐渐积累,数值解偏离解析解。
- RK4法的数值解与解析解几乎重合,显示出高精度。
4.2 误差分析
- 欧拉法误差:随着 x x x 增大,误差逐渐增大。
- RK4法误差:误差保持在极小的范围内,几乎可以忽略。
结论:
- 欧拉法是一种简单但精度较低的数值方法,适用于对精度要求不高的简单问题。
- **四阶龙格库塔法(RK4)**具有更高的精度和更好的稳定性,适用于需要高精度数值解的复杂问题。
5. 优化建议: 使用MATLAB内置ODE求解器
- MATLAB提供了强大的内置ODE求解器(如
ode45
),基于自适应步长的Runge-Kutta方法,可以自动调整步长,提高求解效率和精度。
6. 总结
通过本例,我们实现并比较了欧拉法和四阶龙格库塔法(RK4)在求解常微分方程中的表现。结果显示,RK4法在相同步长下具有显著更高的精度和更好的稳定性,适用于更复杂和精度要求更高的问题。而欧拉法虽然简单易实现,但在精度和稳定性上存在明显不足。
在实际应用中,根据具体问题的需求选择合适的数值方法至关重要。对于需要高精度的科学计算和工程应用,推荐使用RK4或更高阶的Runge-Kutta方法;对于简单或教育用途,欧拉法仍然是一个不错的选择。
订阅与关注
如果你对优化算法、数据科学感兴趣,欢迎订阅我的博客,获取最新的技术文章和研究动态。也可以通过以下方式与我交流:
关于作者
**[Cherngul]**是一位热爱算法与数据科学的技术达人,专注于探索各种优化算法在实际中的应用,希望通过分享,让更多人爱上算法的奇妙世界!
标签
- 四阶龙格库塔法(RK4)
- 欧拉法
- MATLAB
- 数值分析
最后
非常感谢您抽出时间来阅读我的文章!您的意见非常宝贵。文中可能有些地方表达得不够准确或错误,如果您遇到有需要改进或调整的地方。如果有任何问题或建议,请随时告诉我,我会非常感激。再次感谢您的阅读!
版权信息
© 2024 [Cherngul]. 保留所有权利。
开放原子开发者工作坊旨在鼓励更多人参与开源活动,与志同道合的开发者们相互交流开发经验、分享开发心得、获取前沿技术趋势。工作坊有多种形式的开发者活动,如meetup、训练营等,主打技术交流,干货满满,真诚地邀请各位开发者共同参与!
更多推荐
所有评论(0)