用 Euler 法、Heun 法和改进的 Euler 法求解常微分方程

[复制链接]

字体大小: 正常放大

2744557306

1189 主题	4 听众	2934 积分

该用户从未签到

电梯直达

1^#

发表于 2023-12-31 15:58 |只看该作者 |倒序浏览

|招呼Ta 关注Ta

这是一个 MATLAB 函数，用 Euler 法、Heun 法和改进的 Euler 法求解常微分方程，并比较它们的结果。以下是对你的代码的主要部分的解释：
function y = Euler1(a, b, N, af)
h = (b - a) / N;
x(1) = a;
y(1) = af;
yg(1) = af;
yh(1) = af;
jqj(1) = af;

% 迭代计算
for i = 2:N+1
      % Euler 法
      y(i) = y(i-1) + h * f(x(i-1), y(i-1));

      % Heun 法
      yh(i) = yh(i-1) + (h/4) * (f(x(i-1), yh(i-1)) + 3 * f(x(i-1) + 2 * h/3, yh(i-1) + 2 * h * f(x(i-1), yh(i-1)) / 3));

      % 改进 Euler 法
      yg(i) = yg(i-1) + h * (f(x(i-1), y(i-1)) + f(x(i), y(i) + h * f(x(i-1), y(i-1)))) / 2;

      x(i) = a + (i-1) * h;
      jqj(i) = x(i) + exp(-x(i));
end

% 计算误差
er = sum((y - jqj).^2);    % Euler 法误差
erg = sum((yg - jqj).^2); % 改进 Euler 法误差
erh = sum((yh - jqj).^2); % Heun 法误差

% 输出结果和误差
[x', y', yg', yh', jqj']
disp(['Euler法误差: ', num2str(er)]);
disp(['改进Euler法误差: ', num2str(erg)]);
disp(['Heun法误差: ', num2str(erh)]);

% 绘制结果
plot(x, y, 'r', x, yg, 'b', x, yh, 'k', x, jqj, 'g');
legend('Euler法', '改进Euler法', 'Heun法', '精确解');
end

% 待解的常微分方程的右侧函数
function result = f(x, y)
result = -y - exp(-x);
end

这个函数接受四个参数：

1.a 和 b：求解ODE的时间范围。
2.N：迭代步数。
3.af：初始条件。

在函数内部，它使用 Euler 法、Heun 法和改进的 Euler 法分别计算解，并计算了它们与精确解的误差。最后，它输出结果矩阵、误差，并通过 plot 函数将结果绘制出来。
你可以通过调用这个函数并提供合适的参数来运行它，例如：
Euler1(0, 1, 10, 1);

这将演示在给定的时间范围内使用三种不同的数值方法求解简单的常微分方程，并比较它们的结果。

zan