I am getting different outputs when using lsim vs ode45 when solving second order coupled differential equations

I have a model of a spring-mass-damper system that I have already solved numerially using the ode45 function. However, I would like to add controls to the system and thus have made an "equivalent" state space model. I checked the algebra a bunch and cant figure out why I am getting different responses with the lsim function vs the ode45 method. I have attached my code here.

%Define the parameters and system matrices as before
global g c k L m M Io f tau wn Ig K wb q z alpha
g=1.62; %lunar gravity
tau=-0.005; %envelope function 1/time constant
wn=117; %natural frequency of regolith foundation
m=2187500; %regolith foundation mass
M=3000;   %telescope mass (used from paper provided)
Ig=1500;  %Moment of Inertia for telescope Used Ig=1/2*M*R^2 assumed radius of 1m
Io=(M*L^2)+Ig;
c=42485.625; %damping coefficient 
L=10;         %length of telescope 
k=29.96*(10^9); %spring constant of regolith
K=7*(10^10);           %spring constant of aluminum 
wb=482; %angular velocity of seismic input
q=M+m;
z=(M^2)*(L^2);
alpha=Io*(q)-z;
% Define the state-space matrices
% Define the state-space matrices
A = [0, 1, 0, 0; (-Io * k / alpha), (-Io * c / alpha), ((M * L * K - z * g) / alpha), 0; 0, 0, 0, 1; ((M * L * k) / alpha), ((M * L * c) / alpha), (((M * L * g *q) - (K*q)) / alpha), 0];

Error using /
Matrix dimensions must agree.

B = [0; (Io / alpha); 0; (-M * L / alpha)];
C = [1, 0, 0, 0];
C2= [0, 0, 1, 0];
D = [0];
% Define the time vector
t = linspace(0, 900, 1000000);  % Adjust the time vector as needed
% Define the input signal
u= ((k.*exp(tau.*t).*sin(wb.*t))+c.*(((tau.*exp(tau.*t).*sin(wb.*t))+(wb.*exp(tau.*t).*cos(wb.*t)))));
sys=ss(A,B,C,D)
sys2=ss(A,B,C2,D)
y=lsim(sys,u,t);
y2=lsim(sys2,u,t);
figure
subplot(2,1,1)
plot(t,y)
xlabel('Time (s)')
ylabel('X Displacement (m)')
grid on
subplot(2,1,2)
plot(t,y2)
xlabel('Time (s)')
ylabel('$\theta$ Displacement (rad)')
%
%% ODE METHOD
global g c k L m M Io f  tau wn Ig K wb q z

g=1.62; %lunar gravity
tau=-0.005; %envelope function 1/time constant
wn=117; %natural frequency of regolith foundation
m=2187500; %regolith foundation mass
M=3000;   %telescope mass (used from paper provided)
Ig=1500;  %Moment of Inertia for telescope Used Ig=1/2*M*R^2 assumed radius of 1m
Io=(M*L^2)+Ig;
c=42485.625; %damping coefficient 
L=10;         %length of telescope 
k=29.96*(10^9); %spring constant of regolith
K=7*(10^10);           %spring constant of aluminum 
wb=117; %angular velocity of seismic input
q=M+m;
z=(M^2)*L^2;

f= @(time) ((k*exp(tau*time)*sin(wb*time))+c*(((tau*exp(tau*time)*sin(wb*time))+(wb*exp(tau*time)*cos(wb*time)))));

%Initial Conditions
x0=0;
v0=0;
theta0=0;
omega0=0;
IC=[x0,theta0,v0,omega0];

%Time Span (seconds)
t0=0;
tf=700;
tspan=[t0,tf];

% Numerical Integration (Both Portions)

%Linear Portion
[timeL,state_valuesL] = ode45(@linear,tspan,IC);
xL = state_valuesL(:,1);
vL = state_valuesL(:,2);
thetaL = state_valuesL(:,3);
omegaL = state_valuesL(:,4);

% Plotting the Results
figure

subplot(2,1,1)
plot(timeL,xL), ylabel('Displacement (m)')
title('X Displacement vs. Time (Linear)')
grid on
subplot(2,1,2)
plot(timeL,thetaL),ylabel('Angular Dispalcement (rad)')
title('$\theta$ Angular Displacement vs. Time (Linear)')

grid on
%print -depsc 2DOF_482



% Linear Function
function SD_L = linear(T,S)
global g c k L m M Io f tau wn Ig K wb y yprime q z alpha

SD_L = [S(2);
    (-Io*k/alpha)*S(1)+(-Io*c/alpha)*S(2)+((M*L*K-z*g)/alpha)*S(3)+(Io/alpha)*f(T);
    S(4);
    (M*L*k/alpha)*S(1)+(M*L*c/alpha)*S(2)+((M*L*g*q-K*q)/alpha)*S(3)-(M*L/alpha)*f(T)];
    
end

Answer 1

[Sam Chak]

Hi @Ethan Nobre

The reason for the discrepancy between the 'lsim()' solution and the 'ode45()' solution is due to the difference in the seismic angular velocity. In 'lsim()', you use wb = 482, while 'ode45()' uses wb = 117.

%% lsim METHOD

