以下是使用欧拉公式、梯形公式和改进欧拉公式求解微分方程初值问题的Matlab程序，并与精确解进行比较：
```matlab
% 定义微分方程
f = @(t, y) t - y;
% 定义步长和时间区间
h = 0.1;
tspan = [0, 1];
% 定义初始值
y0 = 1;
% 精确解
t_exact = tspan(1):0.01:tspan(2);
y_exact = t_exact + exp(-t_exact);
% 欧拉公式
[t_euler, y_euler] = euler(f, tspan, y0, h);
euler_error = abs(y_euler - interp1(t_exact, y_exact, t_euler));
% 梯形公式
[t_trap, y_trap] = trap(f, tspan, y0, h);
trap_error = abs(y_trap - interp1(t_exact, y_exact, t_trap));
% 改进欧拉公式
[t_impeuler, y_impeuler] = impeuler(f, tspan, y0, h);
impeuler_error = abs(y_impeuler - interp1(t_exact, y_exact, t_impeuler));
% 绘图
plot(t_exact, y_exact, 'k-', t_euler, y_euler, 'r--', t_trap, y_trap, 'g-.', t_impeuler, y_impeuler, 'b:');
legend('Exact', 'Euler', 'Trapezoidal', 'Improved Euler');
xlabel('t');
ylabel('y');
title('Numerical Solutions of y''=t-y');
% 欧拉公式函数
function [t, y] = euler(f, tspan, y0, h)
t = tspan(1):h:tspan(2);
y = zeros(size(t));
y(1) = y0;
for i = 1:length(t)-1
y(i+1) = y(i) + h*f(t(i), y(i));
end
end
% 梯形公式函数
function [t, y] = trap(f, tspan, y0, h)
t = tspan(1):h:tspan(2);
y = zeros(size(t));
y(1) = y0;
for i = 1:length(t)-1
y(i+1) = y(i) + h/2*(f(t(i), y(i)) + f(t(i+1), y(i)+h*f(t(i), y(i))));
end
end
% 改进欧拉公式函数
function [t, y] = impeuler(f, tspan, y0, h)
t = tspan(1):h:tspan(2);
y = zeros(size(t));
y(1) = y0;
for i = 1:length(t)-1
y_temp = y(i) + h*f(t(i), y(i));
y(i+1) = y(i) + h/2*(f(t(i), y(i)) + f(t(i+1), y_temp));
end
end
```
903
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
