Chattering effect using ode23s -- discontinuous state-space matrix

Question

I am attempting to model the following SDOF system with a variable spring and having a sinusoidal input, $$m\ddot{x}+(k_0+k_1)x = m(2\pi f)^2\sin(2\pi ft)$$ where $m$ is the mass, $k_0$ the original spring and $k_1 = k_0$ for $x\dot{x}>0$ and $k_1 = 0$ otherwise.

I was able to code this in Matlab, but for low frequencies, say $f=0.05$, a chattering effect would happen in the plot of $x\dot{x}$ vs. $t$. The following is my code and plots of the result for low frequencies.

f = 0.05;
t_f = 10;
tspan = [0 t_f];
z0 = [0 ; 0];

options = odeset( 'RelTol', 1e-13, 'AbsTol', 1e-13, 'Stats','on' );


[t, z] = ode23s(@(t,z) sys_13(t, z, f), tspan, z0, options);

figure(1)
set(gcf,'units','normalized','outerposition',[0 0 1 1])

subplot(2,1,1);
plot(t,z(:,1),'b');
title(['Time history of $y(t)\,\,(\ddot{u}_{gx}=-(2\pi f)^2\sin(2\pi f t)) \,\, f=\, $' num2str(f) ' $ $'],'interpreter','latex','fontsize',20)
ylabel('$y$','interpreter','latex','fontsize',15)
xlabel('$t$','interpreter','latex','fontsize',15)
xlim(tspan)

grid on
hold on

subplot(2,1,2);
plot(t,z(:,2),'b')
ylabel('$\dot{y}$','interpreter','latex','fontsize',15)
xlabel('$t$','interpreter','latex','fontsize',15)
xlim(tspan)

grid on
hold on

figure(2)
plot(t,z(:,1).*z(:,2),'b')

grid on
grid minor
hold on

base_stiffness = 4;
scale = 1e-5;
k0 = base_stiffness*ones(size(z(:,1),1),1)*scale;
k1 = zeros(size(z(:,1),1),1);

for ii = 1:size(z(:,1),1)
    k1(ii) = ( 1./(1 + exp( -(z(ii,1)*z(ii,2))/(1e-14) )) )*base_stiffness*scale;
end

plot(t, k1)

function dz = sys_13(t, z , f)
dz = zeros(2,1);

m = 1;
k0 = 4;
% tol = 1e-6;
% if abs(z(1)*z(2)) > tol && z(1)*z(2) > 0
%     k1 = 4;
% else
%     k1 = 0;
% end

% if z(1)*z(2) > 0
%     k1 = 4;
% else
%     k1 = 0;
% end

k1 = ( 1/(1 + exp( -(z(1)*z(2))/(1e-14) )) );
finput = ( -(2*pi*f)^2 )*sin( 2*pi*f*t );

dz(1) = z(2);
dz(2) = - finput - ((k0 + k1)/m)*z(1);
end

The plot of the time history of the position and velocity are

The plot of $x\dot{x}$ vs. $t$ is

If I zoom in close to where the plot is zero I see the following behavior

The plots with reduced absolute error and relative error tolerances set to 1e-9 are

Time History plot:

The $x\dot{x}$ plot:

Zoomed in $x\dot{x}$ plot:

I also tried applying a tolerance so to avoid the constantly switching on and off and instead the chattering occurs about the tolerance I set.

Now, I'm not sure what is causing these issues, but my best guess would be that the discontinuity in the state-space matrix $$ A=\begin{bmatrix}0&&1\\-\frac{k_0+k_1}{m}&&0\end{bmatrix}$$ when $k_1$ switches from $0$ to $k_0$ or vice versa is causing some sort of numerical errors in the solver. I was thinking about using an event location to find when $x\dot{x}$ passes through zero and start/stop the solver to avoid the discontinuity.

Thank you for any help or advice!

EDIT: Using the advice from LutzL I was able to remove the chattering, yet the way the additional stiffness $k_1$ is added to the system is gradual, whereas it should be more step like.

The $x\dot{x}$ plot with the additional stiffness value is:

Answer 1

The "chattering" may not just be a numerical artefact: it may be telling you that the solution ceases to exist. Note that the Existence and Uniqueness Theorem just guarantees existence of a solution when the functions in the differential equation are continuous.

Suppose as $t \to t_0$ you approach a point $x = x_0 $, $\dot{x} = 0$, with $\dot{x} > 0$ for $t < t_0$. Suppose $$ \frac{k_0}{m} < -(2\pi f)^2 \sin(2\pi f) < \frac{2k_0}{m} $$ Since $x \dot{x} > 0$, $k_1 = k_0$ and $\ddot{x} = -(2\pi f)^2 \sin(2 \pi f) - 2 k_0/m$, which will be negative. But if you cross $\dot{x} = 0$, $x \dot{x} < 0$, so $k_1 = 0$ and $\ddot{x} = -(2\pi f)^2 \sin(2 \pi f) - k_0/m$, which all of a sudden is positive. The result is that $\dot{x}$ can't actually go from positive to negative here. But it can't stay positive or $0$ here either. Thus the solution ceases to exist.