如何在Matlab中用欧拉法和二阶龙格-库塔法求解一阶微分方程组

错误根源 #

1. 符号函数参数不匹配：用syms f(t)定义了仅接受t单输入的符号函数，但调用时传入了t(i)和y_num(:,i)两个参数，直接触发参数数量不匹配的报错。
2. 微分方程组未正确实现：没有把二阶微分方程转化为可调用的一阶方程组函数。
3. 精确解逻辑混乱：当前exact_sol的定义完全错误，未给出对应微分方程的解析表达式。
4. 数组维度不匹配：y_exact定义为1行结构，但实际需要存储两个变量（y和y'）的解，维度不匹配会导致后续计算报错。

修正后的完整代码 #

核心逻辑说明 #

原二阶微分方程 $y'' = g(t) - 2y' - 2y$，转化为一阶方程组：
令 $y_1 = y$，$y_2 = y'$，则：
$$
\begin{cases}
y_1' = y_2 \
y_2' = g(t) - 2y_2 - 2y_1
\end{cases}
$$
这里默认$g(t)=0$，若你有具体的$g(t)$形式，直接替换即可。

clear;
clc;
clf;

% 定义一阶微分方程组：接受t和y（列向量），返回导数列向量
f = @(t, y) [y(2); -2*y(2) - 2*y(1)];

% 计算参数设置
h = pi/50;
t0 = 0;
tfinal = 2*pi;
t = t0:h:tfinal;
n = length(t) - 1;  % 直接用长度差避免round误差

% 初值：y(0)=1，y'(0)=0
y0 = [1; 0];

% -------------------------- 欧拉法求解 --------------------------
y_euler = zeros(2, n + 1);
y_euler(:, 1) = y0;
for i = 1:n
    y_euler(:, i+1) = y_euler(:, i) + h * f(t(i), y_euler(:, i));
end

% -------------------------- 二阶龙格-库塔法（Heun法）求解 --------------------------
y_rk2 = zeros(2, n + 1);
y_rk2(:, 1) = y0;
for i = 1:n
    k1 = f(t(i), y_rk2(:, i));
    k2 = f(t(i) + h, y_rk2(:, i) + h * k1);
    y_rk2(:, i+1) = y_rk2(:, i) + h * (k1 + k2)/2;
end

% -------------------------- 精确解计算 --------------------------
% 初值y(0)=1, y'(0)=0，g(t)=0时的解析解
exact_y1 = @(t) exp(-t).*(cos(t) + sin(t));
exact_y2 = @(t) -2*exp(-t).*sin(t);

y_exact = zeros(2, n + 1);
y_exact(1, :) = exact_y1(t);
y_exact(2, :) = exact_y2(t);

% -------------------------- 绘图对比 --------------------------
figure('Name','欧拉法与精确解对比');
plot(t, y_euler(1, :), '-k', t, y_exact(1, :), '--r');
xlabel('t'); ylabel('y(t)');
legend('欧拉法数值解','精确解');
title('欧拉法：数值解与精确解对比');
grid on;

figure('Name','二阶龙格-库塔法与精确解对比');
plot(t, y_rk2(1, :), '-b', t, y_exact(1, :), '--r');
xlabel('t'); ylabel('y(t)');
legend('二阶RK法数值解','精确解');
title('二阶龙格-库塔法：数值解与精确解对比');
grid on;

关键修正点 #

• 替换符号函数为匿名函数：f = @(t, y) [y(2); -2*y(2)-2*y(1)] 正确接受t和y两个输入，返回一阶方程组的导数列向量，彻底解决参数不匹配问题。
• 精确解修正：根据初值和微分方程推导了正确的解析表达式，数组维度设置为2行（对应y和y'），匹配数值解结构。
• 二阶龙格-库塔法实现：采用Heun法（平均斜率法），通过计算两个斜率取平均更新数值解，精度远高于欧拉法。
• 维度统一：所有数值解和精确解数组均定义为2行n+1列，避免索引和计算时的维度错误。

内容的提问来源于stack exchange，提问作者candyms801