Solving for nonlinear second order ODE using ode45

Hi,

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