0
$\begingroup$

A white noise process can be simulated using the Matlab command randn(). The numbers will be drawn from a Normal distribution of zero mean and variance 1. Is the mathematics behind generating white noise a stochastic differential equation which is numerically solved by stochastic Euler and Milstein methods ? A web tutorial https://me.ucsb.edu/~moehlis/APC591/tutorials/tutorial7/node3.html

presents the maths. What I could remotely follow is that the Euler method is used to convert a continuous system to discrete, I mean that the time to discrete. But in programming, we just call randn() to generate white noise.

For example, considering a model known as the Lorenz system. The model is a system of three ordinary differential equations now known as the Lorenz equations:

\begin{align} \dot{x} &= \sigma (y - x) \\ \dot{y} &= x (\rho - z) - y \\ \dot{z} &= x y - \beta z \end{align} .

If I want to add white noise to the $x$ variable: $$\dot{x} = \sigma (y - x) + w$$ then how would I do it? Should I generate w using w = randn(1,N) and then use this output to x? Where will I use Euler method and the stochastic differential equation? Why use the stochastic differential equation to generate white noise when we can simply use the programming tools or is it illegal to use the programming tools and the correct method is to generate white noise using the math?

$\endgroup$
2
$\begingroup$

One way to approach your problem is to use an Extended Kalman Filter with a small step size. $x,y,$ and $z$ are your state variables. Runge Kutta works better than Euler. The ODE solvers in Matlab also work, but need to be set up for a single small time step. It looks like at each time step, you propagate the differential equations forward in time, and the perturb the new states with noise, and then propagate those, repeat, ...repeat....

Google ( Lorenz Equation Matlab ) for code examples. You have a Noisy Lorenz system.

$$ \left| \begin{array}{c} x[k+1] \\ y[k+1] \\ z[k +1] \end{array} \right|= \text{ODE STEP} \left(\left| \begin{array}{c} x[k] + \text{noise}_x \\ y[k] + \text{noise}_y \\ z[k]+ \text{noise}_z \end{array} \right| \right) $$

Noisey Van Der Pol Example,

This is an old homework problem with some unnecessary feature like measurement noise, but you can modify it to your needs.

There is a way for the Matlab ODE solver to only iterate from 0 to h instead of calculating 41 sub intervals. but I don't recall how after a few years

You can use a 4th order Runge Kukta https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods

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

You should note that, given your problem statement, you don't really have a pure SDE, The Lorenz is a chaotic system, with a noise perturbation.

function dydt = vanderpoldemo(t,y,Mu)
%VANDERPOLDEMO Defines the van der Pol equation for ODEDEMO.

% Copyright 1984-2002 The MathWorks, Inc.
% $Revision: 1.2 $ $Date: 2002/06/17 13:20:38 $

dydt = [y(2); Mu*(1-y(1)^2)*y(2)-y(1)];

% HW 8 part c
clear all
h=.01;
tspan = [0:h:60];
y0 = [1; 0];
epsilon = [2 1 .2] ;
q=5;
r=1;

y1=zeros(2,length(tspan));
F1=zeros(2,2,length(tspan));
y1(:,1)=y0;
ode = @(t,y) vanderpoldemo(t,y,epsilon(1));
for i=2: length(tspan)
[t,temp] = ode45(ode, [0 h], y1(:,i-1));
y1(:,i) = temp(end,:)'+sqrt(q)*[0 ; randn]*h;
end
y2=zeros(2,length(tspan));
y2(:,1)=y0;
ode = @(t,y) vanderpoldemo(t,y,epsilon(2));
for i=2: length(tspan)
[t,temp] = ode45(ode, [0 h], y2(:,i-1));
y2(:,i) = temp(end,:)'+sqrt(q)*[0 ; randn]*h;
end
y3=zeros(2,length(tspan));
y3(:,1)=y0;
ode = @(t,y) vanderpoldemo(t,y,epsilon(3));
for i=2: length(tspan)
[t,temp] = ode45(ode, [0 h], y3(:,i-1));
y3(:,i) = temp(end,:)'+sqrt(q)*[0 ; randn]*h;
end

figure(1)
plot(tspan,y1(1,:))
[plot_limits]=ginput(2) %select 3 orbits
startp=round(plot_limits(1,1)/h);
endp=round(plot_limits(2,1)/h);
y1=y1(:,startp:endp);
t1=tspan(startp:endp);
h1=plot(y1(1,:),y1(2,:));
set(get(h1,'Parent'),'XLim',[-4 4],'YLim',[-4,4]);
xlabel('y(1)')
ylabel('y(2)') title([' Noisy van der Pol Equation, \epsilon = ',num2str(epsilon(1))])
%
%Plot of the second solution
figure(2)
plot(tspan,y2(1,:))
[plot_limits]=ginput(2) %select 3 orbits
startp=round(plot_limits(1,1)/h);
endp=round(plot_limits(2,1)/h);
y2=y2(:,startp:endp);
t2=tspan(startp:endp);
h2=plot(y2(1,:),y2(2,:));
set(get(h2,'Parent'),'XLim',[-4 4],'YLim',[-4,4]);
xlabel('y(1)')
ylabel('y(2)') title(['Noisy van der Pol Equation, \epsilon = ',num2str(epsilon(2))])
%
% Plot of the third solution
figure(3)
plot(tspan,y3(1,:))
[plot_limits]=ginput(2) %select 3 orbits
startp=round(plot_limits(1,1)/h);
endp=round(plot_limits(2,1)/h);
y3=y3(:,startp:endp);
t3=tspan(startp:endp);
h3=plot(y3(1,:),y3(2,:));
set(get(h3,'Parent'),'XLim',[-4 4],'YLim',[-4,4]);
xlabel('y(1)')
ylabel('y(2)') title(['van der Pol Equation, \epsilon = ',num2str(epsilon(3))])
z1=y1(1,:)+sqrt(r)*randn(size(y1(1,:))); % add measurement noise
z2=y2(1,:)+sqrt(r)*randn(size(y2(1,:)));
z3=y3(1,:)+sqrt(r)*randn(size(y3(1,:)));

$\endgroup$
  • $\begingroup$ No, you are simply adding noise to the lorenz states after you let the ode do it's steps. you need to add noise at each step, and propagate the noisey state $\endgroup$ – Stanley Pawlukiewicz Oct 18 '17 at 20:59
  • $\begingroup$ I cannot follow your steps which you have mentioned. In your answer, the noise is included in the ODE step and is the noise term generated using randn()? Can you please put the code snippet for say just one variable which will help me to get a clear picture? Thank you $\endgroup$ – Srishti M Oct 18 '17 at 21:15
  • $\begingroup$ I help people solve problems, not solve them for them. I'm sorry that you are still confused. randn() produces pseudorandom independent Gaussian numbers. You should also tag homework if it is homework $\endgroup$ – Stanley Pawlukiewicz Oct 18 '17 at 21:24
  • $\begingroup$ I have an example for a noisy van der pol oscillator, somewhere. I'll look for it tomorrow and will add it as an edit to my answer, assuming I can find it, which is likely. $\endgroup$ – Stanley Pawlukiewicz Oct 18 '17 at 22:08
  • $\begingroup$ Van Der Pol example added to answer $\endgroup$ – Stanley Pawlukiewicz Oct 19 '17 at 17:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.