Stochastic Phase Diffusion

Stochastic phase diffusion process is described as

\(x(t) = A\cos(\omega_0 t + \phi(t))\)

where \(\phi(t)\) is the stochastic phase. Assume the phase to be a free diffusion process (a.k.a. Wiener process) defined as

\(\phi \sim \mathcal{N}(0,\sigma^2), \sigma^2 = 2Dt, \dot \phi = sqrt(2D)\xi(t), <\xi(t)> = 0, <\xi(t)\xi(t+\tau)> = \delta(\tau)\)

Initialization

TIME_GRID <- c(0.001,10001,0.) #step, N, time_zero
PARAMETER <- c(10.,100.,1.) #D, omega_zero, amplitude
INITIALIZE <- c(0.) #phi_zero
MAX_TAU <- 100
ITER <- 2
SHOW <- 100
set.seed(1234)

STEP <- TIME_GRID[1]
N <- TIME_GRID[2]
TIME_ZERO <- TIME_GRID[3]
D <- PARAMETER[1]
OMEGA_ZERO <- PARAMETER[2]
AMP <- PARAMETER[3]
PHI_ZERO <- INITIALIZE[1]

Functions

Noise_Generator<-function(D,N,PHI_ZERO,STEP){
     sd <- (2*D)^0.5
     phi_dot <- rnorm(n = N,mean = 0.,sd = sd)
     phi <- c(PHI_ZERO)
     for(i in 1:(N-1)){
          x <- phi_dot[i]*STEP + phi[i]
          phi <- c(phi,x)
     }
     phi
}

Phase_Diffusion<-function(A,O,T,P,N){
     x <- c()
     for(i in 1:N){
          xx <- A*cos(O*T[i] + P[i])
          x <- c(x,xx)
     }
     x
}

Prob<-function(PHI,TIME,D){
     sd <- (2*D*TIME)^0.5
     p <- dnorm(x = PHI,mean = 0.,sd = sd)
     p
}

Process: \(x(t)\) and \(\phi(t)\)

TIME <- seq(from=TIME_ZERO,by=STEP,length=N)
PHI <- Noise_Generator(D,N,PHI_ZERO,STEP)
X <- Phase_Diffusion(AMP,OMEGA_ZERO,TIME,PHI,N)
par(mfrow=c(1,2))
plot(TIME[1:SHOW],PHI[1:SHOW],type="l")
plot(TIME[1:SHOW],X[1:SHOW],type="l")

Autocorrelation: \(R_{xx}(\tau)\)

TAU <- seq(from=STEP,by=STEP,length=MAX_TAU)
a <- 0.5*AMP**2
a <- rep(x = a,length=MAX_TAU)
z <- rep(x = D,length=MAX_TAU)
b <- exp(-z*abs(TAU))
c <- cos(OMEGA_ZERO*TAU)
Rxx_theory <- (0.5*AMP**2)*exp(-D*abs(TAU))*cos(OMEGA_ZERO*TAU)

Rxx <- c()
for(i in 1:MAX_TAU){
     sum <- 0.
     for(j in 1:(N-i)){
          sum <- sum + X[j]*X[j+i]*Prob(PHI[i],TIME[i],D)
     }
     sum <- sum/j
     Rxx <- c(Rxx,sum)
}

par(mfrow=c(1,1))
AA <- 0.5*AMP^2
plot(TAU,Rxx,type="l",main="Aucorrelation function of position: Rxx",ylim = c(-AA,AA))
lines(TAU,Rxx_theory,col="green")
legend(max(TAU)*0.7,AA*0.9,c("Simulation","Theory"),lty=c(1,1),lwd=c(2.5,2.5),col=c("black","green"))

NOTE: theoretically,

\(R_{xx}(\tau) = \frac{A^2}{2}e^{-D\mid\tau\mid}\cos(\omega_0\tau)\)