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

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

• One of your questions is easy to answer: what's the color if the sequence generated by randn? Well, color if noise can easily be seen by the autocorrelation function. Think about what that looks like for a good random number generator! randn is inherently white. – Marcus Müller Nov 25 '17 at 10:58
• Oh and by the way, the Euler-Maruyama Matlab code snippet shows a lot of signs of bad coding practices. Among these: · Re-seeding the RNG at every call (this leads to the same noise being generated every time this function is called). Handling of random state must be up to the user. · unclear what it does. If you find code like this, the first reaction must be to read it, and for now, not use it. I don't really see how you need to throw Euler-Maruyama at the problem if you know the desired spectral shape of noise – take white noise and simply filter it accordingly. – Marcus Müller Nov 25 '17 at 13:15
• @MarcusMüller: Thank you for your feedback. I don't know how to create different colored noises which can be added to a continuous system. Many textbooks say that adding a random process to a continous system creates a stochastic process which needs to be numerically solved by the Euler-Maruyama or other sophisticated methods such as the Milstein method (one such paper is physik.uni-augsburg.de/theo1/hanggi/Papers/160.pdf) – SKM Nov 26 '17 at 2:23
• The link for the Euler-Maruyama method for color noise generation does not say what is the color of the noise. Using randn I will get a white colored noise. But how to generate color noise and add that color noise to a stochastic process is something which I have no clue about, hence I have posted this Questionfor help. Could you please answer the Question if possible?thank you – SKM Nov 26 '17 at 2:23
• I'm not sure what you mean with "needs to be…solved by…"; you're generating or modelling that system when you add noise, you're not solving something; anyway: That system is defined by a and b, which are values you feed in, and since I can't possibly tell you what you use, I can't tell you the color of your noise. They are then used to weigh an autoregressive term (in the for loop, Xtemp = Xtemp + a*Dt*Xtemp+…) and a moving average term (in the same term , b*Winc). As said, this is just a really really awkward way of filtering white noise, and you should probably avoid it. – Marcus Müller Nov 26 '17 at 10:13

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

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.

• Thank you for your reply. I strongly sense that in my problem the noise must be included in the derivatives. But I still don't quite understand how to implement the derivative term and the Euler-Maruyama method. Do you have any example code for a similar problem? If h=0.01;hspan = [0:h:40] is the integration step, and y3 is one of the state variable of the Rossler system to which I have added a noise of variance var =0.01, then is this how I would add noise?for i=2: length(hspan) [t,y] = ode45(ode, [0 h], y3(:,i-1));  y3(:,i) = y(end,:)'+sqrt(var)*[0 ; randn]*h; end  – SKM Nov 28 '17 at 19:13
• Glad you found it helpful, you can upvote and / or accept it via the controls on the right. To come back to your question, in brief, yes. There is some code here. The problem is not producing the random number, the problem is ensuring that you are adding the desired amount to the derivative. Within the ode45 integration step you produce instantaneous values. You have to ensure that these stay within specific limits over time. What other information about the process do you have? – A_A Nov 28 '17 at 19:22
• I will generate colored noise (any spectrum at this point) from a white gaussian noise of a small variance say 0.01. This colored noise is added to all the 3 state equations of the Rossler. I had seen the wiki code but cannot follow. The color generation code using Euler-Maruyama method that I included the link in my question is difficult to understand and I don't know how to apply that method to the Rossler. Since now it is clear that the system would become a stochastic, I don't know how to implement it using Matlab. Could you please extend your help in the implementation. Thank you – SKM Nov 28 '17 at 19:35
• I would like to give it a go but I can't promise it will be any time soon (i.e. tomorrow). I have not employed this integration method before (and maybe I should have at the time). You can find a similar thing in this paper. Look up fig.3 / eq 05 on p 1747. where they work out the coupling coefficients. That model is even more "sensitive" than the Rossler because if you are not careful you might push it into saturation where it stops oscillating altogether. – A_A Nov 28 '17 at 19:57
• Thank you for the paper but I think the Authors are synchronizing using coupling and noise. I have not reached that stage yet although I will be studying about synchronization methods and influence of noise. But before that I need to know how to apply the Euler Maruyama method for just one system. If you have an account in stackexchange, you could answer that question with the implementation. I am thinking of starting a bounty for that question as well. Do you think that the small code snippet which I methioned in the earlier comment is not the correct way for EM?Your help is much appreciated. – SKM Nov 28 '17 at 20:25