0

I am trying to repeat an example I found in a paper.

I have to solve this ODE:

25 a + 15 v + 330000 x = p(t)

where p(t) is a white noise sequence band-limited into the 10-25 Hz range; a is the acceleration, v is the velocity and x the displacement.

Then it says the system is simulated using a Runge-Kutta procedure. The sampling frequency is set to 1000 Hz and a Gaussian white noise is added to the data in such a way that the noise contributes to 5% of the signal r.m.s value (how do I use this last info?).

The main problem is related to the band-limited white noise. I followed the instructions I found here https://dsp.stackexchange.com/questions/9654/how-to-generate-band-limited-gaussian-white-noise and I wrote the following code:

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

% Time infos (for white noise)
T = 1/Fs;
tspan = 0:T:4; % I saw from the plots in the paper that the simulation last 4 seconds
tstart = tspan(1);
tend = tspan (end);

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

Now I have to define the function for the ode, but I do not know how to give it the p (white noise)...I tried this way:

function [y] = p(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

odefun = @(t,u,m,k1,c,p)[u(2); 1/m*(-k1*u(1)-c*u(2)+p(t))];

[t,q,te,qe,ie] = ode45(@(t,u)odefun2(t,u,m,k2,c,@p),tspan,q0,options);

The fact is that the input excitation does not work properly: the natural frequency of the equation is around 14 Hz so I would expect to see the resonance in the response since the white noise is in the range 10-25 Hz.

I had a look at this Q/A also, but I can't still make it works:

How solve a system of ordinary differntial equation with time-dependent parameters

Solving an ODE when the function is given as discrete values -matlab-

Community
  • 1
  • 1
Rhei
  • 127
  • 11
  • Mathematically, there is no such thing as an "ODE with stochastic time dependent input." You have a stochastic differential equation (SDE). If you're interested in accurate measurement of statistical properties, you shouldn't use ODE methods to solve SDEs. – horchler Dec 13 '14 at 16:13
  • But since in the paper they used a Runge-Kutta procedure, does it mean there is one also for SDEs and not only for ODEs? – Rhei Dec 13 '14 at 17:24
  • What paper are you referring to? Sadly it's all too common for inappropriate methods to be used in published (and peer-reviewed) sources. Depending on the system, the noise, and what statistics one is interested in, you can sometimes "get away with it" in that the results are similar as when using a proper SDE solver. "Runge-Kutta procedure" is also unclear as there is such a thing as a [stochastic Runge-Kutta](http://en.wikipedia.org/wiki/Runge–Kutta_method_(SDE)) for SDEs. Despite the name, there are big differences. – horchler Dec 13 '14 at 18:20
  • I think you may need to learn more about your problem. Textbooks rather than random publications found on the Internet may be a better bet in terms of judging methods to to use for your application. If you do need to learn about solving SDEs numerically, the [Euler-Maruyama method](http://en.wikipedia.org/wiki/Euler–Maruyama_method) is a good (and reliable) starting point. Unfortunately, many of these subjects are off-topic for StackOverflow. – horchler Dec 13 '14 at 18:24
  • 1
    Lastly, the [`awgn`](http://www.mathworks.com/help/comm/ref/awgn.html) function may interest you. Or see [this question](http://stackoverflow.com/questions/23690766/proper-way-to-add-noise-to-signal). – horchler Dec 13 '14 at 18:34
  • Thank you for the suggestions, I will have a look. The paper I was talking about is this one: http://orbi.ulg.ac.be/handle/2268/18347 . I am trying to do the same numerical simulation reported in the paper to validate my program for the identification of non-linear structures. The paper in not free, I got it from my university so I cannot share it – Rhei Dec 13 '14 at 19:43
  • For a stiff system like that, if you're going to try to use an ODE solver to replicate what they did in the paper, you should either improve the tolerances for `ode45`, try a stiff solver like `ode15s`, or implement a fixed step-size Runge-Kutta scheme yourself and use sufficiently-small steps (it's never good for replication when they don't even bother to report the specific solver used and step-size if it was a fixed method). – horchler Dec 13 '14 at 19:55
  • Can I not retrieve the step size knowing the sampling frequency? – Rhei Dec 14 '14 at 22:24

0 Answers0