Monday, March 16, 2015

SymPy

I've been using Wolfram Alpha for my integration, derivation and algebraic manipulation needs, but I tried SymPy and I have to say that it's pretty neat. I used it to implement the Gram-Schmidt Process to obtain a set of orthonormal basis functions from a set of non-orthogonal, non-normalized and linearly independent polynomial functions and SymPy has done the job beautifully. If you already know a little Python, the learning curve is quite gentle.
from __future__ import division
from sympy import *
from prettytable import PrettyTable
import functools
import string

x=symbols('x')
init_printing(use_unicode=True)

def print_basis(basis):
    """
    pretty print basis functions
    """
    table = PrettyTable(["e%d" %i for i in range(len(basis))])
    table.padding_width = 1
    table.add_row(map(pretty, basis))
    print table

def print_projection_coefficients(coefficients, significant_digits):
    """
    pretty print projection coefficients
    """
    table = PrettyTable(string.ascii_uppercase[:len(coefficients)])
    table.padding_width = 1
    table.add_row(map(pretty, coefficients))
    table.add_row(map(functools.partial(N, n=significant_digits), coefficients))
    print table

def normalize(f, inner_product):
    """
    normalize input function using input inner product operation
    """
    return f/sqrt(inner_product(f, f))

def project(f, basis):
    """
    project input function in space generated by input basis
    """
    projection = 0
    coefficients = []
    for b in basis:
        coefficient = inner_product(f, b)
        coefficients.append(coefficient)
        projection = projection + coefficient * b
    return projection, coefficients

def gram_schmidt(v_basis, inner_product):
    """
    orthonormalize input basis vectors using gram-schmidt algorithm
    """
    e = normalize(v_basis[0], inner_product)
    e_basis = [e]
    for v in v_basis[1:]:
        f = v
        for e in e_basis:
            f = f - inner_product(e, v) * e
        e_basis.append(normalize(f, inner_product))
    return e_basis

# define non orthonormal basis functions here
v_basis = [1, x, x**2, x**3, x**4]

# define inner_product here
inner_product=lambda f1, f2: integrate(f1*f2,(x, -pi, pi))

# orthonormalize basis using gram-schmidt algorithm
e_basis = gram_schmidt(v_basis, inner_product)
print("orthonormalized basis functions:")
print_basis(e_basis)

# get projection of cos(x) in orthonormal basis
significant_digits=6
approximation, coefficients = project(cos(x), e_basis)
print("projection coefficients of cos(x):")
print_projection_coefficients(coefficients, significant_digits)
print("projection of cos(x):")
pprint(approximation.evalf(significant_digits))

SymPy formats the output nicely in unicode. I used PrettyTable to print the output nicely in table form. I'm very pleased with the result. Together with NumPy and Matplotlib, it makes the perfect combo for scientific programming.

orthonormalized basis functions:
+---------+---------+--------------------+------------------------+--------------------------------------+
|    e0   |    e1   |         e2         |           e3           |                  e4                  |
+---------+---------+--------------------+------------------------+--------------------------------------+
|    ___  |   ___   |          ⎛      2⎞ |          ⎛        2  ⎞ |           ⎛          ⎛      2⎞     ⎞ |
|  ╲╱ 2   | ╲╱ 6 ⋅x |     ____ ⎜ 2   π ⎟ |     ____ ⎜ 3   3⋅π ⋅x⎟ |           ⎜        2 ⎜ 2   π ⎟     ⎟ |
| ─────── | ─────── | 3⋅╲╱ 10 ⋅⎜x  - ──⎟ | 5⋅╲╱ 14 ⋅⎜x  - ──────⎟ |           ⎜     6⋅π ⋅⎜x  - ──⎟    4⎟ |
|     ___ |     3/2 |          ⎝     3 ⎠ |          ⎝       5   ⎠ |       ___ ⎜ 4        ⎝     3 ⎠   π ⎟ |
| 2⋅╲╱ π  |  2⋅π    | ────────────────── | ────────────────────── | 105⋅╲╱ 2 ⋅⎜x  - ────────────── - ──⎟ |
|         |         |          5/2       |            7/2         |           ⎝           7          5 ⎠ |
|         |         |       4⋅π          |         4⋅π            | ──────────────────────────────────── |
|         |         |                    |                        |                   9/2                |
|         |         |                    |                        |               16⋅π                   |
+---------+---------+--------------------+------------------------+--------------------------------------+
projection coefficients of cos(x):
+---+---+------------+---+----------------------------------------------------------+
| A | B |     C      | D |                            E                             |
+---+---+------------+---+----------------------------------------------------------+
| 0 | 0 |      ____  | 0 |             ⎛            3⎞             ⎛      3       ⎞ |
|   |   | -3⋅╲╱ 10   |   |         ___ ⎜        16⋅π ⎟         ___ ⎜  16⋅π        ⎟ |
|   |   | ────────── |   |   105⋅╲╱ 2 ⋅⎜-24⋅π + ─────⎟   105⋅╲╱ 2 ⋅⎜- ───── + 24⋅π⎟ |
|   |   |     3/2    |   |             ⎝          7  ⎠             ⎝    7         ⎠ |
|   |   |    π       |   | - ───────────────────────── + ────────────────────────── |
|   |   |            |   |                9/2                         9/2           |
|   |   |            |   |            16⋅π                        16⋅π              |
| 0 | 0 |  -1.70371  | 0 |                         0.486659                         |
+---+---+------------+---+----------------------------------------------------------+
projection of cos(x):
           4             2           
0.0261599⋅x  - 0.452289⋅x  + 0.978327

Saturday, February 7, 2015

Testing D3.Js

Saturday, October 25, 2014

Bifurcation Diagram for the Logistic Map Using D3.JS

$Map$
$x_{n+1} = r \cdot x_n \cdot (1 - x_n)$
$r_{min}$ $r_{max}$ $r$ $x_{min}$ $x_{max}$ $x_0$
$r_{step}$ $size$ $count$ $k$ $n$

Sunday, October 12, 2014

Interactive CobWeb Diagram for the Logistic Map Using D3.JS and DataTables

$Map$
$x_{n+1} = r \cdot x_n \cdot (1 - x_n)$
$r$
$x_0$
$n$

$\text{CobWeb}$
$\text{Time Series}$

Saturday, December 15, 2012

Testing Polls...

If you choose an answer to this question at random, what is the chance you will be correct?

Thursday, December 13, 2012

Testing Pygments...

from __future__ import division #makes float division default

g = 6.672e-11 # m**3 / (Kg * s**2); Gravitational constant

potential_energy = lambda m1, m2, d: -((g*m1*m2)/d)

print 'Potential Energy: %e J' % potential_energy(1, 2, 3)