如何在MATLAB中绘制二阶微分方程解的图像及代码修复

原问题 #

需要绘制以下二阶微分方程的解曲线：

φ(t) + (d²φ(t)/dt²) - ((d²φ(t)/dt²) * m * sin(π) * a) = m * g * sin(π) * a

已知参数：

• a = v·t，v=1
• m=0.1
• t范围：0 < t < 1/v（即0<t<1）
• g=9.81

原代码问题分析 #

第一段符号解法代码问题 #

1. 先定义符号变量t，随后用linspace重新赋值t，覆盖了符号变量，导致ODE构建失败
2. 绘图时误用变量名ySol，实际解是phiSol
3. 关键问题：sin(π)=0，代入后方程直接简化为φ(t)=0，解为常数，失去绘图意义，推测是角度输入错误（比如应为sin(π/2)或其他非零角度）

第二段数值解法代码问题 #

1. ode45仅支持一阶微分方程组，不能直接使用二阶导数diff(phi,t,2)，必须先降阶为一阶方程组
2. 匿名函数中t是标量输入，但外部定义的t是数组，导致a = v*t无法匹配
3. 初始条件仅提供1个，二阶微分方程需要2个初始条件（φ(0)和φ’(0)）
4. 同样存在sin(π)=0导致方程无意义的问题

修正后的代码示例 #

假设将sin(π)替换为sin(π/2)（即1），以下给出两种可行解法：

方法1：符号解法 #

syms phi(t)
dphi = diff(phi);
d2phi = diff(phi, 2);
g = 9.81;
v = 1;
m = 0.1;
% 定义a为符号表达式，避免覆盖符号t
a_sym = v * t;
% 替换sin(pi)为sin(pi/2)，修正方程
ode = phi + d2phi - (d2phi * m * sin(pi/2) * a_sym) == m * g * sin(pi/2) * a_sym;
% 设定初始条件
cond1 = phi(0) == 0;
cond2 = dphi(0) == 1;
phiSol = dsolve(ode, cond1, cond2);
% 转换为数值函数并绘图
fphi = matlabFunction(phiSol);
t_vals = linspace(0, 1/v, 1000);
phi_vals = fphi(t_vals);
figure
plot(t_vals, phi_vals)
xlabel('t')
ylabel('φ(t)')
title('二阶微分方程符号解曲线')

方法2：数值解法（降阶为一阶方程组） #

二阶方程降阶：令y1 = φ(t)，y2 = φ’(t)，则：

• y1’ = y2
• y2’ = [mgsin(π/2)vt - y1] / [1 - m*sin(π/2)vt]

g = 9.81;
v = 1;
m = 0.1;
tspan = [0, 1/v];
% 初始条件：y1(0)=0, y2(0)=1
y0 = [0; 1];
% 定义一阶微分方程组
odefun = @(t, y) [y(2); (m*g*sin(pi/2)*v*t - y(1))/(1 - m*sin(pi/2)*v*t)];
% 求解
[t_vals, y_vals] = ode45(odefun, tspan, y0);
% 绘图
figure
plot(t_vals, y_vals(:,1))
xlabel('t')
ylabel('φ(t)')
title('二阶微分方程数值解曲线')

注意事项 #

如果sin(π)是正确输入，那么方程简化为φ(t)=0，解为常数0，绘图会是一条水平线，此时只需直接绘制plot(t_vals, zeros(size(t_vals)))即可。

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