Problem in solving four system of differential equation using dsolve command. I'm getting result dependent on other parameters

I am making half car model with tire striffness. In that case I've to solve 4 system of differential equation. I used dsolve function but I'm getting the results dependent of other parameter.

here I'm getting all the results in terms of 'z' . even z is showing itself in terms of z please help in finding out what I've done wrong.

here is the code which I've written :

%% 1.) -Parameter definition
% Masses and inertias
m                     = 1000;                   % Mass of the body [kg]
I                   = 2000;                   % Inertia of the body [kg*m^2]
m1 = 25;
m2= 25;

% Stiffness and damping values
k_f               = 60000;                  % Stiffness coefficient of spring [N/m]
c_f               = 5000;                   % Damping coefficient of damper [Ns/m]
k_r               = 60000;                  % Stiffness coefficient of spring [N/m]
c_r               = 5000;                   % Damping coefficient of damper [Ns/m]
k = 250;

% Lengths center of gravity to front and rear end
l_f                  = 2.5;                    % Distance of the right spring-damper to the center of mass [m]
l_r                  = 2.5;                    % Distance of the left spring-damper to the center of mass [m]

% Time and initial conditions
time                  = 0:0.005:1;                  % Time [s]
x_0                   = 0.1;                        % Initial Condition displacement [m]
x_dot_0               = 2;                          % Initial Condition velocity [m/s]
phi_0                 = 0;                          % Initial Condition angle [rad]
phi_dot_0             = 0;                          % Initial Condition angle velocity [rad/s]
y_0                   = 0.2;
y_dot_0               = 0;
z_0                   = 0.1;
z_dot_0               = 1;
    
%% 2.) Computing
%% 2.) Solving
    % define symbolic variables (dsolve)
syms x(t) y(t) z(t) phi(t)
Dx = diff(x,1);
D2x = diff(x,2);
Dy = diff(y,1);
D2y = diff(y,2);
Dz = diff(z,1);
D2z = diff(z,2);
Dphi = diff(phi,1);
D2phi = diff(phi,2);

% Solve equation
[phi,x,y,z] = dsolve([m*D2x + (c_r+c_f)*Dx + (c_f*l_f - c_r*l_r)*Dphi - c_r*Dz - c_f*Dy + (k_r+k_f)*x + (k_f*l_f - k_r*l_r)*phi - k_r*z - k_f*y ==0,...
    I*D2phi + (c_f*l_f - c_r*l_r)*Dx + (c_f*l_f*l_f + c_r*l_r*l_r)*Dphi + c_r*l_r*Dz - c_f*l_f*Dy + (k_f*l_f - k_r*l_r)*x + (k_f*l_f*l_f + k_r*l_r*l_r)*phi + k_r*l_r*z - k_f*l_r*y == 0,...
    m1*D2y + c_f*Dy - c_f*Dx - c_f*l_f*Dphi + (k_f+k)*y - k_f*x - k_f*l_r*phi == 0,...
    m2*D2z + c_r*Dz - c_r*Dx + c_r*l_r*Dphi + (k_r+k)*z - k_r*x + k_r*l_r*phi == 0],...
    [x(0)==x_0  ,Dx(0) == x_dot_0, y(0)==y_0  ,Dy(0) == y_dot_0, z(0)==z_0  ,Dz(0) == z_dot_0, phi(0) == phi_0, Dphi(0) == phi_dot_0],'t');

Answer 1

[Jason Louison]

Hello Ravi,

The "dsolve" command is only used for solving for the exact analytical solutions of linear and simple nonlinear ODE systems. It is not meant to handle complex nonlinear systems such as yours. The sensible way to solve this system would be numerical integration, and MATLAB has built in ODE solvers that can be used for this; the most commonly used being ode45. In order to do so, you must first convert the system of ODE's to a first order system, which can be done with the "odeToVectorField" command. Using this new array of ODE's, you then create a function hande for the ode45 command to use using the "matlabFunction" command. Then, specifiy the start and end time of the simulation, as well as your initial conditions in a vector for ode45 to use.

I have modified your code to accomodate this method of solving. See below:

