我正在尝试模拟具有以下方程描述的双线性刚度的系统：

$25\ddot{x}+15\dot{x}+330000x = p(t)$如果$x < 0.00072$

$25\ddot{x}+15\dot{x}+930000x = p(t)$如果$x >= 0.00072$

$p(t)$是我在另一个问题How to solve an ode with stochastic time-dependent input中已经讨论过的白噪声。过渡点称为$d$在代码中。

首先我定义输入和方程：

% design FIR filter to filter noise in the range [10-25] Hz
Fs = 1000; % sampling frequency

function [y] = nb(t)
    b = fir1(64, [10/(Fs/2) 25/(Fs/2)]);
    % generate Gaussian (normally-distributed) white noise
    noise = randn(length(t), 1);
    % apply filter to yield bandlimited noise
    y = filter(b,1,noise);
end

% State space functions ===================================================
odefun1 = @(t,u,m,k1,c,nb)[u(2); 1/m*(-k1*u(1)-c*u(2)+nb(t))];
odefun2 = @(t,u,m,k2,c,nb)[u(2); 1/m*(-k2*u(1)-(k1-k2)*d-c*u(2)+nb(t))];

然后我选择起始方程

% Time infos
T = 1/Fs;
tspan = 0:T:7.7;
tstart = tspan(1);
tend = tspan(end);

% Initial conditions
q0 = [0;0];

% Choose starting equation
if q0(1,1) > d
    flag1 = 1;
else
    flag1 = 0;
end

并启动集成循环：

% Options for ODE solver
options = odeset('RelTol',1e-6,'AbsTol',1e-8,'Events',@events);

% Accumulators for ODE solve output
tout = tstart;  % global solution time
qout = q0.';    % global solution position
teout = [];     % time at which events occur
qeout = [];     % position at which events occur
ieout = [];     % flags which event trigger the switch

% Main loop
while tout(end) < tend
    if flag1 == 1
        [t,q,te,qe,ie] = ode45(@(t,u)odefun2(t,u,m,k2,c,@nb),tspan,q0(:,1),options);
        flag1 = 0;
    else
        [t,q,te,qe,ie] = ode45(@(t,u)odefun1(t,u,m,k1,c,@nb),tspan,q0(:,1),options);
        flag1 = 1;
    end
% accumulate data into accumulators vectors
nt = length(t);
tout = [tout; t(2:nt)];
qout = [qout; q(2:nt,:)];
teout = [teout; te]; % Events at tstart are never reported.
qeout = [qeout; qe];
ieout = [ieout; ie];
tspan = [t(end), tend]; % reset time span,'where you just ended to tend'
q0 = [q(end,1);q(end,2)]; % reset IC as where you just ended
end

最后我定义了事件：

function [value,isterminal,direction] = events(t,q)
    value = [q(1)-d,q(1)-d]; % Detect zero of event functions
    isterminal = [1,1]; % Stop the integration
    direction = [-1,1]; % Direction
end 

这是正确的方法吗？因为当我以这种方式绘制力与位移矢量时，我没有得到刚度变化：

q2dot = zeros(size(qout(:,2))); % integrate veloocity to get acceleration
q2dot(1) = 0;
for i = 2:length(qout(:,2))
    h = tout(i)-tout(i-1);
    q2dot(i) = (qout(i,2)-qout(i-1,2))/h; 
end

RF = nb(t)-m*q2dot; % restoring force

编辑

我开始认为问题在于我定义白噪声并在方程中使用它的方式，因为如果我使用形式的强制函数$A \sin(\omega t)$（A足够高）我得到了斜率的变化