Switch between 2 equations to be given to ODE using events

Question

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

Answer 1

You're on the right track are correct to be using events for this purpose, but you're doing several things wrong that will result in an incorrect output.

1. You alway use the same initial conditions. I don't know your system, but it looks like you're trying to continuously integrate over time. When you switch from one set of ODEs to the other you need to specify the new initial conditions as the last state of of the previous equations.

2. You're doing something very similar with the tspan vector. This is particularly important as it looks like your ODEs are time dependent.

3. You're overwriting the values of t, q, etc. output by ode45 when you switch. Or is that not shown (% accumulate data into accumulators vectors comment)? It hard to help based on incomplete non-runnable code.

4. Based on the large stiffness values in your second order equations, you may want to look into using a stiff solver. ode15s is usually a good starting point to compare to what you got with a general solver like ode45. Cleve Moler recently wrote an excellent blog post detailing aspects of stiffness and the methods available in Matlab.

5. If you care about detailed statistical properties of your results, you probably shouldn't be using any adaptive scheme or solver designed for ODEs if you're you're working with noise. Using such methods may change the color of your noise and/or the resultant probability distributions of your outputs. You have have a stochastic differential equation (SDE). I refer you to this paper, which warns agains using ODE methods for SDEs. There are many cases where one can get away with what you're doing, but others that will lead incorrect results.

Though the example doesn't use events (the switching is based on time rather than state), see my answer here on StackOverflow.com/Matlab for how you might do 1–3 above.