我正在使用 Matlab 模拟汽车悬架系统。
具体来说，我必须使用拉格朗日方法推导出运动方程，然后使用 ode 45 来求解它。但是，在使用时，odetoVectorField我不断收到错误：系统包含一个非线性方程diff(theta_rod_1(t), t)。系统必须是准线性的：最高导数必须线性地进入微分方程。
有没有办法解决这个问题？
以下是部分代码：
%Defining relations between variables
th_M = asin( ( L_rod*cos(th_rod_1) - L_rod*cos(th_rod_2)) / (LMx - 2*(L1+L2) ) );
L_spring_1 = sqrt(L_rod^2 + L2^2 - 2*L_rod*L2*cos(pi/2 + th_M - th_rod_1) ) ; %%Law of cosine
th_spring_1 = pi/2 - th_M - acos((L2^2 - L_rod^2 + L_spring_1^2)/(2*L2*L_spring_1));
L_spring_2 = sqrt(L_rod^2 + L2^2 - 2*L_rod*L2*cos(pi/2 - th_M - th_rod_2) );
th_spring_2 = pi/2 + th_M - acos((L2^2 - L_rod^2 + L_spring_2^2)/(2*L2*L_spring_2));
dx = LMx/2 * th_M;
dth_M = (L_rod*sin(th_rod_1)*Dth_rod_1 - L_rod*sin(th_rod_2)*Dth_rod_2)/((1 - (L_rod*cos(th_rod_1) - L_rod*cos(th_rod_2))^2/(2*L1 + 2*L2 - LMx)^2)^(1/2)*(2*L1 + 2*L2 - LMx));
ddx = (LMx*(L_rod*sin(th_rod_1)*Dth_rod_1 - L_rod*sin(th_rod_2)*Dth_rod_2))/(2*(1 - (L_rod*cos(th_rod_1) - L_rod*cos(th_rod_2))^2/(2*L1 + 2*L2 - LMx)^2)^(1/2)*(2*L1 + 2*L2 - LMx));
%% Kinematics
% Define the displacements of the double spring in Cartesian coordinates
xm1 = (x0+L1+L2)*cos(th_M) - L_rod*sin(th_rod_1); %x-pos of the left wheel, the cos(th_M) projects the length onto the x-direction
ym1 = 0;
xm2 = (x0+LMx-L1-L2)*cos(th_M) + L_rod*sin(th_rod_2);
ym2 = 0;
xM = x0;
yM = L_rod*cos(th_rod_1) + LMy/2*cos(th_M) - (LMx/2 -L2-L1)*sin(th_M); %yM is defined as the height of mass M center of gravity
%Find the velocities by differentiating the displacements with respect to
%time (using diff function)
vxm1 = ((L_rod*cos(th_rod_1) - L_rod*cos(th_rod_2))*(L_rod*sin(th_rod_1)*Dth_rod_1 - L_rod*sin(th_rod_2)*Dth_rod_2)*(L1 + L2 + x0))/((1 - (L_rod*cos(th_rod_1) - L_rod*cos(th_rod_2))^2/(2*L1 + 2*L2 - LMx)^2)^(1/2)*(2*L1 + 2*L2 - LMx)^2) - L_rod*cos(th_rod_1)*Dth_rod_1
vym1 = 0;
vxm2 = L_rod*cos(th_rod_2)*Dth_rod_2 - ((L_rod*cos(th_rod_1) - L_rod*cos(th_rod_2))*(L_rod*sin(th_rod_1)*Dth_rod_1 - L_rod*sin(th_rod_2)*Dth_rod_2)*(L1 + L2 - LMx - x0))/((1 - (L_rod*cos(th_rod_1) - L_rod*cos(th_rod_2))^2/(2*L1 + 2*L2 - LMx)^2)^(1/2)*(2*L1 + 2*L2 - LMx)^2);
vym2 = 0
vxM=0;
vyM = ((L_rod*sin(th_rod_1)*Dth_rod_1 - L_rod*sin(th_rod_2)*Dth_rod_2)*(L1 + L2 - LMx/2))/(2*L1 + 2*L2 - LMx) - L_rod*sin(th_rod_1)*Dth_rod_1 + (LMy*(L_rod*cos(th_rod_1) - L_rod*cos(th_rod_2))*(L_rod*sin(th_rod_1)*Dth_rod_1 - L_rod*sin(th_rod_2)*Dth_rod_2))/(2*(1 - (L_rod*cos(th_rod_1) - L_rod*cos(th_rod_2))^2/(2*L1 + 2*L2 - LMx)^2)^(1/2)*(2*L1 + 2*L2 - LMx)^2);
%% Lagrange
% Kinetic energy
T = 0.5*M*ddx^2 + 0.5*J0*dth_M^2 + 0.5*m1*vxm1^2 + 0.5*m2*vxm2^2 ;
%
T = simplify(T);
% Potential energy
V = m1*g*ym1 + m2*g*ym2 + M*g*yM + 0.5*k1*(L_spring_1-L_spring_rest)^2 + 0.5*k2*(L_spring_2-L_spring_rest)^2 ;
V = simplify (V);
%Dissipation energy
R = 0.5*b*(vxm1*sin(th_spring_1))^2 + 0.5*b*(vxm2*sin(th_spring_2))^2;
Lagrange = T-V;
%Generalized coordinates:
q = [th_rod_1, th_rod_2];
dq = [Dth_rod_1, Dth_rod_2];
ddq = [DDth_rod_2, DDth_rod_2];
DL_Dq = jacobian(Lagrange,q')'; % Operator ' is transpose
DL_Ddq = jacobian(Lagrange,dq);
DDL_DtDdq = jacobian(DL_Ddq',[q, dq]) * [dq, ddq]';
DR_Ddq = jacobian(R,dq')'; %Damping term of Lagrange_Euler equation
EoM = DL_Dq - DDL_DtDdq ; %Euler-Lagrange equation with damping
%- DR_Ddq
EoM = EoM == 0;
%simplify(EoM(1));
%%
syms theta_rod_1(t) theta_rod_2(t)
EoM=subs(EoM,[th_rod_1, Dth_rod_1, DDth_rod_1, th_rod_2, Dth_rod_2, DDth_rod_2], ...
[theta_rod_1, diff(theta_rod_1), diff(theta_rod_1,2), theta_rod_2, diff(theta_rod_2), diff(theta_rod_2,2)]);
%% Compute with value
% m1 m2 M g k1 k2 b J0 x0 L1 L2 LMy LMx L_rod L_spring_rest
x0=0;
L1 = 1;
L2 = 3;
LMx = 10;
LMy = 5;
L_spring_rest = 2;
L_rod = 3;
m1 = 1;
m2 = 1;
M = 5;
g = 9.8;
k1 = 500;
k2 = 700;
b = 100;
J0 = 50;
eqn1 = subs(EoM(1))
eqn2 = subs(EoM(2));
%The two equations are nonlinear second-order differential equations.
%To solve these equations, convert them to first-order differential equations
[V,S] = odeToVectorField(eqn1,eqn2);