CompSciWeek13

Numerical Integration + Binary Output Format

• Van der Pol Oscillator: <math>\ddot u - \mu (1-u^2) \dot u + u = 0</math>
  • Implementations in OpenOffice Spreadsheet, GNU ODE, and Python
  • Perturbed initial conditions, Lyapunov exponents, and large deviation principles
• Working well with others: text
• Working with machines: binary
  • This is preferable when you have lots and lots of data
  • The relevant numpy methods are x.tofile("file") and x = fromfile("file") -- but save("file", x) and x = load("file") are preferred
  • inspecting binary formats with od
• Lecture Slides

GNU ODE Integrator

vanderpol.ode

<source lang="haskell"> mu = %MU

u = 1.0 y = 0.0

u' = y y' = mu*(1-u*u)*y - u

print t, u, y, y! every 10 from 50 step 0, 100 </source>

You can do the substitution with sed and run all dynamics with:

 for((i=0;i<10;i++)); do
   sed -e "s|%MU|$i|" vanderpol.ode >vanderpol-$i.ode
   ode -f vanderpol-$i.ode >out-$i.dat </dev/null
 done

Plotting with Gnuplot

<source lang="gnuplot"> set term postscript enh color lw 3 "Times-Roman" 28 set out 'van.eps'

set xlabel "u" set ylabel "y" plot 'out-0.dat' u 2:3 w l title "0", \

    'out-1.dat' u 2:3 w l title "1", \
    'out-4.dat' u 2:3 w l title "{/Symbol m} = 4"

</source>

Run with:

gnuplot plot_van.gplot

Python Integrator

First, form the heart of the algorithm: <source lang="python"> def step(N=1):

   for i in xrange(N):
       v += (mu*(1.0-x**2)*v - x) * dt
       x += v * dt

</source> Note that this will work if x and v are vectors as well as scalars. That means we can run a bunch of trajectories with different initial starting points in one fell swoop.

Next, wrap it inside a class: <source lang="python"> 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 state(self, N=1):

1. must change all variables (e.g. "name") to "self.name"

</source>

Next, instantiate and run the class:

<source lang="python"> st = State(5.0, 1.0, 2.0, 0.01) st.step() print st.x st.step() print st.x </source>

Next, add a nice printing function to "print the state," <source lang="python">

   def __str__(self):
       s = ""
       for x,v in zip(self.x, self.v):
           s += " %f %f"%(x, v)
       return s

</source> This prints (x1 v1 x2 v2 ...) in a single, long line.

OK, so now we're ready to run lots of tests. Let's lift the python program to a command-line program. Add this to the header: <source lang="python">

1. !/usr/bin/env python

import sys from numpy import * from numpy.random import standard_normal norm = standard_normal </source>

Add this to the end of the file: <source lang="python"> def main(argv):

   print argv

if __name__=="__main__":

   main(sys.argv)

</source>

From the shell, run <source lang="bash"> chmod +x vanderpol.py ./vanderpol.py a test argument </source>

Note that you get to provide parameters to your program at the command-line. Unfortunately, they all have type str. Let's formalize them, <source lang="python"> def main(argv):

   assert len(argv) == 5, "Usage: %s x0 p0 R mu"%(argv[0])
   x0, p0, R, mu = map(float, argv[1:5])
   st = State(x0,p0,mu,0.01)
   for i in range(500):
       st.step(10)
       print st

</source>

Finally, we have a new tool that can be run to generate reproducible simulations of the van der Pol oscillator. Next, we need to consider how to store the data it generates for later analysis and use.

Here's an example using a 500x2x5 tensor: <source lang="python"> T = zeros((500,2,5)) # steps by (x or v) by points for i in range(500):

   st.step(10)
   T[i,0] = st.x
   T[i,1] = st.v

</source>

This can be "flattened" to a 500x10 tensor, <source lang="python"> T = reshape(T, (500,10)) # or reshape(T, (-1,10)) "figures out" the 500 size </source>

Working with binary

Example: ELF header

od -t x1 -t c -N 8 /bin/bash

Binary comes in lots of units

8 bits	= 1 byte
4 bits	= 1 nibble
1 byte	= 1 ascii character
2 bytes	= 1 short = 1 word
4 bytes	= 1 32-bit int = 1 float
8 bytes	= 1 64-bit int = 1 double = 1 float complex
16 bytes	= 1 long double = 1 double complex

Byte ordering on most modern 64-bit processors is little-endian (Intel, AMD). For more info, see [1].

In an imaginary system where base 10 numbers are encoded by 4-bit nibbles, a 4-digit representation of the number 123 is

3 = 0011
2 = 0010
1 = 0001
0 = 0000

The binary together would read:

0011 0010 0001 0000

Since we have to use 4 digits, the most significant digit is set to all zeros, and goes last. The little end (least significant byte) goes first. Notice that within the 4-bit nibbles, the least-significant bit is still on the right. So the number is encoded using this left-to right (within a digit) and bottom to top (digits within a word) order.

Real binary encoding is just like the above, except that each 'digit' is a byte, so every 8 bits are in order, but the bytes are reversed.

Basis Functions

• Constructing B-splines - tensor method
• Representing the differentiation operator
• Solving PDEs using the implicit method
• Lecture Slides