MATLAB求解含阶跃函数的ODE：结果异常无法绘图的解决咨询

DevOps 小满

2026-8-6

解决含阶跃函数常微分方程的MATLAB绘图问题

问题根源

原代码通过symsum生成无穷项阶跃函数序列，导致拉普拉斯变换后的表达式包含无穷级数，ilaplace无法生成可直接用于绘图的闭式解，最终结果无法计算数值进行可视化。

方案1：数值求解法（推荐，高效稳定）

直接使用MATLAB的ode45求解器，将分段激励定义为数值函数，无需复杂符号运算：

% 定义二阶ODE：d2y/dt2 + 0.1*dy/dt + y = f(t)
odefun = @(t,y) [y(2); -0.1*y(2) - y(1) + f(t)];

% 定义分段方波激励f(t)，周期2π，t≥0时交替取1和-1
function val = f(t)
    period = 2*pi;
    phase = mod(t, period);
    val = phase < pi ? 1 : -1;
end

% 设置求解区间（取10个周期足够观察稳态）和初始条件
tspan = [0, 10*pi];
y0 = [0; 0]; % y(0)=0, y'(0)=0

% 求解ODE
[t, y] = ode45(odefun, tspan, y0);

% 绘图
figure;
plot(t, y(:,1), 'LineWidth', 1.5);
xlabel('t');
ylabel('y(t)');
title('ODE的数值解');
grid on;

方案2：改进符号法（有限项近似）

若坚持使用符号运算，将无穷级数截断为有限项（可根据精度需求调整项数），使ilaplace生成可计算的表达式：

syms t s Y y(t) k
Dy = diff(y,t);
D2y = diff(Dy,t);

% 截断无穷级数为前10项（N值越大精度越高，计算时间越长）
N = 10;
f = heaviside(t) + 2*symsum((-1)^k*heaviside(t - k*pi), k, 1, N);

% 拉普拉斯变换求解
eqn = D2y + 0.1*Dy + y == f;
leqn = laplace(eqn, t, s);
LT_Y = subs(leqn, laplace(y,t,s), Y);
LT_Y = subs(LT_Y, y(0), 0);
LT_Y = subs(LT_Y, subs(diff(y(t),t), t, 0), 0);

Y = solve(LT_Y, Y);
y_sym = ilaplace(Y, s, t);

% 转为数值函数并绘图
y_fun = matlabFunction(y_sym);
t_vals = linspace(0, 10*pi, 1000);
y_vals = y_fun(t_vals);

figure;
plot(t_vals, y_vals, 'LineWidth', 1.5);
xlabel('t');
ylabel('y(t)');
title('ODE的符号近似解');
grid on;