Skip to content

Instantly share code, notes, and snippets.

@zackmdavis
Created May 21, 2012 07:39
Show Gist options
  • Save zackmdavis/2760987 to your computer and use it in GitHub Desktop.
Save zackmdavis/2760987 to your computer and use it in GitHub Desktop.
Stokes intrigue
#!/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]))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment