我正在尝试使用可变弹簧和正弦输入对以下单自由度系统进行建模，

$$m\ddot{x}+(k_0+k_1)x = m(2\pi f)^2\sin(2\pi ft)$$ 在哪里$m$是质量，$k_0$原来的弹簧和$k_1 = k_0$为了$x\dot{x}>0$和$k_1 = 0$否则。

我能够在 Matlab 中对此进行编码，但是对于低频，比如说$f=0.05$, 的情节中会发生颤抖效应$x\dot{x}$对比$t$. 以下是我的代码和低频结果图。

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

k0 = 4;
k1 = 4;
ktot = k0+k1
k = calcKSmooth(z(:,1).*z(:,2), k0, ktot )*1e-5;

plot(t,k)


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

m = 1;
k0 = 4;
k1 = 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
ktot = k0 + k1;
k = calcKSmooth(z(1)*z(2) , k0 , ktot);

finput = ( -(2*pi*f)^2 )*sin( 2*pi*f*t );

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

function k = calcKSmooth(xxdot, k0, ktot)
  k1 = k0; k2 = ktot; c = 0;
  r = 1e6;
  ecr = exp(c*r);
  erx = exp(r*xxdot);
  k = (k1*ecr + k2*erx)./(ecr+erx);
end

位置和速度的时间历程图是

的情节$x\dot{x}$对比$t$是

如果我放大靠近绘图为零的位置，我会看到以下行为

我还尝试应用一个容差，以避免不断地打开和关闭，而是发生关于我设置的容差的喋喋不休。

现在，我不确定是什么导致了这些问题，但我最好的猜测是状态空间矩阵中的不连续性

$$A=\begin{bmatrix}0&&1\\-\frac{k_0+k_1}{m}&&0\end{bmatrix}$$ 什么时候$k_1$从切换$0$到$k_0$反之亦然，导致求解器中出现某种数值错误。我正在考虑使用活动地点来查找时间$x\dot{x}$通过零并启动/停止求解器以避免不连续性。

感谢您的任何帮助或建议！

编辑：使用 Bill Greene 提供的答案，消除了颤振，但向系统添加刚度的方式是渐进的，但物理上的关闭或打开应该更突然，很像阶跃函数。

的情节$x\dot{x}$与$k=k_0+k_1$叠加的是：