1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
|
#!/usr/bin/env python
"""
This example shows how to work with coordinate transformations, curvilinear
coordinates and a little bit with differential geometry.
It takes polar, cylindrical, spherical, rotating disk coordinates and others
and calculates all kinds of interesting properties, like Jacobian, metric
tensor, Laplace operator, ...
"""
from sympy import var, sin, cos, pprint, Matrix, eye, trigsimp, Eq, \
Function, simplify, sinh, cosh, expand
def laplace(f, g_inv, g_det, X):
"""
Calculates Laplace(f), using the inverse metric g_inv, the determinant of
the metric g_det, all in variables X.
"""
r = 0
for i in range(len(X)):
for j in range(len(X)):
r += g_inv[i, j]*f.diff(X[i]).diff(X[j])
for sigma in range(len(X)):
for alpha in range(len(X)):
r += g_det.diff(X[sigma]) * g_inv[sigma, alpha] * \
f.diff(X[alpha]) / (2*g_det)
return r
def transform(name, X, Y, g_correct=None, recursive=False):
"""
Transforms from cartesian coordinates X to any curvilinear coordinates Y.
It printing useful information, like Jacobian, metric tensor, determinant
of metric, Laplace operator in the new coordinates, ...
g_correct ... if not None, it will be taken as the metric --- this is
useful if sympy's trigsimp() is not powerful enough to
simplify the metric so that it is usable for later
calculation. Leave it as None, only if the metric that
transform() prints is not simplified, you can help it by
specifying the correct one.
recursive ... apply recursive trigonometric simplification (use only when
needed, as it is an expensive operation)
"""
print("_"*80)
print("Transformation:", name)
for x, y in zip(X, Y):
pprint(Eq(y, x))
J = X.jacobian(Y)
print("Jacobian:")
pprint(J)
g = J.T*eye(J.shape[0])*J
g = g.applyfunc(expand)
#g = g.applyfunc(trigsimp)
print("metric tensor g_{ij}:")
pprint(g)
if g_correct is not None:
g = g_correct
print("metric tensor g_{ij} specified by hand:")
pprint(g)
print("inverse metric tensor g^{ij}:")
g_inv = g.inv(method="ADJ")
g_inv = g_inv.applyfunc(simplify)
pprint(g_inv)
print("det g_{ij}:")
g_det = g.det()
pprint(g_det)
f = Function("f")(*list(Y))
print("Laplace:")
pprint(laplace(f, g_inv, g_det, Y))
def main():
var("mu nu rho theta phi sigma tau a t x y z w")
transform("polar", Matrix([rho*cos(phi), rho*sin(phi)]), [rho, phi])
transform("cylindrical", Matrix([rho*cos(phi), rho*sin(phi), z]),
[rho, phi, z])
transform("spherical",
Matrix([rho*sin(theta)*cos(phi), rho*sin(theta)*sin(phi),
rho*cos(theta)]),
[rho, theta, phi],
recursive=True
)
transform("rotating disk",
Matrix([t, x*cos(w*t) - y*sin(w*t), x*sin(w*t) + y*cos(w*t), z]),
[t, x, y, z])
transform("parabolic",
Matrix([sigma*tau, (tau**2 - sigma**2)/2]),
[sigma, tau])
# too complex:
#transform("bipolar",
# Matrix([a*sinh(tau)/(cosh(tau)-cos(sigma)),
# a*sin(sigma)/(cosh(tau)-cos(sigma))]),
# [sigma, tau]
# )
transform("elliptic",
Matrix([a*cosh(mu)*cos(nu), a*sinh(mu)*sin(nu)]),
[mu, nu]
)
if __name__ == "__main__":
main()
|
