MTH 337: Week 6

In [1]:
%pylab inline

Populating the interactive namespace from numpy and matplotlib

"Mean" finds the arithmetic mean of all the values in an array

In [2]:
x = arange(1, 11)
print x

[ 1  2  3  4  5  6  7  8  9 10]
In [3]:
mean(x)
Out[3]:
5.5

Differential Equations

We suppose that a system is described by a state $\vec{x}$ and the differential equation $\dfrac{d\vec{x}}{dt} = \vec{f}(\vec{x}(t), t)$ as follows: $$\vec{x} = \begin{bmatrix}x_1 \\ x_2 \\ \vdots \\ x_n \end{bmatrix}, \hspace{.5in} \dfrac{d\vec{x}}{dt} = \begin{bmatrix}\frac{d x_1}{dt} \\ \frac{d x_2}{dt} \\ \vdots \\ \frac{d x_n}{dt} \end{bmatrix} = \vec{f}(\vec{x}(t), t) = \begin{bmatrix}f_1(\vec{x}(t), t) \\ f_2(\vec{x}(t), t) \\ \vdots \\ f_n(\vec{x}(t), t)\end{bmatrix}$$

A Planet Orbiting a Star in a Plane

The state of the system is given by its position (2 components, $x$ and $y$) and velocity (2 components, $v_x$ and $v_y$).

From Newton's Second Law, and the Law of Universal Gravitation, the differential equation can be derived:

$$\vec{x} = \begin{bmatrix}x \\ y \\ v_x\\ v_y \end{bmatrix}, \hspace{.5in} \dfrac{d\vec{x}}{dt} = \begin{bmatrix} v_x \\ v_y \\ \dfrac{-G M x}{(\sqrt{x^2 + y^2})^3} \\ \dfrac{-G M y}{(\sqrt{x^2 + y^2})^3} \end{bmatrix}$$

Euler's Method

Euler's method gives a way to calculate the new state of the system after some small time $h$:

$\vec{x}(t + h) = \vec{x}(t) + h \dfrac{d\vec{x}}{dt}$

To implement Euler's Method for the planetary system, we:

  • Use a 1D NumPy array containing 4 elements for the state.
  • Define a function that returns the derivative of the state.
  • Keep track of the changing state through time in a 2D NumPy array.
  • Finally, slice out the $x$ and $y$ coordinates for the location, and plot them.
In [4]:
# Take the derivative of the state, using Newton's Laws.
# For simplicity, assume that the units have been changed so that:
# - the gravitational constant G = 1
# - the mass of the sun, M = 1
def diff_P(P):
    G = 1.0
    M = 1.0
    x, y, vx, vy = P
    r3 = (sqrt(x**2 + y**2))**3
    ax = -G*M*x/r3
    ay = -G*M*y/r3
    return array([vx, vy, ax, ay])
In [5]:
def plot_euler(steps, h, x=3., y=0., vx=0, vy=.5):
    orbit = empty((steps, 4))
    P = array([x, y, vx, vy])
    for i in xrange(steps):
        orbit[i] = P
        P += h*diff_P(P)
    plot(0, 0, 'yo', ms=6)
    plot(orbit[:, 0], orbit[:, 1], 'r')
    title("Euler's Method, h = {}, {} steps taken".format(h, steps))

Test it out...

In [6]:
plot_euler(100, .1)

The orbit is not connecting back with itself...

In [7]:
plot_euler(500, .1)

Using more steps, it's apparent that the errors in Euler's Method are accumulating

In [8]:
plot_euler(5000, .1)

Using a smaller step size decreases the problem, but doesn't eliminate it

In [9]:
figure(figsize=(10, 5))
subplot('121')
plot_euler(5000, .1)
subplot('122')
plot_euler(50000, .01)

Heun's Method - an "Improved Euler's Method"

Since the changing derivative is the cause of the problem, we try to use a better estimate of the mean derivative between t and t + h.

First calculate the new state using Euler's method:

$\tilde{x} = \vec{x}(t) + h \dfrac{d\vec{x}}{dt}$

Use $\tilde{x}$ to calculate a better estimate of the new state incorporating the changing derivative:

$\vec{x}(t + h) = \vec{x}(t) + \dfrac{h}{2}\left(\dfrac{d\vec{x}}{dt}\bigg|_{\vec{x}(t)} + \dfrac{d\vec{x}}{dt}\bigg|_{\tilde{x}}\right)$

In [10]:
def plot_heun(steps, h, x=3., y=0., vx=0, vy=.5):
    orbit = empty((steps, 4))
    P = array([x, y, vx, vy])
    for i in xrange(steps):
        orbit[i] = P
        Ptilde = P + h*diff_P(P)
        P += h*(diff_P(P) + diff_P(Ptilde))/2.
    plot(0, 0, 'yo', ms=6)
    plot(orbit[:, 0], orbit[:, 1], 'r')
    title("Heun's Method, h = {}, {} steps taken".format(h, steps))

The orbit now connects to itself, eliminating most of the error from Euler's Method.

In [11]:
figure(figsize=(10, 5))
subplot('121')
plot_euler(500, .1)
subplot('122')
plot_heun(500, .1)