Implementing a finite difference method in Mathematica

Question

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.

Answer 1

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]
  enddo
  x[0,i] = x[nx,i]
  x[nx+1,i] = x[1,i]
enddo

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

Answer 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.