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
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
|
/* ***** This package may be broken. Please use the VECT
package on the SHARE directory. - JPG 6/5/78 ***** */
vector:list$
for f in ["grad","div","curl","laplacian","curlgrad","graddiv",
"divcurl","curlcurl"] do prefix(f,112,expr,expr)$
infix("cross",112,112,expr,expr,expr)$
infix("dotdel",108,108,expr,expr,expr)$
nofix("christoffel",expr)$
dimension():=if dimension=3 then [1,2,3] else [1,2]$
type(arg):=if list= /* check required form of answer */
(if listp(arg) then /* operation performed on a vector */
(for element in arg do /* check each argument */
if listp(element) /* an argument is a list */
then return(list)) /* return a list */
else vector) /* operation performed on a scalar */
then result /* return a list */
else apply('matrix,[result]) /* return a matrix */$
_coordsystem&& coordsystem(sys):=
(if sys=rectangular then
(coodvar:[x,y],
scalefactor:[1,1])
else if sys=polar then
(coordvar:[r,th],
scalefactor:[1,r])
else if sys=cartesian then
(coordvar:[x,y,z],
scalefactor:[1,1,1])
else if sys=cylindrical then
(coordvar:[r,ph,z],
scalefactor:[1,r,1])
else if sys=spherical then
(coordvar:[r,th,ph],
scalefactor:[1,r,r*sin(th)])
else (coordvar:read("coordinate variables"),
scalefactor:read("scale factors")),
dimension:length(coordvar),
coordsystem:sys)$
coordsystem(cartesian)$
_cross&& (a cross b) := if dimension=3 then block([result],
result:[a[2]*b[3]-a[3]*b[2],
a[3]*b[1]-a[1]*b[3],
a[1]*b[2]-a[2]*b[1]],
type([a,b]))
else /* 2 dimensional case */
if nonscalarp(a) then
(if nonscalarp(b) then /* vector x vector */
a[1]*b[2]-a[2]*b[1]
else block([result], /* vector x scalar */
result:[a[2]*b,
-a[1]*b],
type([a])))
else block([result], /* scalar x vector */
result:[-a*b[2],
a*b[1]],
type([b]))$
_grad&& (grad s) := block([result],
result:map(lambda([i],
diff(s,coordvar[i])/scalefactor[i]),
dimension()),
type(vector))$
_div&& (div v) := if dimension=3 then
(diff(scalefactor[2]*scalefactor[3]*v[1],coordvar[1])+
diff(scalefactor[3]*scalefactor[1]*v[2],coordvar[2])+
diff(scalefactor[1]*scalefactor[2]*v[3],coordvar[3]))
/scalefactor[1]/scalefactor[2]/scalefactor[3]
else /* 2 dimensional case */
(diff(scalefactor[2]*v[1],coordvar[1])
+diff(scalefactor[1]*v[2],coordvar[2]))
/scalefactor[1]/scalefactor[2]$
_curl&& (curl a) := if dimension=3 then block([result],
result:[(diff(scalefactor[3]*a[3],coordvar[2])
-diff(scalefactor[2]*a[2],coordvar[3]))
/scalefactor[2]/scalefactor[3],
(diff(scalefactor[1]*a[1],coordvar[3])
-diff(scalefactor[3]*a[3],coordvar[1]))
/scalefactor[3]/scalefactor[1],
(diff(scalefactor[2]*a[2],coordvar[1])
-diff(scalefactor[1]*a[1],coordvar[2]))
/scalefactor[1]/scalefactor[2]],
type([a]))
else /* 2 dimensional case */
if nonscalarp(a) then block([result],
result:(diff(scalefactor[2]*a[2],coordvar[1])
-diff(scalefactor[1]*a[1],coordvar[2]))
/scalefactor[1]/scalefactor[2],
type([a]))
else block([result], /* scalar argument */
result:[diff(a,coordvar[2])/scalefactor[2],
-diff(a,coordvar[1])/scalefactor[1]],
type(vector))$
_laplacian&& (laplacian a) := if nonscalarp(a) then grad div a -curl curl a
else if dimension=3 then
(diff(diff(a,coordvar[1])*scalefactor[2]
*scalefactor[3]/scalefactor[1],coordvar[1])
+diff(diff(a,coordvar[2])*scalefactor[3]
*scalefactor[1]/scalefactor[2],coordvar[2])
+diff(diff(a,coordvar[3])*scalefactor[1]
*scalefactor[2]/scalefactor[3],coordvar[3]))
/scalefactor[1]/scalefactor[2]/scalefactor[3]
else /* 2 dimensional case */
(diff(diff(a,coordvar[1])*scalefactor[2]
/scalefactor[1],coordvar[1])
+diff(diff(a,coordvar[2])*scalefactor[1]
/scalefactor[2],coordvar[2]))/scalefactor[1]/scalefactor[2]$
_dotdel&& (v dotdel b) := if nonscalarp(b) then block([result],
result:if last(properties(christsym))=declared\ array
then /* use christoffel symbols */
map(lambda([j],
sum((diff(b[i]*scalefactor[j],
coordvar[i])-sum(b[k]*scalefactor[k]
*christsym[k,j,i],k,1,dimension))
*v[i]/scalefactor[i],i,1,dimension)
/scalefactor[j]),dimension())
else /* vector b, no christoffel symbols */
map(lambda[j],
sum(diff(b[i]*scalefactor[j],coordvar[i])
*v[i]/scalefactor[i],i,1,dimension)
/scalefactor[j],dimension()),
type([v,b]))
else block([result], /* scalar b case */
result:if last(properties(christsym))=declared\ array
then /* use christoffel symbols */
sum((diff(b,coordvar[i])-b*christsym[1,1,i])
*v[i]/scalefactor[i],i,1,dimension)
else /* no christoffel symbols */
sum(diff(b,coordvar[i])*v[i]
/scalefactor[i],i,1,dimension),
type([v]))$
_christoffel&& christoffel := (array(christsym,3,3,3),
christsym[i,j,k]:=0,
for i thru 3 do
(christsym[i,i,i]:diff(scalefactor[i],coordvar[i])
/scalefactor[i],
for j thru 3 do if j#i then
(christsym[i,j,i]:christsym[i,i,j]:diff(scalefactor[i],
coordvar[j])/scalefactor[i],
christsym[j,i,i]:-diff(scalefactor[i],
coordvar[j])*scalefactor[i]/scalefactor[j]^2)))$
(curlgrad s) := 0$
(graddiv v) := block([result],
result:div v,
result:map(lambda([i],
diff(result,coordvar[i])/scalefactor[i]),
dimension()),
type(vector))$
(divcurl v) := 0$
|
