Lecture 2 - Random processes

by Leonid Pryadko

Interpolated functions

Consider in more detail the structure of a numerical solution of a differential equation

sol=NDSolve[{y'[x]+y[x]==0,y[0]==1},y[x],{x,0,1}]

FullForm[y[x]/.sol[[1]]]

Such a function can be constructed directly using a table of values

tabl=Table[{x,Sin[x]},{x,0.,Pi,Pi/10.}]

This defines an interpolating function based on the table of values

f[x_]:=Interpolation[tabl][x]

You can use it directly

f[0.3]-Sin[0.3]

Plot[{Sin[x],f[x]},{x,0,Pi}]

Here is  the plot of a difference:

Plot[Sin[x]-f[x],{x,0,Pi}]

You can use interpolating functions anywhere in place of ordinary numerical functions

NIntegrate[f[x],{x,0,Pi}]

Accuracy can be improved if we take more points

tabl=Table[{x,Sin[x]},{x,0.,Pi,Pi/100.}]; 
f[x_]:=Interpolation[tabl][x];
NIntegrate[f[x],{x,0,Pi}]

The error is very good:

2-%

Random functions

Random[NormalDistribution[0,1]]

This defines normally distributed random number generator.

tabl=Table[Random[NormalDistribution[0,1]],{x,10}]

Here are the corresponding  Mean and Variance

{Mean[tabl],Variance[tabl]}

Try a bit more points

tabl=Table[Random[NormalDistribution[0,1]],{x,10000}];{Mean[tabl],Variance[tabl]}

We can now try to create a randomly distributed function.  First, try to use a table of random numbers

tabl=Table[{x,Random[NormalDistribution[0,1]]},{x,0,1,0.02}];
f[x_]:=Interpolation[tabl][x];
Plot[f[x],{x,0,1}]

This is not a function one could use reliably because the interpolation is a bit unpredictable.  Instead, let us smoothen these values by forming a convolution with a Gaussian (or any other suitable function):
∫ G(t-t') f(t') d t' = (Subscript[f, smooth])(t)
The corresponding correlator is easily found for ⟨f(t') f(t'')⟩=δ(t'-t'')
g(x-y)=⟨ (Subscript[f, smooth])(x)  (Subscript[f, smooth])(y)⟩ = ∫ G(x-t') d t' ∫ G(y-t'')  d t' ⟨f(t') f(t'')⟩ =∫ G(x-t') G(y-t')  d t'.

G[t_,τ_]:=(2/Pi)^(1/4)((Exp[-t^2/(τ^2)])/(Sqrt[τ]));(* Define the weight function *)

Integrate[G[(Subscript[t, 1]),τ] G[t-(Subscript[t, 1]),τ],{(Subscript[t, 1]),-Infinity,Infinity},Assumptions->τ>0]

First, generate a table of values for the weighting function

Gn=Table[ G[i* 0.01,0.05],{i,-25,25}];

ListPlot[Gn]

Make a table of random numbers

x=Table[Random[NormalDistribution[0,1]],{i,1000+51+1}];

ListPlot[x]

Convolution of the two lists

y=Sqrt[0.01] ListConvolve[Gn,x];

ListPlot[y]

Corresponding time values

t=Table[0.01(i-1),{i,Length[y]}];

And the table of pairs:

lis=Transpose[{t,y}] ; lis//Shallow

Corresponding random function is smooth on short time scale

f=Interpolation[lis]; Plot[f[x],{x,.0,.1}]

This is the corresponding correlation function

NIntegrate[f[x]f[x],{x,0,10},MaxPoints->100000,MaxRecursion->1000,MinRecursion->10,WorkingPrecision->40,PrecisionGoal->5,Compiled->False]

fff[y_]:=NIntegrate[f[x]f[y+x],{x,0,10-y},MaxPoints->10000,MaxRecursion->1000,MinRecursion->10,WorkingPrecision->40,PrecisionGoal->5]/(10-y)

fff[0.1]

outtab=Table[{x,fff[x]},{x,0.,.1,0.02}]

ListPlot[outtab]

Plot[Exp[-y^2/(2 0.05^2)],{y,0,0.1}]

Show[%,%%]

Put everything together as a simple function

ndCorrTab[Dt_,dt_,tau_]:=Block[{Gn,M=Ceiling[ Dt/dt], N=5tau/dt,x,y,t}, (* arguments: 
Dt: time interval, dt: time step, tau is the correlation time *)
(* variables: 
M: how many total values to output, 
N is the number of values for the weighting function *)
Gn=Table[ G[i* dt,tau],{i,-N,N}]; (* Table of the weighting function values *)
x=Table[Random[NormalDistribution[0,1]],{i,M+2N+1}]; (* random numbers *)
y=Sqrt[dt] ListConvolve[Gn,x]; (* convolution between the two; Sqrt[dt] needed for normalization *)
t=Table[dt (i-1),{i,Length[y]}]; (* corresponding time values *)
Transpose[{t,y}] (* return the list of pairs (t,y) *)
];

ndCorrTab[1,0.10,0.1]

This defines a random function over a time interval [0,100] with the time step 0.05 and the rather long correlation time τ=1.0

f=Interpolation[ndCorrTab[100,0.05,1.0]];

f[0.5]

Plot[f[t],{t,0,10}]

Stochastic differential equation

Clear values for x,y,t

x=.;y=.;t=.

sol=NDSolve[{y''[t]+0.1 y'[t]+y[t]==f[t],y[0]==0,y'[0]==0},y[t],{t,0,100}]

Plot[Evaluate[{y[t]^2,y[t]}/.sol[[1]]],{t,0,100},PlotStyle->{{Dashing[{0.02}]},{}}]

Same equation without friction shows linearly increasing square of the amplitude, not surprisingly!

sol=NDSolve[{y''[t]+0.0 y'[t]+y[t]==f[t],y[0]==0,y'[0]==0},y[t],{t,0,100}]

Plot[Evaluate[{y[t]^2,y[t]}/.sol[[1]]],{t,0,100},PlotStyle->{{Dashing[{0.02}]},{}},PlotRange->All]