Stokes intrigue

gistfile1.py

#!/usr/share/sage-4.8-linux-64bit-ubuntu_10.04.3_lts-x86_64-Linux/sage

from math import cos, sin, exp, pi, e
from sage.all import *

# Terrible ad hoc code follows
x1 = var('x1')
x2 = var('x2')
x3 = var('x3')
x4 = var('x4')
x5 = var('x5')
x6 = var('x6')
Xs = [x1, x2, x3, x4, x5, x6]
# &c. as needed

def identity_matrix(n):
    M = [[0 for j in range(n)] for i in range(n)]
    for i in range(n):
        M[i][i] = 1
    return M

def rotation(*angles):
    n = len(angles) + 1
    Ms = []
    for (i, a) in enumerate(angles):
        M = identity_matrix(n)
        M[i][i] = cos(a)
        M[i][i+1] = sin(a)
        M[i+1][i] = -sin(a)
        M[i+1][i+1] = cos(a)
        Ms.append(matrix(M))
    R = matrix(identity_matrix(n))
    for M in Ms:
        R *= M
    return R

def custom_exterior(R, F):
    n = R.nrows()
    V = ['.' for i in range(n)]
    for i in range(n):
        V[i] = R[i].dot_product(diff(F, Xs[i]))
    return sum(V)

def divergence(F):
    return sum([diff(F[i], eval('x'+str(i+1))) for i in range(len(F))])

def boundary_integral(F, Ss, angles):
    # TODO

def region_integral(F, R, *intervals):
    n = len(intervals)
    #
    # I'm aware that the following is really terrible ad hoc code, but
    # I don't understand how to define a Sage function with an
    # arbitrary number of arguments
    #
    # (e.g., this approach doesn't work, presumably because 'eval'
    # evaluates code rather than treating code-as-code:)
    # args = str(Xs[:n])[1:-1]
    # r(eval(args)) = divergence(F)
    # 
    # (but this is wrong, too:)
    # exec 'r(' + args + ') = divergence(F)'
    if n == 1:
        r(x1) = custom_exterior(R, F)
    elif n == 2:
        r(x1, x2) = custom_exterior(R, F)
    elif n == 3:
        r(x1, x2, x3) = custom_exterior(R, F)
    elif n == 4:
        r(x1, x2, x3, x4) = custom_exterior(R, F)
    elif n == 5:
        r(x1, x2, x3, x4, x5) = custom_exterior(R, F)
    elif n == 5:
        r(x1, x2, x3, x4, x5, x6) = custom_exterior(R, F)
    # &c. as needed
    for (i, I) in enumerate(intervals):        
        r = r.integrate(eval('x'+str(i+1)), I[0], I[1])
    return r


# basic test
F = vector([exp(x1)*sin(x2), exp(x1)*cos(x2), x2*x3**2])
R = rotation(0, 0)

print("custom ext. deriv. ", custom_exterior(R, F))
print("div ", divergence(F))
print("integral of custom. ext ", region_integral(F, R, [0,1], [0,1], [0,2]))