SDOF spring mass damper system with a nonlinear spring

I am solving a linear spring mass damper system with the following equation:

The code I used was this:

clear all
clc
% damped resonant
% spring mass system
% 
y0=[0;0]; % [init_vel: init_disp]
tspan=[0 40];
k=9; 
m=4;
F=10;
c=0;
w=sqrt(k/m)/2;

[Td,Yd]=ode45(@(t,y) damped_spring_mass(F,m,k,w,c,t,y),tspan,y0);
plot(Td,Yd(:,1));

function dy=damped_spring_mass(F,m,k,w,c,t,y)
x=y(1); 
v=y(2);
dy=zeros(2,1);
dy(1)=v;
dy(2)=(F*cos(w*t)-k*x-c*v)/m;
end

The code was downloaded from code, contributed by author. This worked fine until I tried to do something like this:

by changing the last line of function into this:

dy(2)=(F*cos(w*t)-k*x^2-c*v)/m;

It showed a warning and the graph was wrong

Warning: Failure at t=4.169703e+00.  Unable to meet integration tolerances without reducing the step size below the smallest value allowed (1.421085e-14) at
time t. 
> In ode45 (line 360)
In code1 (line 14) 

I also tried this

But this one gave a good result.

Whenever exponent of x is even, I was getting warnings and graphs were wrong.

Can someone help me with the problem here? Thank you.

Answer 1

[Sam Chak]

Hi @jithin mathew

You've got the error message simply because the system is unstable. In other words, the nonlinearity causes the solution to blow up. For simplicity, let's consider that there is no external force, and just analyze the nonlinearities of the mass–spring systems according to Lyapunov stability theorem.

Nonlinear system Type-1:

, for and .

The result implies that the nonlinear system is unstable, and the states will diverge when both and (confirmed by your plot of x above and my results below.

Nonlinear system Type-2:

, for and .

The result implies that the nonlinear system is marginally stable.

Script for simulating the nonlinear mass–spring system Type-1:

k = 9; 
m = 4;
F = 10;
c = 0;
w = sqrt(k/m)/2;

tspan = linspace(0, 3, 1001)';
fun = @(t, x) [x(2); ...
             - k*x(1)^2 - c*x(2) + F*cos(w*t)];
[t, x] = ode45(fun, tspan, [1, 0]);
plot(t, x)

Result: