What is the way of solving the noise added differential equations in matlab?

I want to study the stoichastic resonance in duffing oscillator. The duffing oscillator is given as,

,

where are constant parameter and f is the forcing amplitude and is the white gaussian noise and D is the varience. I want to see the effect of noise on the dyanmics of the system. I tried solving the above equation using ode45 but it's not giving any results......

clear all;
clc;
clf;
%% ----------------INPUT PARAMETERS---------------------------
=0.33; alpha=0.5; w0=sqrt(-1); beta=1; w=0.96; D=10;
x0 = [-1 0 0];                        %initial condition
tspan = [0:0.1:100];          %time duration
n_rec= 300;                        %recording time
options = odeset ('RelTol',1e-7,'AbsTol',1e-6);  %acuuracy set
[t,x] = ode45('duffing_sr',tspan,x0, options);     %solving using ode45
[row, col]=size(x);
%Storing the data after discarding initial transients
x1 = x(n_rec:row,1);
x2 = x(n_rec:row,2);
x3 = x(n_rec:row,3);
n = length(x1); % Length of x1 after discarding the initial transients
t_n=t(n_rec:row);
plot(t_n,x1);             %plotting
%---------------------------Function-----------------------------------------------
function dxdt = duffing_sr(t,x)
global  f alpha w0 beta w D ;
dxdt = zeros(3,1);
dxdt(1) =  x(2);
dxdt(2) = f*sin(x(3))-alpha*x(2)-w0^2*x(1)-beta*x(1)^3+sqrt(D)*(randn()-0.5);
dxdt(3) = w;
end

****************************************************************************************

I have no idea how to solve noise added differential equations.

Answer 1

[Alan Stevens]

Try to avoid the (mis)use of global!

%% ----------------INPUT PARAMETERS---------------------------
f=0.33; alpha=0.5; w0=sqrt(-1); beta=1; w=0.96; D=10;
x0 = [-1 0];                        %initial condition
tspan = 0:0.1:100;          %time duration
n_rec= 300;                        %recording time
%options = odeset ('RelTol',1e-7,'AbsTol',1e-6);  %acuuracy set
[t,x] = ode45(@(t,x) duffing_sr(t,x,f,alpha,beta,w,D),tspan,x0);     %solving using ode45
[row, col]=size(x);
%Storing the data after discarding initial transients
x1 = x(n_rec:row,1);
x2 = x(n_rec:row,2);

n = length(x1); % Length of x1 after discarding the initial transients
t_n=t(n_rec:row);
plot(t_n,x1);             %plotting
%---------------------------Function---------------------------------------
function dxdt = duffing_sr(t,x,f,alpha,beta,w,D)
     dxdt = [x(2);
             f*sin(w*t)-alpha*x(2)-(-1)*x(1)-beta*x(1)^3+sqrt(D)*(randn-0.5)];
end

[Paul]

Considering two steps of the ode45 solver,

first from t1 to t1 + dt1

second from t1 + dt1 to t1 + dt1 + dt2.

As I undersand it the solver calls duffing_sr at t + dt1 as part of the integration over the first step, and then calls duffing_sr at t + dt1 as part of the integration over the second step. On these calls, the state input, x, will be different. But is it important to ensure that the randn input is the same on those calls? IOTW, should the noise input should be the same at t + dt1 regardless of if it's being applied for the dxdt computation at the end of the first step or for the dxdt computation at the start of the second step?

Same question would apply if ode45 computes dxdt at the same times, like t1 + dt1/2, when computing deriatives within a step.

[Ms. Sita]

Thank you Alan Stevens for correcting the code.

[Jürgen]

I worry a bit about this approach, because ode45 accesses the derivative multiple times in one time step.

My intuition is, that each call of the derivative function generates a new random number at each sub-time step of the calculation and that will at least theoretically influence the outcome of the computation.

Somehow the algorithm has to ensure, that the random-generator is only activated once for each time-step.

Am I wrong?

[Torsten]

ode45 and the other ODE solvers are only suited to solve deterministic systems of ordinary differential equations. Any stochastic inputs during the solution process are unacceptable due to the numerical methods they use. These methods are usually based on n-fold differentiability of the derivatives function.

[haif S]

Hello, I also encountered the same problem, how did you solve it, is the ode solver unable to solve this problem

[Bruno Luong]

+1 to Torsen comment. Do not use ODE45 or friends. It assumes then the noise is constant over one interval of discretization which is not correct.

The correct way seems to be this https://www.mathworks.com/help/finance/sde.html

[Paul]

Simulink, which I think uses ode solvers similar to (if not exactly the same as) ODE45 and friends, suggests using an approach that approximates the noise as being constant over an interval that's small relative to the dynamics of the system, which will require the user having some knowledge of the system dynamics: Band-Limited White Noise.