clearvars
%% 1.) -Parameter definition
% Masses and inertias
m                     = 1000;                   % Mass of the body [kg]
I                   = 2000;                   % Inertia of the body [kg*m^2]
m1 = 25;
m2= 25;
% Stiffness and damping values
k_f               = 60000;                  % Stiffness coefficient of spring [N/m]
c_f               = 5000;                   % Damping coefficient of damper [Ns/m]
k_r               = 60000;                  % Stiffness coefficient of spring [N/m]
c_r               = 5000;                   % Damping coefficient of damper [Ns/m]
k = 250;
% Lengths center of gravity to front and rear end
l_f                  = 2.5;                    % Distance of the right spring-damper to the center of mass [m]
l_r                  = 2.5;                    % Distance of the left spring-damper to the center of mass [m]
% Time and initial conditions
t_start               = 0;                          % Simulation Start Time [s]
t_end                 = 15;                         % Simulation End Time [s]
x_0                   = 0.1;                        % Initial Condition displacement [m]
x_dot_0               = 2;                          % Initial Condition velocity [m/s]
phi_0                 = 0;                          % Initial Condition angle [rad]
phi_dot_0             = 0;                          % Initial Condition angle velocity [rad/s]
y_0                   = 0.2;
y_dot_0               = 0;
z_0                   = 0.1;
z_dot_0               = 1;
    
% define symbolic variables (dsolve)

syms x(t) y(t) z(t) phi(t)
Dx = diff(x,1);
D2x = diff(x,2);
Dy = diff(y,1);
D2y = diff(y,2);
Dz = diff(z,1);
D2z = diff(z,2);
Dphi = diff(phi,1);
D2phi = diff(phi,2);

%Ordinary Differential Equations 

eq_phi = m*D2x + (c_r+c_f)*Dx + (c_f*l_f - c_r*l_r)*Dphi - c_r*Dz - c_f*Dy + (k_r+k_f)*x + (k_f*l_f - k_r*l_r)*phi - k_r*z - k_f*y == 0;
eq_x = I*D2phi + (c_f*l_f - c_r*l_r)*Dx + (c_f*l_f*l_f + c_r*l_r*l_r)*Dphi + c_r*l_r*Dz - c_f*l_f*Dy + (k_f*l_f - k_r*l_r)*x + (k_f*l_f*l_f + k_r*l_r*l_r)*phi + k_r*l_r*z - k_f*l_r*y == 0;
eq_y = 1*D2y + c_f*Dy - c_f*Dx - c_f*l_f*Dphi + (k_f+k)*y - k_f*x - k_f*l_r*phi == 0;
eq_z = m2*D2z + c_r*Dz - c_r*Dx + c_r*l_r*Dphi + (k_r+k)*z - k_r*x + k_r*l_r*phi == 0;

eqs = [eq_phi eq_x eq_y eq_z]; %Ordinary Differential Equation Array

A = odeToVectorField(eqs); %Convert ODE's to first order systems

M = matlabFunction(A,'vars',{'t','Y'}); %Create Function Handle for ode45

initCond = [phi_0 phi_dot_0 x_0 x_dot_0 y_0 y_dot_0 z_0 z_dot_0]; %Initial Conditon Vector 

%ODE Solver: Function Handle, time Arrary, and Intial Condition Array:

sols = ode45(M,[t_start t_end],initCond); 

%To visualize the results in a smoother, more defined graph, we can use
%a time array of specified resolution

t = linspace(t_start,t_end,500); %define t as an array, with 500 points between t_start and t_end

y = deval(sols,t); %evaluate y at time t as an array function of ODE solver 

%Plotting All State Variables 

figure(1)

subplot(2,2,1)
plot(t,y(3,:),t,y(4,:))
legend('x','dx/dt')
title('Solution of State Variables (x)')
xlabel('Time (s)')
ylabel('Solutions (m or m/s)')

subplot(2,2,2)
plot(t,y(5,:),t,y(6,:))
legend('y','dy/dt')
title('Solution of State Variables (y)')
xlabel('Time (s)')
ylabel('Solutions (units)')

subplot(2,2,3)
plot(t,y(7,:),t,y(8,:))
legend('z','dz/dt')
title('Solution of State Variables (z)')
xlabel('Time (s)')
ylabel('Solutions (units)')

subplot(2,2,4)
plot(t,y(1,:),t,y(2,:))
legend('\phi','d\phi/dt')
title('Solution of State Variables (phi)')
xlabel('Time (s)')
ylabel('Solutions (rad or rad/s)')