Euler Approximations

euler_approximations.py

import math

def f(t, x):
    """An example first-order ODE, in the form x' = f(t, x).

    Arguments:
    t -- the independent variable
    x -- the dependent variable
    """
    return math.exp(t) * math.sin(x)

def steps(t_init, t_final, delta_t):
    """Calculate the number of steps needed to reach t_final.

    Arguments:
    t_init  -- the initial value of the variable
    t_final -- the final value of the variable
    delta_t -- the step size
    """
    return round( (t_final - t_init) / float(delta_t) )

def euler(f, t_init, x_init, t_final, delta_t):
    """
    Compute the array of Euler approximations of the solution of a
    first-order ODE.

    Arguments:
    f       -- a first-order ODE, in the form x' = f(t, x)
    t_init  -- the initial value of the independent variable
    x_init  -- the initial value of the dependent variable
    t_final -- the final value of the independent variable
    delta_t -- the step size
    """
    n = steps(t_init, t_final, delta_t)
    dt = float(delta_t)
    t = float(t_init)
    x = float(x_init)
    approx = [ (t, x) ]
    while len(approx) <= n:
        (t, x) = (t + dt, x + f(t, x) * dt)
        approx.append( (t, x) )
        print( (t, x) ) # added for debugging
    return approx

def runge_kutta(f, t_init, x_init, t_final, delta_t):
    """
    Compute the array of Runge-Kutta approximations of the solution of a
    first-order ODE.

    Arguments:
    f       -- a first-order ODE, in the form x' = f(t, x)
    t_init  -- the initial value of the independent variable
    x_init  -- the initial value of the dependent variable
    t_final -- the final value of the independent variable
    delta_t -- the step size
    """
    n = steps(t_init, t_final, delta_t)
    dt = float(delta_t)
    t = float(t_init)
    x = float(x_init)
    approx = [ (t, x) ]
    while len(approx) <= n:
        (t, x) = (t + dt, runge_kutta_step(f, dt, t, x))
        approx.append( (t, x) )
        print( (t, x) ) # added for debugging
    return approx

def runge_kutta_step(f, t_current, x_current, delta_t):
    """
    Calculate the next value of the independent variable in a
    Runge-Kutta iteration.

    Arguments:
    f         -- a first-order ODE, in the form x' = f(t, x)
    t_current -- the current value of the dependent variable
    x_current -- the current value of the dependent variable
    delta_t   -- the step size
    """
    dt = float(delta_t)
    t = float(t_current)
    x = float(x_current)
    m = f(t, x)
    y = x + m * dt / 2
    n = f(t + dt / 2, y)
    z = x + n * dt / 2
    p = f(t + dt / 2, z)
    w = x + p
    q = f(t + dt, w)
    return x + (m + 2 * n + 2 * p + q) * dt / 6

EulerApproximations.hs

module EulersMethod where

type R = Double
type ODE = (R, R) -> R

next :: ODE -> R -> (R, R) -> (R, R)
next y' dx (x, y) = (x + dx, y + y' (x, y) * dx)

eulers :: ODE -> R -> (R, R) -> [(R, R)]
eulers y' dx (x, y) = (x, y) : eulers y' dx (next y' dx (x, y))

delta :: R -> R -> Int -> R
delta x0 xn n = (xn - x0) / fromIntegral n

euler :: ODE -> (R, R) -> R -> Int -> [(R, R)]
euler y' (x0, y0) xn n = take (n + 1) $ eulers y' (delta x0 xn n) (x0, y0)