How to generate colored Gaussian noise and adding it to a ODE system - Do I need Euler-Maruyama method?

Question

In the tutorial, when white noise process is added to ordinary differential equations (ODE), the ODE becomes a stochastic process. Then the stochastic process needs to be solved using Euler Maruyama method and not ODE.

I am having a hard time understanding how to generate and add colored noise in the form of process noise to a continous system such as the Rossler system. randn() is used to generate White Noise in Matlab and I directly added it to the equations. But I don't know if this is the correct way to add noise to a continous time system or not. This is what I did to add White Noise to the equations:

noise_x = randn(1);
noise_y = randn(1);
noise_z = randn(1);
y = [
     -x(2)-x(3) + noise_x;            % dx/dt
     x(1)+a*x(2) + noise_y;           % dy/dt
     b+x(3)*(x(1)-c) + noise_z;       % dz/dt
] ; 

color noise process using Euler-Maruyama Method (EM) seems to be the way to solve the stochastic system after the noise is added. But I cannot understand what is the color of the noise from this method. Looking at the code, there is no way to know if there is any color. After the color noise is generated by calling the function function [W,t]= colornoise(a,b,T,N) how to use the output W as an additive term in the Rossler system. This is not clear to me. Can somebody please provide a small code snippet on how to add White noise and color noise such as pink, red, violet to the Rossler system.

The following are my questions:

Question 1) How to generate colored noise in general?

Question 2) How to add the colored noise in the form of process noise to the Rossler system or any continuos time system and which method to use for solving the system -- ODE or Euler Maruyama? Is the noise added Rossler system stochastic or not?

The code below is for integrating the Rossler system using ODE45 solver.

clc;
clear all
%%%% Number of variable and initial conditions:

nbvar=3; 
xini=ones(1,nbvar)/10;

%%%% Time parameters:

trans=100;
tend=500;
tstep=0.01;


b=0.2;   % (default value for chaos) 


ttrans = [0:tstep:trans];
tspan = [0:tstep:tend];

x0=xini;
option = [];
[t x] = ode45(@dxdt,tspan,x0,option,b);

plot3(x(:,1),x(:,2),x(:,3));
xlabel('X','fontsize',18);
ylabel('Y','fontsize',18);
zlabel('Z','fontsize',18);
box on;


% ===================================================================
% dxdt
% ===================================================================

function y = dxdt(t,x,b)

%%% parameters

a=0.2;
%b=0.2;
c=5.7;

%%% equations

y = [
     -x(2)-x(3);            % dx/dt
     x(1)+a*x(2);           % dy/dt
     b+x(3)*(x(1)-c);       % dz/dt
] ; 

Answer 1

Question 1) How to generate colored noise in general?

Produce white noise and filter it.

Question 2) How to add the colored noise in the form of process noise to the Rossler system or any continuos time system

This is a little bit more complicated. The addition of noise is a simple matter of $y = x + q$ where $q$ is the noise signal. The key problem here is how much noise to add? This is important for linear cases (such as $y$ previously) because random number generators are producing numbers in some fixed interval ($[0 \ldots 1]$, $[-1 \ldots 1]$ or other) which will be entirely different than the interval the signal $x$ occupies.

In this case, assuming that the signal is stationary, then given some $x$, generate a $q$, scale $q$ so that it has the same power as $x$ and then mix them together at the desired ratio. For example $y = (1 - a) \cdot x + a \cdot q$. Here, the contribution of noise is directly proportional to $a$.

This is nice because you can now use $a$ later on in your experiments as a variable that shows how much noise was added to the signal.

But, when you are dealing with non-linear systems, this becomes even more important (and difficult). The reason for this is that, depending on the system, you cannot assume that it is stationary. Or, you can, but only at a narrow range of operation and furthermore, you are going to have to make sure that the system stays in this region.

The added difficulty is in two places:

1) You are dealing with differential equations that get integrated

• Therefore small spurious excitations here and there can accumulate to large numbers quickly

2) You are dealing with non-linear systems

• Therefore, combined with #1 it becomes likely that you might push the system into a region that its behaviour changes. For example, you might be interested to use the system when it produces nice smooth waveforms but inadvertently you might tip it over to the region where it starts producing spiking oscillations.

The reason I am mentioning this is because you refer to "process noise". So, you can take the Rossler model and integrate it and look at its signature attractor as an experiment or, you might want to use it as an actual model of a process, in which case, $x,y,z$ now have units, $dt$ is tied to your sampling frequency or other physical limit and also your noise now has a physical meaning too.

I don't know what these are in your case, but you have to understand them in order to be able to add the noise, especially if you have to add the noise in the derivatives (are you sure about this?).

The derivative is not the signal but rather the rate of change of the signal with respect to time. This is incredibly powerful. If you were to produce some $x$ by integrating the $\frac{dx}{dt} = b$, you would get a ramp in $x$. If your derivative now was a ramp itself with $\frac{dx}{dt} = c \cdot t + b$, your $x$ would vary exponentially.

Imagine then, how much you are impacting $x$ by specifying some random factor in its $\frac{dx}{dt}$.

If you don't have to add the noise at the derivatives, it would be slightly easier, because you simply integrate your differential equations to produce some output $z$ and then add noise to this specific $z$ as if it was our $x$ from above.

...and which method to use for solving the system -- ODE or Euler Maruyama? Is the noise added Rossler system stochastic or not?

If you have to add the noise to the derivatives then the system becomes stochastic and it would be better to use the Euler - Maruyama. The difference of Euler - Maruyama from the plain Euler integration is in the inclusion of a derivative term for the random process as well. Of course, if you differntiate white noise, you still get white noise, so from a practical point of view, there is no difference there because you still just have to add a random number. BUT this random number has to be appropriately scaled for its amplitude and then be multiplied by the integrating step too.

Hope this helps.