The Kitchin Research Group

Matlab post

Recursive functions are functions that call themselves repeatedly until some exit condition is met. Today we look at a classic example of recursive function for computing a factorial. The factorial of a non-negative integer n is denoted n!, and is defined as the product of all positive integers less than or equal to n.

The key ideas in defining a recursive function is that there needs to be some logic to identify when to terminate the function. Then, you need logic that calls the function again, but with a smaller part of the problem. Here we recursively call the function with n-1 until it gets called with n=0. 0! is defined to be 1.

def recursive_factorial(n):
    '''compute the factorial recursively. Note if you put a negative
    number in, this function will never end. We also do not check if
    n is an integer.'''
    if n == 0:
        return 1
    else:
        return n * recursive_factorial(n - 1)

print recursive_factorial(5)

120

from scipy.misc import factorial
print factorial(5)

120.0

n = 5
factorial_loop = 1
for i in range(1, n + 1):
    factorial_loop *= i

print factorial_loop

120

Copyright (C) 2013 by John Kitchin. See the License for information about copying.

org-mode source

Matlab post

Problem statement

find the integral of a function f(x) from a to b i.e.

$$\int_a^b f(x) dx$$

In python we use numerical quadrature to achieve this with the scipy.integrate.quad command.

as a specific example, lets integrate

$$y=x^2$$

from x=0 to x=1. You should be able to work out that the answer is 1/3.

from scipy.integrate import quad

def integrand(x):
    return x**2

ans, err = quad(integrand, 0, 1)
print ans

0.333333333333

from scipy.integrate import dblquad
import numpy as np

def integrand(y, x):
    'y must be the first argument, and x the second.'
    return y * np.sin(x) + x * np.cos(y)

ans, err = dblquad(integrand, np.pi, 2*np.pi,
                   lambda x: 0,
                   lambda x: np.pi)
print ans

-9.86960440109

from scipy.integrate import tplquad
import numpy as np

def integrand(z, y, x):
    return y * np.sin(x) + z * np.cos(x)

ans, err = tplquad(integrand,
                   0, np.pi,  # x limits
                   lambda x: 0,
                   lambda x: 1, # y limits
                   lambda x,y: -1,
                   lambda x,y: 1) # z limits

print ans

2.0

Copyright (C) 2013 by John Kitchin. See the License for information about copying.

org-mode source

Matlab post

The odesolvers in scipy can only solve first order ODEs, or systems of first order ODES. To solve a second order ODE, we must convert it by changes of variables to a system of first order ODES. We consider the Van der Pol oscillator here:

$$\frac{d^2x}{dt^2} - \mu(1-x^2)\frac{dx}{dt} + x = 0$$

\(\mu\) is a constant. If we let \(y=x - x^3/3\) http://en.wikipedia.org/wiki/Van_der_Pol_oscillator, then we arrive at this set of equations:

$$\frac{dx}{dt} = \mu(x-1/3x^3-y)$$

$$\frac{dy}{dt} = \mu/x$$

here is how we solve this set of equations. Let \(\mu=1\).

from scipy.integrate import odeint
import numpy as np

mu = 1.0

def vanderpol(X, t):
    x = X[0]
    y = X[1]
    dxdt = mu * (x - 1./3.*x**3 - y)
    dydt = x / mu
    return [dxdt, dydt]

X0 = [1, 2]
t = np.linspace(0, 40, 250)

sol = odeint(vanderpol, X0, t)

import matplotlib.pyplot as plt
x = sol[:, 0]
y = sol[:, 1]

plt.plot(t,x, t, y)
plt.xlabel('t')
plt.legend(('x', 'y'))
plt.savefig('images/vanderpol-1.png')

# phase portrait
plt.figure()
plt.plot(x,y)
plt.plot(x[0], y[0], 'ro')
plt.xlabel('x')
plt.ylabel('y')
plt.savefig('images/vanderpol-2.png')

Here is the phase portrait. You can see that a limit cycle is approached, indicating periodicity in the solution.

Copyright (C) 2013 by John Kitchin. See the License for information about copying.

org-mode source

Matlab post

Polynomials can be represented as a list of coefficients. For example, the polynomial \(4*x^3 + 3*x^2 -2*x + 10 = 0\) can be represented as [4, 3, -2, 10]. Here are some ways to create a polynomial object, and evaluate it.

import numpy as np

ppar = [4, 3, -2, 10]
p = np.poly1d(ppar)

print p(3)
print np.polyval(ppar, 3)

x = 3
print 4*x**3 + 3*x**2 -2*x + 10

139
139
139

numpy makes it easy to get the derivative and integral of a polynomial.

Consider: \(y = 2x^2 - 1\). We know the derivative is \(4x\). Here we compute the derivative and evaluate it at x=4.

import numpy as np

p = np.poly1d([2, 0, -1])
p2 = np.polyder(p)
print p2
print p2(4)

 
4 x
16

The integral of the previous polynomial is \(\frac{2}{3} x^3 - x + c\). We assume \(C=0\). Let us compute the integral \(\int_2^4 2x^2 - 1 dx\).

import numpy as np

p = np.poly1d([2, 0, -1])
p2 = np.polyint(p)
print p2
print p2(4) - p2(2)

        3
0.6667 x - 1 x
35.3333333333

One reason to use polynomials is the ease of finding all of the roots using numpy.roots.

import numpy as np
print np.roots([2, 0, -1]) # roots are +- sqrt(2)

# note that imaginary roots exist, e.g. x^2 + 1 = 0 has two roots, +-i
p = np.poly1d([1, 0, 1])
print np.roots(p)

[ 0.70710678 -0.70710678]
[ 0.+1.j  0.-1.j]

There are applications of polynomials in thermodynamics. The van der waal equation is a cubic polynomial \(f(V) = V^3 - \frac{p n b + n R T}{p} V^2 + \frac{n^2 a}{p}V - \frac{n^3 a b}{p} = 0\), where \(a\) and \(b\) are constants, \(p\) is the pressure, \(R\) is the gas constant, \(T\) is an absolute temperature and \(n\) is the number of moles. The roots of this equation tell you the volume of the gas at those conditions.

import numpy as np
# numerical values of the constants
a = 3.49e4
b = 1.45
p = 679.7   # pressure in psi
T = 683     # T in Rankine
n = 1.136   # lb-moles
R = 10.73       # ft^3 * psi /R / lb-mol

ppar = [1.0, -(p*n*b+n*R*T)/p, n**2*a/p,  -n**3*a*b/p];
print np.roots(ppar)

[ 5.09432376+0.j          4.40066810+1.43502848j  4.40066810-1.43502848j]

Note that only one root is real (and even then, we have to interpet 0.j as not being imaginary. Also, in a cubic polynomial, there can only be two imaginary roots). In this case that means there is only one phase present.

Copyright (C) 2013 by John Kitchin. See the License for information about copying.

org-mode source