How to solve nonlinear second order ODE in Matlab? [duplicate]

I am working on simulating a car suspension system using Matlab.

Specifically, I have to derive equation of motion using the Lagrange method and then use ode 45 to solve it. However, while using odetoVectorField I keep getting the error: System contains a nonlinear equation in diff(theta_rod_1(t), t). The system must be quasi-linear: highest derivatives must enter the differential equations linearly.

Is there anyway to get around this?

Here is part of the code:

%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);

• Would Computational Science be a better home for this question? Jul 10 at 13:09
• This program was a simulation of a mechanical system so I thought I should put it here, but you are right
– Luận Trần Công
Jul 10 at 14:59
• I’m voting to close this question because it is a question about programming a simulation and not a physics concept. Jul 11 at 2:47