CompSciHW9

Introduction to Scientific Computing, HW 9. Due Friday, Apr. 15, 2016.

1) Use the van der Pol oscillator class to simulate 60 second (6000 steps) trajectories with 3 initial conditions, 
  x, v = [4, 3], [4, 4], and [4, 5]
  Here, v = dot x = dx/dt

2) Add a method to calculate the energy of a harmonic oscillator,
  E(x, v) = (x^2 + v^2)/2

3) Plot E(t) for all 3 simulations on the same axes.  Remember to scale t = arange(6000)*dt.

4) Is there a pattern that would let you know when the oscillator is near the steady-state?

Remember, we built the van der Pol integrator on CompSciWeek13. The complete code is below,

<source lang="python">

1. !/usr/bin/env python

import sys from numpy import * from numpy.random import standard_normal norm = standard_normal

class State:

   def __init__(self, x, v, mu, dt):
       assert x.shape == v.shape
       self.mu = mu
       self.dt = dt
       self.x = x
       self.v = v

   def step(self, N=1):
       for i in xrange(N):
          self.v += (self.mu*(1.0-self.x**2)*self.v - self.x)*self.dt
          self.x += self.v*self.dt
   def __str__(self):
       s = ""
       for x,v in zip(self.x, self.v):
           s += " %f %f"%(x, v)
       return s

def main(argv):

   assert len(argv) == 2, "Usage: %s <mu>"%argv[0]
   mu = float(argv[1])
   print mu
   dt = 0.01
   st = State(array([2.0]), array([0.0]), mu, dt)
   st.step(10)
   print st

if __name__=="__main__":

   main(sys.argv)

</source>