Because it approximates a continuous-time signal a discrete-time signal at a fixed sampling period, this approach also mitigates the need to keep track of what the noise values are at minor time steps in the solver that may be hit more than once.

[Bruno Luong]

Not only the noise needs to be known but the local stiffness of the ode should have an impact.

If the noise is assumed to be constant, the discrete amplitude should be adapted to the noise pdf AND also the time step and the siffness. And possibly sde solver does that a certain degree, not ode45 and friends.

People just do this kind of discretization without any theoretical about how the system would behave.

That's how garbage parpers get published.

[Paul]

If by "local stiffness of the ode" you mean what I called "the system dynamics," then we are in agreement on this point.

If you look at the Algorithms section of the linked doc page, I think it discusses what you're referring to as adapting the discrete amplitude, which you'll see is dependent on the correlation time, which is the sample time.

The block doesn't adapt any parameter to the stiffness (at least what I think you mean by stiffness) of the system. It relies on the user to speficy a correlation time (time step) to be small enough "relative to the natural bandwidth of the system."

Not familar with the SDE solver you linked, so can't comment on that.

[Jürgen]

I think my statement is the same as that of Bruno Luong, Paul, Torsten... .

Many solvers use variable time-steps to cope for accuracy.

As a quick workaround I made a modified version of ode45 to provide a identical random noise for all the substeps that the solver uses to compute the current progress of the system.

My workaround version has a slightly different signature:

function varargout = ode45_NoiseJM(ode,diffInj,tspan,y0,options,varargin)

The original mathworks function has 2 loops:

% THE MAIN LOOP (~ line 270)
% LOOP FOR ADVANCING ONE STEP. (some lines later)
The first loop steps along the whole integration interval and the second loop copes for accuracy and thus uses autoadjusted time-steps (h).
% LOOP FOR ADVANCING ONE STEP.
nofailed = true;                      % no failed attempts
diffRates = diffInj(h); % generates a random injection for time interval h
% with function val = diffInj(h) looks like
% val = val_RMS_sqHz*sqrt(1/h);
while true
    y2 = y + h .* (b11.*f1 );
    t2 = t + h .* a2;
    f2 = ode(t2, y2)+diffRates;
    % and identical for the following f3, f4,.., f6 of the RK-scheme

I think this approach applies

• the same random injection for all RK-subsections
• and still allows the solver to use variable steps
• and scales the noise injection wiht the correct frequency-density (the *sqrt(1/h))

So I think wiht this approach I get random noise effects correctly in interaction wiht my system behavior and get approximately correct spectral shaping e.g..

I guess this is what is also done in more professional modules referred as stochastic ode mehtods.

I would appreciate some Mathworks Expert to comment and bring me back on the right track if I am completely wrong.

Also I am interested now in combining delayed ode with stochastic ode. But did not yet find such a minimum invasive way to modify dde23.

Here help is appreciated.

My area of work is currently the interaction of statistical processes with mechanical systems (MEMS) and I want to be as accurate as possible to correctly cope for interaction effects of statistical noise and system behavior and unfortunately also solver-properties interfere here.

[Bruno Luong]

The eventual issue I see when using ode45 is that ode45 adapts the time step based on an internal error estimator.

It was designed for smooth solution and non stiff ode. I have no idea how the random noise interact with the solver in this aspect.

I would not bet on the accuracy of the numerical solution and correctness of using ode45 for stochastic system. I don't know simulink and did not read Paul's document. Myself I never practice ode with stoschastic system. I hear people talk about noise "bandwidth" like we integrate of a linear system. My objection is that for non-linear system it is pretty hard to predict how the noise would behave. In non-linear system such as Navier Stokes there is a cascade of energy transmission between different scales, and no one quite fully understand yet this equation. I would imagine something like that could happen on a non-linear stiff ode too.

Anyway people still solve numerically Navier-Stokes equation eventhough they have no idea how good their numerical scheme are!

For those you are interested there is a submission oh FEX SDETools

[Paul]

Hi Bruno,

Moving more to the philosophical side ....

The output of any numeric solver of system of differential equations, or differential algebraic equations, or hybrid differential-discrete equations, or whatever, is always an approximation to the tue solution, and the equations themselves are often approximations to the underlying phyics of the problem (if representing a physical system).

Therefore, it's important for the modeler to have an understanding of the system to be able to have any confidence in the results of the simulation based on any any solve and always ask the question: "do the outputs of the solver make sense?" The less able the modeler is to answer that question, the higher the risk that the outputs of the solver will be accepted when they shouldn't be. But, if the system under simulation, even if noninear, is well understood, pehaps can be compared to test data for certain cases, etc., then simulating such systems with random process inputs may be feasible.

