Code reorganization

[ronojoy]

@vavachan, please do the following reorganization of your code

• write the Wilkie stochastic integrator as a another function, or better, as a class.
• write the function f(x, v) which returns the RHS of the ode system in a separate file and pass this to integrator class. Take a look at odespy for an example of this.
• write all the remaining code in an Ipython notebook, together with comments on what is being done
• use the acor package to compute the correlation function

[vavachan]

I have seperated the program into 3 , wilkie_stochastic.py ,function.py, const.py. wilkie_stochastic.py has the class wilkie_stoch() In function.py, functions A(X),B(X),F(X) and B_der(X) return matrices corresponding to matrices in equation no (7) in wilkie paper when written in matrix form. const.py has all the constants in one place. Since I am using the Computer center, I could not use acor or ipython,or acor but i will figure it out soon.

[ronojoy]

@vavachan, first study the odespy ordinary differential equation API example :

def f(u, t):
    """2x2 system for a van der Pool oscillator."""
    return [u[1], 3.*(1. - u[0]*u[0])*u[1] - u[0]]

import odespy, numpy
solver = odespy.Vode(f, rtol=0.0, atol=1e-6,
                     adams_or_bdf='adams', order=10)
solver.set_initial_condition([2.0, 0.0])
t_points = numpy.linspace(0, 30, 150)
u, t = solver.solve(t_points)

u0 = u[:,0]
from matplotlib.pyplot import *
plot(t, u0)
show()

It does the following things

• defines the function to be integrated
• sets up a solver object, which is passed the function to be integrated
• sets initial conditions
• defines the time points at which the solution is desired
• calls ''solve'' method to get the actual solution

This is the type of API I want you to construct, even if you have only one solver method now (i.e. the Wilkie integrator). Note that the main difference between an ode and a sde is that the sde has an additional "diffusion" term. The general sde is of the form

dx = (drift) * dt + (diffusion) * dW

When dW or the drift is zero, the sde reduces to an ode. Therefore, the API call for an sde should take the following form.

def drift(x, t):
    """ Drift for inertial Brownian motion in a harmonic potential"""
    return [-x[0], -x[1]]

def diffusion(x, t):
  """ Diffusion for inertial Brownian motion in a harmonic potential"""
   # write code to return 2x2 diffusion matrix    
   return D

import sde             # assume that you call your module sde
import numpy

# assume that the module has a class called ito 
# for solving the ito sde dx = (drift)*dt  + diffusion*dW

solver = sde.ito(drift, diffusion, method='wilkie-rk4', ) 
solver.set_initial_condition([2.0, 0.0])
t_points = numpy.linspace(0, 30, 150)
x, t = solver.solve(t_points)

x0 = x[:,0]
from matplotlib.pyplot import *
plot(t, u0)
show()

Rewrite your class file to reflect this structure. Have a look at the odespy codebase if you need to.

[ronojoy]

@vavachan, you can implement the simplest Euler-Maruyama method instead of Wilkie's integrator. Please have a look at the last chapter of Gardiner for an overview of numerical methods for stochastic differential equations.

[vavachan]

I have uploaded two files stoch_api.py and sde.py. I have some confusion regarding what matrices to sent to ito class, I have written the details as a comment above the stoch_api.py. For runge-kutta scheme the jacobian of B matrix is needed so I have made a 3rd matrix, to be sent to class ito.
Also the position(t) is not smooth when found by euler-maruyama scheme.

[ronojoy]

sde.py is moving in the right direction. Lets just stick to Euler-Maruyama for now, as we do not need to enter into the complication of the derivative of B for the problem of Brownian motion near a wall.

Regarding your comments in stoch_api.py, recall that the diffusion is a matrix, not a vector. So, there is

# B_vv(x,v)=sqrt(2*gamma*KB*T)/m
# B_xv = B_vx = 0
# B_xx(x,v)=0

The noise is a contraction of this matrix with the Wiener increment vector, that is

f_ = B_ij dW_j

For the E-M method, you only need the Wiener increments. For higher order methods, you will need stochastic double integrals, which is why I suggest we do this after the simpler E-M method is benchmarked to satisfaction.

Try reducing dt when integrating. The ratio of dt to the momentum relaxation time should be at least 10. Try to match the numerical results with the analytical expressions given in Chandrasekhar's Review of Modern Physics article on stochastic processes.

[vavachan]

I modified the program to include more parameters like more than 1 noise term, I have the example 3 given in joshua wilkie's paper,as default equation, it has 3 noise terms and 2 variables.

For the langevin equations, for both schemes MSD saturates to the theoretical value, but the fluctuations around it don't die out,(for both underdamped and overdamped parameters) similarly for correlation functions.

But both schemes seem to work very well for two examples given in joshua wilkie and another equation I found in a paper.

![Uploading em_omega_2_gamma_5.png…]() ![Uploading rk_omega_2_gamma_5.png…]()

[ronojoy]

@vavachan not sure what you mean by "fluctuations around it don't die out" - please upload all figures for the Langevin equation in a new issue thread.