I am trying to simulate a system with bilinear stiffness described by the following equations:
$25\ddot{x}+15\dot{x}+330000x = p(t)$ if $x < 0.00072$
$25\ddot{x}+15\dot{x}+930000x = p(t)$ if $x >= 0.00072$
$p(t)$ is the white noise I already talked about in another question How to solve an ode with stochastic time-dependent input . The transition point is called $d$ in the code.
First I define the input and the equations:
% 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))];
Then I select the starting equation
% 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
And start the loop for integration:
% 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
And finally I define the events:
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
Is this the right way to do it? Because I do not get a stiffness change when I plot the force vs displacement vectors in this way:
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
EDIT
I start thinking the problem is in the way I defined the white noise and I used it in the equations, because if I use a forcing function of the form $A \sin(\omega t)$ (with A high enough) I get the change in slope