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);

% 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;
    flag1 = 0;

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;
        [t,q,te,qe,ie] = ode45(@(t,u)odefun1(t,u,m,k1,c,@nb),tspan,q0(:,1),options);
        flag1 = 1;
% 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

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

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; 

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


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


1 Answer 1


You're on the right track and 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 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.

  • $\begingroup$ I added the whole code I am running $\endgroup$
    – Rhei
    Nov 21, 2014 at 7:46
  • $\begingroup$ By the way, the improvement in the results obatined using ode15s instead of ode45 is really astonishing $\endgroup$
    – Rhei
    Nov 23, 2014 at 9:13
  • $\begingroup$ You could add the point that the components of the state work on completely different scales. The ODE solvers of matlab are more geared towards producing nice looking plots for tutorial examples, so the least one can do is equalize the scales. This could be done by changing the time scale by a factor of 1000. Or set for the second state variable $1000u=25\dot x+15x$ so that $\dot x=40u-0.6x$, $\dot u=-330x+0.001p(t)$ etc. is a little better balanced. $\endgroup$ Dec 11, 2020 at 13:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.