# # See the exercise # in http://www.physics.cornell.edu/~myers/teaching/ComputationalMethods/ComputerExercises/CardiacDynamics/CardiacDynamics.html # import scipy, scipy.optimize, time from DynamicLattice import DynamicLattice DISPLAYINTERVAL=10 def del2_5(A, dx): """del2_5(A, dx) returns a finite-difference approximation of the laplacian of the array A, with lattice spacing dx, using the five-point stencil: 0 1 0 1 -4 1 0 1 0 The no-flow boundary conditions mean that the sites on the boundary (the "ghost cells") are not really part of the simulation: they exist only to set the normal derivatives to zero. Our second derivatives will be zero along these ghost-cells, for convenience. You'll need to know about slicing. A slice of a list or array is a contiguous chunk of the array. For a one-dimensional list v, for example, v[3:6] = [v[3],v[4],v[5]], and for our matrix A, A[3:5,4:6] = [[A[3,4],A[3,5]],[A[4,4],A[4,5]]]. One can set the values of a whole slice at once: v[3:5] = [3.2,2.2] will set those two elements. For all arrays and lists, a negative argument denotes counting backward from the end of the list, so v[-1] is the last element of v. In a slice, and omitted argument is interpreted as continuing to the end of the list, so v[3:] = [v[3], v[4], ...] and v[:] is the entire list. Compute the laplacian by shifting the array A appropriately (up, down, left, and right) to implement the stencil. The resulting array returned should represent the laplacian of the array A on all of the interior points of the array A[1:-1, 1:-1]: your function should set del2 to zeros with the shape of A, and then del2[1:-1,1:-1] = (A[1:-1,2:] + ...)/(dx*dx) """ nx, ny = A.shape del2 = scipy.zeros((nx, ny), float) del2[1:-1, 1:-1] = (A[1:-1,2:] + A[1:-1,:-2] + \ A[2:,1:-1] + A[:-2,1:-1] - 4.*A[1:-1,1:-1])/(dx*dx) return del2 def del2_9(A, dx): """del2_9(A, dx) returns a finite-difference approximation of the laplacian of the array A, with lattice spacing dx, using the nine-point stencil: 1/6 2/3 1/6 2/3 -10/3 2/3 1/6 2/3 1/6 Compute the laplacian by shifting the array A appropriately (up, down, left, and right) to implement the stencil. The resulting array returned should represent the laplacian of the array A on all of the interior points of the array A[1:-1, 1:-1] """ nx, ny = A.shape del2 = scipy.zeros((nx, ny), float) del2[1:-1, 1:-1] = ((2./3.)*A[1:-1,2:] + (2./3.)*A[1:-1,:-2] + \ (2./3.)*A[2:,1:-1] + (2./3.)*A[:-2,1:-1] - \ (10./3.)*A[1:-1,1:-1] + \ (1./6.)*A[2:,2:] + (1./6.)*A[2:,:-2] + \ (1./6.)*A[:-2,2:] + (1./6.)*A[:-2,:-2])/(dx*dx) return del2 del2 = del2_9 def FindFixedPoint(gamma, beta): f = lambda v, gamma, beta: (v-(v**3)/3.)-((1./gamma)*(v+beta)) Vstar = scipy.optimize.brentq(f, -2., 2., args=(gamma, beta)) Wstar = ((1./gamma)*(Vstar+beta)) return Vstar, Wstar def CopyGhost(A): """ Given a 2D array A, create a zero-derivative boundary condition at all four edges of the array. This is done by changing the values along the boundary rows and columns to be equal to the second-from-the-boundary rows and columns. For example, to ensure the last column of A (in python, A[:,-1]) has a zero derivative, you copy the second-to-the-last column: A[:,-1] = A[:,-2] """ A[0,:] = A[1,:] A[:,0] = A[:,1] A[-1,:] = A[-2,:] A[:,-1] = A[:,-2] class FitzNag2D: """FitzNag2D is a class to support the simulation of the FitzHugh-Nagumo equations on a 2D square grid.""" def __init__(self, N=100, dx=1.0, eps=0.2, gamma=0.8, beta=0.7): """Initialize a FitzNag2D2D class, consisting of NxN points on a square grid with lattice spacing dx, and model parameters eps, gamma and beta; the FitzHugh-Nagumo equations are: dV_dt = laplacian(v) + (1./eps) * (v - (1./3.)*v**3 - w) dW_dt = eps*(v - gamma*w + beta) The __init__ method will need to store the values of the specified parameters, and create the fields v and w (as scipy arrays of the appropriate size). v and w should be initialized uniformly to their values at the fixed point (resting state), Vstar and Wstar, except for the initialization of a pulse of height self.pulseheight within a square of size 10 at the center of the simulation grid. The __init__ method should also set a member variable, e.g., self.pulseheight, that will indicate the height of the pulse generated by interactive mouse events. Animated displays will be taken care of by the DynamicLattice class which has been imported. The FitzNag2D class should initialize an instance of DynamicLattice within the __init__ method: self.DL = DynamicLattice((N, N), zmin=-2.0, zmax=2.0) This particular instance will animate dynamics on an NxN grid, creating grayscale images of a specified field based on the field value in the interval (zmin=-2.0, zmax=2.0). """ self.DL = DynamicLattice((N, N), zmin=-2.0, zmax=2.0) self.N = N self.dx = dx self.L = N*dx self.eps = eps self.gamma = gamma self.beta = beta self.pulseheight = 3.0 Vstar, Wstar = FindFixedPoint(gamma, beta) self.V = Vstar*scipy.ones((N, N), float) self.W = Wstar*scipy.ones((N, N), float) self.V[(N/2-5):(N/2+5), (N/2-5):(N/2+5)] = self.pulseheight self.t = 0. def rhs(self): """self.rhs() sets the instantaneous value of the right-hand-side of the FitzHugh-Nagumo equations, by creating two scipy arrays self.dV_dt and self.dW_dt, which describe the time evolution of the v and w fields, respectively. self.rhs() will be called by the step() method as part of the Euler time-stepping scheme.""" self.dV_dt = del2(self.V, self.dx) + \ (1./self.eps) * (self.V - (1./3.)*self.V**3 - self.W) self.dW_dt = self.eps*(self.V - self.gamma*self.W + self.beta) def step(self, dt): """self.step(dt) increments the fields v and w by dt. The first and last rows and columns of our arrays are "ghost cells", which are added to the boundaries in order to make it convenient to enforce the boundary conditions. Our no-flow boundary condition wants the normal derivative (derivative perpendicular to the boundary) to equal zero. step(dt) first uses CopyGhost to copy the boundary (ghost) cells to ensure the zero--derivative (no-flow) boundary condition. It then implements an Euler step, by calling self.rhs() to set the current values of the fields dV_dt and dW_dt, and then incrementing V by dt*dV_dt and similarly for W. You may wish to make dV_dt and dW_dt member variables of the FitzNag2D class. """ CopyGhost(self.V) self.rhs() self.V += dt * self.dV_dt self.W += dt * self.dW_dt def run(self, T, dt): """self.run(T, dt) integrates the FitzHugh-Nagumo equations for a time interval T, by taking steps of size dt. self.run() should also increment an overall time counter by dt. Every DISPLAYINTERVAL steps (e.g., 10), interaction with the self.DL object should be undertaken, to update the display of the v field, and to determine if a pulse box has been selected with the mouse. The display can be updated with the self.DL.display() method, which takes an array to be displayed as an argument (e.g., self.V). An additional nicety is to set the title of the self.DL window with the current simulation time, using the self.DL.setTitle() method. The self.DL object supports the query IsBoxSelected(), to ascertain whether a region (box) of the grid has been selected with the mouse. If a box has been selected, the self.DL.GetMouseBox() can be called to unpack a set of grid points (x0, y0, x1, y1) delineating the selected box. The v field within the selected box should be set to the pre-selected pulseheight. """ count=0 t0 = self.t while self.t < t0 + T: if count%DISPLAYINTERVAL==0: if self.DL.IsBoxSelected(): x0, y0, x1, y1 = self.DL.GetMouseBox() self.V[x0:x1, y0:y1] = self.pulseheight self.DL.setTitle('t = %f' % self.t) self.DL.display(self.V) time.sleep(0.05) self.step(dt) self.t += dt count +=1 def demo(): print("Cardiac Dynamics:") print(" Fitzhugh Nagumo in 2D") print(" Click and drag a rectangle for electroshock") f = FitzNag2D(N=100, dx=1.) f.run(300., 0.1) return f if __name__ == '__main__': f = demo() # Copyright (C) Cornell University # All rights reserved. # Apache License, Version 2.0