I am trying to iterate the following equation $$ x_{k}(n+1)=x_k (n)-\epsilon (x_{k+1}(n)-2x_k(n) +x_{k-1}(n))+\sqrt{\epsilon}\; \eta_{k}(n) $$ where $n$ denotes which time step I'm on and $k$ is the location on the string with periodic boundary conditions, $\epsilon$ is the size of the time step and $\eta$ is a random variate from a gaussian distribution with mean 0 and variance 1. $\eta$ is re-picked after each time step and is varies from location to location. I generated my initial $x$ randomly with

x[0] = RandomReal[{0, π}, 10]

which gives me 10 real numbers between 0 and $\pi$, which are my 1D lattice site variables. I then define $\eta$ by

η[0] = RandomVariate[NormalDistribution[0, 1], 10]

which gives me 10 random numbers from my distribution. I am taking the time step to be $\epsilon =1/100$ Then to define my iteration equation I write

x[[j]][n_ + 1] := x[[j]][n] - ϵ (x[[j + 1]][n] - 2 x[[j]][n] + x[[j + 1]][n]) + Sqrt[ϵ] η[[j]][n]

but this failed to iterate when I tried when just looking at a sigle site and its neighbors. If someone could be so kind as to help me put in periodic boundary conditions, as well as to get the iteration to work, I would be grateful. I am using Mathematica 8. Thank you.


2 Answers 2


In Mathematica I recommend do this with Sparse Matrix-Vector products as follows. First I will define a circulant matrix $C$ (with the periodic BCs) which appears on the RHS of your equation.

In matrix notation your equation reads \begin{equation} x(n+1) = x(n) + \epsilon C x(n) + \sqrt{\epsilon} \eta(n) \end{equation} or \begin{equation} x(n+1) = (I + \epsilon C) x(n) + \sqrt{\epsilon} \eta(n) \end{equation}

So, here are the matrix definitions. The classical (-1 2 -1) stiffness matrix is:

StiffnessMatrix[n_] := SparseArray[{Band[{1, 1}] -> 2., Band[{1, 2}] -> -1., 
   Band[{2, 1}] -> -1.}, {n, n}, 0.];

From this definition the circulant matrix is constructed with a simple adjustment.

CirculantMatrix[n_] := Block[{K = StiffnessMatrix[n]}, 
   ReplacePart[K, {{1, n} -> -1., {n, 1} -> -1.}]];

For a given $n$ (size of your system) we the define

C = CirculantMatrix[n];

The identity matrix can be computed with IdentityMatrix[n]. Conversely, a sparse implementation reads

Eye[n_] := SparseArray[{{i_, i_} -> 1.}, {n, n}, 0.];

So, to implement $A = (I + \epsilon C)$ we write

A = Eye[n] + eps C;

We are now ready to iterate:

x = x + A.x + sqe eta;

where $sqe = Sqrt[eps]$.

The following implementation is for a noise $\eta$ which is independent of position. You can adjust it for your needs.

Iteration[nmax_, eps_, size_] := With[{
   A = Eye[size] + eps CirculantMatrix[size],
   x0 = RandomReal[{0, \[Pi]}, size],
   sqe = N@Sqrt[eps]
    NestList[A.# + sqe RandomVariate[NormalDistribution[0, 1], size] &, x0, nmax]];

Input is the number of iterations nmax, the $\epsilon$ guy eps and the size of your system (which in you example is 10). Output is a matrix where each entry (e.g. sol[[1]]) is the state of your system at a given time. If you want only the last output, either use Last in your output or replace NestList with Nest.


It's been a while since I coded something in mathematica so I'll write in fortranish pseudocode instead

ntimesteps = 20
nx = 10
x = array(0:nx+1, 0:ntimesteps)
x[1:nx,0] = RandomReal(0,pi)
// boundary conditions:
x[0,0] = x[nx,0]
x[nx+1,0] = x[1,0]
do i=1,ntimesteps
  eta[1:nx] = RandomNormalDistribution(0,1)
  do j=1,nx
    x[j][i] = x[j,i-1] - eps(x[j-1,i-1]-2x[j,i-1]+x[j+1,i-1]) + sqrt(eps)eta[j]
  x[0,i] = x[nx,i]
  x[nx+1,i] = x[1,i]

Note that in the formula of the iteration you gave there are two mistakes: First use n_ instead of n_+1, secondly you have two j+1 instead of one j+1 and one j-1


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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