I don't know about Navier-Stokes, so can't comment on whether those equations can be reliably simulated with random process inputs, but I suspect that there are systems for which reasonable, approximate solutions can be obtained.

[Bruno Luong]

@Paul The "undestanding" a given equation in math is essentially proving existence, uniqueness of the solution, then proving the stability, then derive a numerical scheme and prove the convergence of the scheme and eventually find a bound of the error (wrt to discretization).

This is essentially basic numerical analysis of applied math person does.

Of course if nice- noise is introduced to a linear system with diffusion term, everything will work out.

Back to SDE: Using ode45 a classical Runge-Kutta scheme for stochastic system is far from understanding the tool. ODE45 is for classical ODE only.

The sde is entirely another beast. Here is how it works Runge-Kutta SDE, also take a look at the references on the bottom, and this paper. You'll see that the stochastic process must be rooted inside the solver, and it has little to do with MATLAB ode45.

May be it is not clear for few of you to insist but it is simply wrong to apply pseudo-noise on top of classical runge kutta scheme, since there is no knowledge of which prediction/correction and time-step. There are none theoretical justification of using ode45 for sde.

[Torsten]

@Jürgen

As far as I understand the numerical solution of SDEs, the numerical scheme to use is in general a fixed-step scheme and strongly depends on the SDE you are trying to solve:

https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_method_(SDE)

So I don't think that somehow randomly manipulating the Runge-Kutta steps will result in something useful and interpretable.

[Jürgen]

@Torsten

I 100% agree to what you said in your very first comment.

https://de.mathworks.com/matlabcentral/answers/1646055-what-is-the-way-of-solving-the-noise-added-differential-equations-in-matlab#comment_2695379

Especially with '...stochastic inputs during the solution process are unacceptable...'. This is exactely what I told my friends about that situation. One can see it even wiht tic-toc that the solver runs in chaos and takes much longer time.

I thought, with the workaround I found a way to place then random injection at a position where it does not corrupt the process. Do you agree with this?

For the extension to SDE I also see that they work and stay understandable if one restricts the algorithm to fixed step. I found no comparable simple way and have only a intuition that it might also work for variable step: But then my internal picture of it tells me that it related fixed or even variable past times with variable current time... and inevitably results in some kind of complicated implicit solution-process.

At the moment it somehow see it as a obscured picture it, but it is very confusing, what might be necessary to code for this.

[Bruno Luong]

• Try to use SDETools
• Use MATLAB SDE in Financial Toolbox
• Implement the Euler scheme in the Runge Kutta SDE wiki page to have a feeling
• Find a paper on the subject and implement it (you have few cited along the discussion)
• But NOT using ODE45

[Torsten]

I thought, with the workaround I found a way to place then random injection at a position where it does not corrupt the process. Do you agree with this?

But what stochastic differential equation do you solve with these manipulations of the Runge-Kutta steps ? You will get a result, but how to interprete it ? You make changes to the numerical scheme, but in reality, you have to make changes to the differential equation. After this, you have to derive or - if already existent - use an integration scheme suited for this stochastic differential equation.

[Paul]

@Bruno Luong

I hope you don't think I'm one of the "few of you to insist" "to apply pseudo-noise on top of classical runge kutta scheme" insofar as my very first comment pointed out one of the pitfalls of such an approach, i.e., calling randn (or rand, or whatever) on every call to the odefun. I can definitely see other problems with using that approach with ode45, et. al.

Thanks for the SDE links.

[Bruno Luong]

I'm relieved ;-)

[Paul]

@Bruno Luong or anyone else ...

Are you aware of or can provide reference to a nonlinear differential equation that has a known output in repsonse to a random input? By "known output" I mean that it can be described in some closed form fashion, perhaps by its probability density function or Fourier transform.

[Bruno Luong]

The most known is Kolmogorov's K41 theory. I do not follow the latest development, but you can start with such paper on simple Burgers 1D model https://hal.science/hal-03110850/document

[Jürgen]

The noise injection in the example is directly added in the ODE function for the Duffing oscillator. I don't feel comfortable with that in the case of an adaptive step solver like ode45. The randn function is activated at every call of the ODE function, which, if I am correct, happens four times per valid time step for ode45. My intuition is that this likely reduces the noise compared to the intended sqrt(D) value. I believe it should only be activated once: when the ode45 solver chooses the total step. However, this means one must modify the ode45 code. I considered using 'OutputFcn' to trigger the random generation, but this also does not work, as the valid duration of the step is only known after the step. The correct computation of the noise injection requires knowledge in advance to compute the magnitude of the random inhomogeneity injection, which must be derived from its spectral density.