2630.sagews

This worksheet implements the Euler-Maruyama method for SDEs (see http://en.wikipedia.org/wiki/Euler-Maruyama_method for more details) and presents an interactive demo of a simple dynamical system exhibiting a Hopf bifurcation with noise (taken from http://johncarlosbaez.wordpress.com/2010/12/24/this-weeks-finds-week-308/).

Other uses of the Euler-Maruyama method can be found in the Sage interact/diffeq demo.  In fact, this worksheet was heavily inspired by that demo.

def euler_maruyama(a, b, x0, h, N):
    """
    Implements the Euler-Maruyama method for SDEs in 2D of the form 
    dX = a(X)dt + b(X)dW, where W is a vector of independent 
    Brownian motions.
    
    INPUT:
        
        - ``a``, ``b`` - vector-valued functions 
        - ``x0`` - initial conditions (2D vector)
        - ``h`` - step size of the method
        - ``N``  - number of steps to take
        
    OUTPUT: List of integration points.
    """
    
    output = [x0]
    for _ in range(N):
        W = vector([normalvariate(0, sqrt(h)), normalvariate(0, sqrt(h))])
        output.append(output[-1] + h*a(output[-1]) + b(output[-1]).pairwise_product(W))
    return output

The interactive demo allows you to control the parameter $\beta$ (determining the limit cycle) and the amplitude of the noise $\ell$.  When the noise is absent ($\ell = 0$), the origin is a stable fixed point when $\beta < 0$.  As $\beta$ becomes positive, the origin turns into an unstable fixed point and a stable limit cycle is born at $r = \sqrt{\beta}$. 

When noise is added to the system, the limit cycle can still be discerned for $\beta > 0$. 

For those of us schooled in classical ODEs, it is worth keeping in mind that these trajectories are continuous but not differentiable.  The EM solver computes only a finite number of points on each trajectory and then plots a straight line between the individual points, giving the appearance of a curve which is "somewhat smooth" for large step sizes.  The step size can of course be lowered, revealing the intricate features of the trajectories at the expense of longer computations.

@interact
def hopf_sde(beta=slider(srange(-2, 2, .1), default=1.0), ell=slider(srange(0, 1, .1), default=0.1)):      
    h = .02 # Step size      
    hopf_a = lambda x : vector([-x[1], x[0]]) + beta*x - x.norm()^2*x
    hopf_b = lambda x : vector([ell, ell])
    output1 = euler_maruyama(hopf_a, hopf_b, vector([0.0, 0.1]), h, int(10/h))
    output2 = euler_maruyama(hopf_a, hopf_b, vector([0.0, 2.0]), h, int(10/h))
    p =  list_plot([p.list() for p in output1], plotjoined=True, thickness=1, color='red')
    p += list_plot([p.list() for p in output2], plotjoined=True, thickness=1, color='green')
    p.show(xmin=-2, xmax=2, ymin=-2, ymax=2)

beta  ell  [removed] [removed] [removed] [removed]