This is a summer project I did in 2018 summer in Imperial College London. The main task is to simulate the following logistic stochastic delay differential equation.

where and are parameters of the deterministic delay logistic equation, and .

We were also askd to simulate the system when a pullback in time is set, in order to find pullback attractors.

I am going to explain some of the algorithms I used in the scripts.

A pullback means we have t taking negative values. For example, if pullback is 200, then the system starts at t = -200. We want to make sure that the noise realisation should be the same at each time point, so we can adjust the value of pullback and explore pullback attractors.

Thus, in order to allow pullback to vary without altering the noise realisation, we should not generate noise forward from the very initial time point (e.g. t = -200). Instead, we should generate two streams of noise, starting from t = 0.

Stream A runs forward from t = 0 to .

Stream B runs backward from t = 0 to .

NOTICE : we need two different seeds to generate those streams, otherwise there is a symmetry in noise.

This makes sure that once the random seed is fixed, we can adjust pullback and the final time as we wish without altering the noise.

The philosophy of numeric integration is to discretise time, and use summation to replace integration.

Two integration schemes are used for integrations. In most of the scripts, Euler-Maruyama method are used to save computing time. Heun's Method is only used when stated in the title of the scripts. More sophisticated integration schemes like Runge-Kutta, requires fractional time step, which I found infeasible for stochastic delay differential equation.

The Euler-Maruyama method is basically a stochastic version of the Euler's method for deterministic equation. Under Euler-Maruyama method, our equation becomes

follows a normal distribution with variance . Thus, is realised by drawing a sample from the normal distribution with variance . In fact, in my implementation, a smaller times step called tDelta is set, and , where is an integer. Now

Heun's Method is supposed to be more accurate than Euler's Method for integrating the deterministic equation, but it is more time- consuming. The scheme for integrating the random variable is the same as Euler-Maruyama, as Heun's Method only improves evaluation of the deterministic gradient.

For the sake of simplicity, suppose our equation is

Then, under Heun's method, we first use the same technique as Euler's method to find the value of X(t+dt) by

However, this is only an intermediate value. The purpose is to use this to evaluate and then use the average gradient to evaluate again.

The Lyapunov specturm is very useful in determining types of attractors. Negative Lyapunov spectra mean stable attractors, and positive Lyapunov spectra mean unstable attractors, and Lyapunov spectra with both positive and negative values mean strange attractors.

The dynamics we are interested in is

We use finite points to represent .

For the purpose of more robust computation, we use points, which are evenly spaced in the time domain, to represent the function, where is set to guarantee that is an integer. Now the delay function is approximated by a -dimensional object. Thus, this -dimensional object can be represented by a -dimensional vector, indicating the values the object takes at each time point. (i.e. when , where is an integer and )

Since we are dealing with -dimension, we choose the canonical basis . Aligning them together gives us the identity matrix. To estimate the Lyapunov spectrum, we are going to see how this canonical basis develops when time span is very large, under the linearisation of the of stochastic delay differential equation. Since this system is autonomous, the dynamics of the tangent equation can be written in the form

Thus starting with , we have . The Lyapunov spectrum is a vector :

where is the vector of all real parts of eigenvalues of .

Working from first principles, the tangent equation is found to be

Now by implementing the integration scheme, we can find , and estimate the Lyapunov Spectrum by choosing a large positive t .

The invariant measure could be understood as a (probability) distribution of particles in the dynamical system when time tends to infinity. Two methods are used to cross-validate the calculation of the invariant measure.

The Ulam's method uses Markov Chain.

Suppose we have a discrete finite state process, then the invariant measure could be calculated by computing the eigenvalues and eigenvectors of the corresponding Markov matrix. The normalised eigenvector corresponding to eigenvalue 1, is the required invariant measure.

If we have , we could partition the interval into N equally spaced subintervals. Then, considering each subinterval as a single state, we have a N-state Markov process.

And the Markov matrix of this process is estimated by:

However the process we are dealing with is

What should we do?

We have to map onto , so we can define a Markov process that we can work with. A natural choice is an equivalence relation on :

Thus is the required mapping which maps onto

We have to restrict to a finite interval . This is totally for the purpose of numerical simulation. We choose:

b = the largest value a path has ever visited

a = the smallest value a path has ever visited

Now, we have reduced our problem to what we have discussed at the very beginning of this section. So then we estimate the the Markov matrix and compute its normalised eigenvector with eigenvalue 1.

The Birkhoff's Theorem states that given is a probability space and a measurable function such that is invariant with respect to f then we have:

This means that if is finite, then

We use the same mapping to reduce our infinite dimensional object to one dimension with finite states as what we have done in Ulam's method.

Apart from the Logistic equation, the same scripts were also written for other two equations, which are the Predator-Pray equation, and Wreight's conjecture.

Predator-Pray:

Wreight's conjecture:

That is all.