Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- # -*- coding: utf-8 -*-
- """
- Spyder Editor
- This is a temporary script file.
- """
- #!python
- """
- brownian() implements one dimensional Brownian motion (i.e. the Wiener process).
- """
- # File: brownian.py
- from math import sqrt
- from scipy.stats import norm
- import numpy as np
- def brownian(x0, n, dt, delta, out=None):
- """\
- Generate an instance of Brownian motion (i.e. the Wiener process):
- X(t) = X(0) + N(0, delta**2 * t; 0, t)
- where N(a,b; t0, t1) is a normally distributed random variable with mean a a
- nd
- variance b. The parameters t0 and t1 make explicit the statistical
- independence of N on different time intervals; that is, if [t0, t1) and
- [t2, t3) are disjoint intervals, then N(a, b; t0, t1) and N(a, b; t2, t3)
- are independent.
- Written as an iteration scheme,
- X(t + dt) = X(t) + N(0, delta**2 * dt; t, t+dt)
- If `x0` is an array (or array-like), each value in `x0` is treated as
- an initial condition, and the value returned is a numpy array with one
- more dimension than `x0`.
- Arguments
- ---------
- x0 : float or numpy array (or something that can be converted to a numpy arr
- ay
- using numpy.asarray(x0)).
- The initial condition(s) (i.e. position(s)) of the Brownian motion.
- n : int
- The number of steps to take.
- dt : float
- The time step.
- delta : float
- delta determines the "speed" of the Brownian motion. The random variabl
- e
- of the position at time t, X(t), has a normal distribution whose mean is
- the position at time t=0 and whose variance is delta**2*t.
- out : numpy array or None
- If `out` is not None, it specifies the array in which to put the
- result. If `out` is None, a new numpy array is created and returned.
- Returns
- -------
- A numpy array of floats with shape `x0.shape + (n,)`.
- Note that the initial value `x0` is not included in the returned array.
- """
- x0 = np.asarray(x0)
- # For each element of x0, generate a sample of n numbers from a
- # normal distribution.
- r = norm.rvs(size=x0.shape + (n,), scale=delta*sqrt(dt))
- # If `out` was not given, create an output array.
- if out is None:
- out = np.empty(r.shape)
- # This computes the Brownian motion by forming the cumulative sum of
- # the random samples.
- np.cumsum(r, axis=-1, out=out)
- # Add the initial condition.
- out += np.expand_dims(x0, axis=-1)
- return out
- -------------------------------------------------------------------------------------------------------------------------
- # -*- coding: utf-8 -*-
- """
- Created on Wed Jul 1 02:34:50 2015
- @author: Raleigh
- """
- import numpy
- from pylab import plot, xlabel, ylabel, title, grid, show
- from brownian import brownian
- #This is the code that performs the main iteration of the Euler Marayuma process, from t=0 to t=10.
- def main():
- # The Wiener process parameter.
- delta = 2
- # Total time.
- T = 10.0
- # Number of steps.
- N = 500
- # Time step size
- dt = T/N
- # Number of realizations to generate.
- m = 1000
- # Create an empty array to store the realizations.
- x = numpy.empty((m,N+1))
- # Initial values of x.
- x[:, 0] = 0
- brownian(x[:,0], N, dt, delta, out=x[:,1:])
- t = numpy.linspace(0.0, N*dt, N+1)
- for k in range(m):
- plot(t, x[k])
- xlabel('t', fontsize=16)
- ylabel('x', fontsize=16)
- title('Brownian Motion of 1000 particles with x_0 = 0')
- grid(True)
- show()
- if __name__ == "__main__":
- main()
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement