In the tutorial, when white noise process is added to ordinary differential equations (ODE), the ODE becomes 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) shows an implementation of generating colored noise. 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.
Question 1) When to use the EM method and do I need it for my problem?
Question 2) Why can't I simply use randn() and not use EM method which is so complicated.
The code below is for integrating the Rosslersystem using ODE45 solver. Which part do I add the noise?
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
] ;