What is stochastic differential equation and its need?

Question

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?

Answer 1

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,:)));