% Define the parameters and system matrices as before 
g       = 1.62;         % lunar gravity
tau     = -0.005;       % envelope function 1/time constant
wn      = 117;          % natural frequency of regolith foundation
m       = 2187500;      % regolith foundation mass
M       = 3000;         % telescope mass (used from paper provided)
Ig      = 1500;         % Moment of Inertia for telescope Used Ig=1/2*M*R^2 assumed radius of 1m
c       = 42485.625;    % damping coefficient 
L       = 10;           % length of telescope 
k       = 29.96*(10^9); % spring constant of regolith
K       = 7*(10^10);    % spring constant of aluminum 
wb      = 482;          % angular velocity of seismic input

% derivations
q       = M + m;
z       = (M^2)*(L^2);
Io      = (M*L^2) + Ig;
alpha   = Io*(q) - z;

% Define the state-space matrices
a21 = - Io*k/alpha;
a22 = - Io*c/alpha;
a23 = (M*L*K - z*g)/alpha;
a24 = 0;
a41 = M*L*k/alpha;
a42 = M*L*c/alpha;
a43 = (M*L*g*q - K*q)/alpha;
a44 = 0;

A = [  0,   1,   0,   0; 
     a21, a22, a23, a24; 
       0,   0,   0,   1; 
     a41, a42, a43, a44];

B = [0; 
     Io/alpha; 
     0; 
    -M*L/alpha];

C = [1, 0, 0, 0;
     0, 0, 1, 0];

D = [0];

% Define the time vector
t = linspace(0, 0.25, 5000001);  % Adjust the time vector as needed

% Define the input signal (the same as f in ode45)
u = k*exp(tau*t).*sin(wb*t) + c*(tau*exp(tau*t).*sin(wb*t) + wb*exp(tau*t).*cos(wb*t));

% State-space
sys = ss(A, B, C, D);
y0  = [0 0 0 0];
y   = lsim(sys, u, t, y0);

% Plotting the Results
figure(1)
subplot(2,1,1)
plot(t, y(:,1)), ylabel('Displacement (m)')
title('lsim: Time response of Displacement, x(t)')
grid on
subplot(2,1,2)
plot(t, y(:,2)), ylabel('Angular Dispalcement (rad)')
title('lsim: Time response of Angular displacement, \theta(t)')
grid on

%% ode45 METHOD

% Time Span (seconds)
t0      = 0;
tf      = 0.25;
tspan   = linspace(t0, tf, 5000001);

% Initial Condition
x0      = 0;
v0      = 0;
theta0  = 0;
omega0  = 0;
IC      = [x0, theta0, v0, omega0];

% Numerical Integration (Both Portions)

%Linear Portion
options = odeset('RelTol', 1e-8, 'AbsTol', 1e-10);
[t, y]  = ode45(@odefcn, tspan, IC, options);
xL      = y(:,1);
vL      = y(:,2);
thetaL  = y(:,3);
omegaL  = y(:,4);

% Plotting the Results
figure(2)
subplot(2,1,1)
plot(t, xL), ylabel('Displacement (m)')
title('ode45: Time response of Displacement, x(t)')
grid on
subplot(2,1,2)
plot(t, thetaL), ylabel('Angular Dispalcement (rad)')
title('ode45: Time response of Angular displacement, \theta(t)')
grid on

% Linear Function
function dydt = odefcn(t, y)
    % parameters
    g       = 1.62;         % lunar gravity
    tau     = -0.005;       % envelope function 1/time constant
    wn      = 117;          % natural frequency of regolith foundation
    m       = 2187500;      % regolith foundation mass
    M       = 3000;         % telescope mass (used from paper provided)
    Ig      = 1500;         % Moment of Inertia for telescope Used Ig=1/2*M*R^2 assumed radius of 1m
    c       = 42485.625;    % damping coefficient 
    L       = 10;           % length of telescope 
    k       = 29.96*(10^9); % spring constant of regolith
    K       = 7*(10^10);    % spring constant of aluminum 
    wb      = 482;          % angular velocity of seismic input

    % derivations
    q       = M + m;
    z       = (M^2)*L^2;
    Io      = (M*L^2) + Ig;
    alpha   = Io*(q) - z;
    f       = k*exp(tau*t)*sin(wb*t) + c*(tau*exp(tau*t)*sin(wb*t) + wb*exp(tau*t)*cos(wb*t));
    
    % Equations of motion
    dydt    = zeros(4, 1);
    dydt(1) = y(2);
    dydt(2) = - Io*k/alpha*y(1) -  Io*c/alpha*y(2) +   (M*L*K - z*g)/alpha*y(3) +  Io/alpha*f;
    dydt(3) = y(4);
    dydt(4) =  M*L*k/alpha*y(1) + M*L*c/alpha*y(2) + (M*L*g*q - K*q)/alpha*y(3) - M*L/alpha*f;
end

[Sam Chak]

Yes @Ethan Nobre, you're absolutely right. The system has a very large stiffness term, and the external input signal is also highly oscillatory. If the tolerances are not properly set, the ode45 solver will automatically determine them by default.

If you find the above solution helpful, please consider clicking 'Accept' ✔ on the answer and voting 👍 for it. Thanks a bunch!

Also check out the common issues that you might encounter when solving ODEs:

https://www.mathworks.com/help/matlab/math/troubleshoot-common-ode-problems.html