用MATLAB求解抛体运动二阶微分方程组时odeToVectorField报错如何解决

报错原因 #

• 核心错误：odeToVectorField要求输入的微分方程必须显式表示为最高阶导数等于右侧表达式的准线性形式，你的代码中仅定义了方程的右侧项，没有将其与二阶导数d2x、d2y做等式绑定，符号引擎无法识别这是二阶微分方程，因此抛出非准线性方程的报错。
• 潜在问题1：atan(dy/dx)的写法存在x方向速度趋近于0时的除零风险，且可以通过三角恒等式直接化简，降低计算复杂度。
• 潜在问题2：初始条件中保留了u、theta两个未赋值的符号变量，且ode数值求解器要求输入的初始条件为数值向量，不能直接传入符号条件数组。

修正方案 #

核心修改点 #

1. 将方程显式写为d2x == 右侧表达式、d2y == 右侧表达式的标准形式
2. 利用三角恒等式化简阻力项：
   • cos(atan(dy/dx)) * (dx² + dy²) = dx * sqrt(dx² + dy²)
   • sin(atan(dy/dx)) * (dx² + dy²) = dy * sqrt(dx² + dy²)
3. 给u、theta赋值具体数值，将初始条件转换为符合ode求解要求的数值状态向量

修正后代码 #

syms x(t) y(t) a b c d u theta
% 定义导数项
dx = diff(x,t);
dy = diff(y,t);
d2x = diff(x,t,2);
d2y = diff(y,t,2);

% 赋值常数（可根据实际需求调整）
a = 1;
b = 2;
c = 3;
d = 4;
u = 10; % 初始速度V0，自行调整数值
theta = pi/6; % 初始发射角30度，自行调整数值

% 显式定义微分方程
eq1 = d2x == -a * dx * sqrt(dx^2 + dy^2);
eq2 = d2y == -b * (c + d * dy * sqrt(dx^2 + dy^2));

% 转换为向量场
V = odeToVectorField([eq1, eq2]);
M = matlabFunction(V, 'vars', {'t','Y'});

% 初始状态向量：[x(0), dx(0), y(0), dy(0)]，顺序与odeToVectorField输出对应
y0 = [0, u*cos(theta), 0, u*sin(theta)];
interval = [0 5];  % 求解时间区间
ySol = ode23(M, interval, y0);

% 可选：运行完成后可调用下方代码绘制结果
% fplot(@(t)deval(ySol,t,1), interval, 'DisplayName','x(t)')
% hold on
% fplot(@(t)deval(ySol,t,3), interval, 'DisplayName','y(t)')
% legend

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