Sat Jun 29 13:37:28 PDT 1991
REDUCE 3.4, 15-Jul-91 ...
1: 1:
2: 2:
3: 3: showtime;
Time: 0 ms
comment some examples of the FOR statement;
comment summing the squares of the even positive integers
through 50;
for i:=2 step 2 until 50 sum i**2;
22100
comment to set w to the factorial of 10;
w := for i:=1:10 product i;
W := 3628800
comment alternatively, we could set the elements a(i) of the
array a to the factorial of i by the statements;
array a(10);
a(0):=1$
for i:=1:10 do a(i):=i*a(i-1);
comment the above version of the FOR statement does not return
an algebraic value, but we can now use these array
elements as factorials in expressions, e. g.;
1+a(5);
121
comment we could have printed the values of each a(i)
as they were computed by writing the FOR statement as;
for i:=1:10 do write a(i):= i*a(i-1);
A(1) := 1
A(2) := 2
A(3) := 6
A(4) := 24
A(5) := 120
A(6) := 720
A(7) := 5040
A(8) := 40320
A(9) := 362880
A(10) := 3628800
comment another way to use factorials would be to introduce an
operator FAC by an integer procedure as follows;
integer procedure fac (n);
begin integer m;
m:=1;
l1: if n=0 then return m;
m:=m*n;
n:=n-1;
go to l1
end;
FAC
comment we can now use fac as an operator in expressions, e. g.;
z**2+fac(4)-2*fac 2*y;
2
- 4*Y + Z + 24
comment note in the above example that the parentheses around
the arguments of FAC may be omitted since it is a unary operator;
comment the following examples illustrate the solution of some
complete problems;
comment the f and g series (ref Sconzo, P., Leschack, A. R. and
Tobey, R. G., Astronomical Journal, Vol 70 (May 1965);
deps:= -sig*(mu+2*eps)$
dmu:= -3*mu*sig$
dsig:= eps-2*sig**2$
f1:= 1$
g1:= 0$
for i:= 1:8 do
<<f2:= -mu*g1 + deps*df(f1,eps) + dmu*df(f1,mu) + dsig*df(f1,sig);
write "F(",i,") := ",f2;
g2:= f1 + deps*df(g1,eps) + dmu*df(g1,mu) + dsig*df(g1,sig);
write "G(",i,") := ",g2;
f1:=f2;
g1:=g2>>;
F(1) := 0
G(1) := 1
F(2) := - MU
G(2) := 0
F(3) := 3*SIG*MU
G(3) := - MU
2
F(4) := MU*(3*EPS - 15*SIG + MU)
G(4) := 6*SIG*MU
2
F(5) := 15*SIG*MU*( - 3*EPS + 7*SIG - MU)
2
G(5) := MU*(9*EPS - 45*SIG + MU)
2 2 4
F(6) := MU*( - 45*EPS + 630*EPS*SIG - 24*EPS*MU - 945*SIG
2 2
+ 210*SIG *MU - MU )
2
G(6) := 30*SIG*MU*( - 6*EPS + 14*SIG - MU)
2 2 4
F(7) := 63*SIG*MU*(25*EPS - 150*EPS*SIG + 14*EPS*MU + 165*SIG
2 2
- 50*SIG *MU + MU )
2 2 4
G(7) := MU*( - 225*EPS + 3150*EPS*SIG - 54*EPS*MU - 4725*SIG
2 2
+ 630*SIG *MU - MU )
3 2 2 2
F(8) := MU*(1575*EPS - 42525*EPS *SIG + 1107*EPS *MU
4 2 2
+ 155925*EPS*SIG - 24570*EPS*SIG *MU + 117*EPS*MU
6 4 2 2 3
- 135135*SIG + 51975*SIG *MU - 2205*SIG *MU + MU )
2 2 4
G(8) := 126*SIG*MU*(75*EPS - 450*EPS*SIG + 24*EPS*MU + 495*SIG
2 2
- 100*SIG *MU + MU )
comment a problem in Fourier analysis;
factor cos,sin;
on list;
(a1*cos(wt) + a3*cos(3*wt) + b1*sin(wt) + b3*sin(3*wt))**3
where {cos(~x)*cos(~y) => (cos(x+y)+cos(x-y))/2,
cos(~x)*sin(~y) => (sin(x+y)-sin(x-y))/2,
sin(~x)*sin(~y) => (cos(x-y)-cos(x+y))/2,
cos(~x)**2 => (1+cos(2*x))/2,
sin(~x)**2 => (1-cos(2*x))/2};
2
(COS(9*WT)*A3*(A3
2
-3*B3 )
2
+3*COS(7*WT)*(A1*A3
2
-A1*B3
-2*A3*B1*B3)
2
+3*COS(5*WT)*(A1 *A3
2
+A1*A3
-2*A1*B1*B3
2
-A1*B3
2
-A3*B1
+2*A3*B1*B3)
3
+COS(3*WT)*(A1
2
+6*A1 *A3
2
-3*A1*B1
3
+3*A3
2
+6*A3*B1
2
+3*A3*B3 )
3
+3*COS(WT)*(A1
2
+A1 *A3
2
+2*A1*A3
2
+A1*B1
+2*A1*B1*B3
2
+2*A1*B3
2
-A3*B1 )
2
+SIN(9*WT)*B3*(3*A3
2
-B3 )
+3*SIN(7*WT)*(2*A1*A3*B3
2
+A3 *B1
2
-B1*B3 )
2
+3*SIN(5*WT)*(A1 *B3
+2*A1*A3*B1
+2*A1*A3*B3
2
-A3 *B1
2
-B1 *B3
2
+B1*B3 )
2
+SIN(3*WT)*(3*A1 *B1
2
+6*A1 *B3
2
+3*A3 *B3
3
-B1
2
+6*B1 *B3
3
+3*B3 )
2
+3*SIN(WT)*(A1 *B1
2
+A1 *B3
-2*A1*A3*B1
2
+2*A3 *B1
3
+B1
2
-B1 *B3
2
+2*B1*B3 ))/4
remfac cos,sin;
off list;
comment end of Fourier analysis example;
comment the following program, written in collaboration with David
Barton and John Fitch, solves a problem in general relativity. it
will compute the Einstein tensor from any given metric;
on nero;
comment here we introduce the covariant and contravariant metrics;
operator p1,q1,x;
array gg(3,3),h(3,3);
gg(0,0):=e**(q1(x(1)))$
gg(1,1):=-e**(p1(x(1)))$
gg(2,2):=-x(1)**2$
gg(3,3):=-x(1)**2*sin(x(2))**2$
for i:=0:3 do h(i,i):=1/gg(i,i);
comment generate Christoffel symbols and store in arrays cs1 and cs2;
array cs1(3,3,3),cs2(3,3,3);
for i:=0:3 do for j:=i:3 do
<<for k:=0:3 do
cs1(j,i,k) := cs1(i,j,k):=(df(gg(i,k),x(j))+df(gg(j,k),x(i))
-df(gg(i,j),x(k)))/2;
for k:=0:3 do cs2(j,i,k):= cs2(i,j,k) := for p := 0:3
sum h(k,p)*cs1(i,j,p)>>;
comment now compute the Riemann tensor and store in r(i,j,k,l);
array r(3,3,3,3);
for i:=0:3 do for j:=i+1:3 do for k:=i:3 do
for l:=k+1:if k=i then j else 3 do
<<r(j,i,l,k) := r(i,j,k,l) := for q := 0:3
sum gg(i,q)*(df(cs2(k,j,q),x(l))-df(cs2(j,l,q),x(k))
+ for p:=0:3 sum (cs2(p,l,q)*cs2(k,j,p)
-cs2(p,k,q)*cs2(l,j,p)));
r(i,j,l,k) := -r(i,j,k,l);
r(j,i,k,l) := -r(i,j,k,l);
if i neq k or j>l
then <<r(k,l,i,j) := r(l,k,j,i) := r(i,j,k,l);
r(l,k,i,j) := -r(i,j,k,l);
r(k,l,j,i) := -r(i,j,k,l)>>>>;
comment now compute and print the Ricci tensor;
array ricci(3,3);
for i:=0:3 do for j:=0:3 do
write ricci(j,i) := ricci(i,j) := for p := 0:3 sum for q := 0:3
sum h(p,q)*r(q,i,p,j);
Q1(X(1))
RICCI(0,0) := RICCI(0,0) := (E *(
X(1)*DF(P1(X(1)),X(1))*DF(Q1(X(1)),X(1))
2
- 2*X(1)*DF(Q1(X(1)),X(1),2) - X(1)*DF(Q1(X(1)),X(1))
P1(X(1))
- 4*DF(Q1(X(1)),X(1))))/(4*E *X(1))
RICCI(1,1) := RICCI(1,1) := (
- X(1)*DF(P1(X(1)),X(1))*DF(Q1(X(1)),X(1))
2
+ 2*X(1)*DF(Q1(X(1)),X(1),2) + X(1)*DF(Q1(X(1)),X(1))
- 4*DF(P1(X(1)),X(1)))/(4*X(1))
RICCI(2,2) := RICCI(2,2) := ( - X(1)*DF(P1(X(1)),X(1))
P1(X(1)) P1(X(1))
+ X(1)*DF(Q1(X(1)),X(1)) - 2*E + 2)/(2*E )
2
RICCI(3,3) := RICCI(3,3) := (SIN(X(2)) *( - X(1)*DF(P1(X(1)),X(1))
P1(X(1)) P1(X(1))
+ X(1)*DF(Q1(X(1)),X(1)) - 2*E + 2))/(2*E )
comment now compute and print the Ricci scalar;
rs := for i:= 0:3 sum for j:= 0:3 sum h(i,j)*ricci(i,j);
2
RS := (X(1) *DF(P1(X(1)),X(1))*DF(Q1(X(1)),X(1))
2 2 2
- 2*X(1) *DF(Q1(X(1)),X(1),2) - X(1) *DF(Q1(X(1)),X(1))
+ 4*X(1)*DF(P1(X(1)),X(1)) - 4*X(1)*DF(Q1(X(1)),X(1))
P1(X(1)) P1(X(1)) 2
+ 4*E - 4)/(2*E *X(1) )
comment finally compute and print the Einstein tensor;
array einstein(3,3);
for i:=0:3 do for j:=0:3 do
write einstein(i,j):=ricci(i,j)-rs*gg(i,j)/2;
EINSTEIN(0,0) :=
Q1(X(1)) P1(X(1))
E *( - X(1)*DF(P1(X(1)),X(1)) - E + 1)
-------------------------------------------------------
P1(X(1)) 2
E *X(1)
P1(X(1))
- X(1)*DF(Q1(X(1)),X(1)) + E - 1
EINSTEIN(1,1) := -------------------------------------------
2
X(1)
EINSTEIN(2,2) := (X(1)*(X(1)*DF(P1(X(1)),X(1))*DF(Q1(X(1)),X(1))
- 2*X(1)*DF(Q1(X(1)),X(1),2)
2
- X(1)*DF(Q1(X(1)),X(1))
+ 2*DF(P1(X(1)),X(1)) - 2*DF(Q1(X(1)),X(1))))
P1(X(1))
/(4*E )
2
EINSTEIN(3,3) := (X(1)*SIN(X(2)) *(
X(1)*DF(P1(X(1)),X(1))*DF(Q1(X(1)),X(1))
- 2*X(1)*DF(Q1(X(1)),X(1),2)
2
- X(1)*DF(Q1(X(1)),X(1))
+ 2*DF(P1(X(1)),X(1)) - 2*DF(Q1(X(1)),X(1))))
P1(X(1))
/(4*E )
comment end of Einstein tensor program;
clear gg,h,cs1,cs2,r,ricci,einstein;
comment an example using the matrix facility;
matrix xx,yy,zz;
let xx= mat((a11,a12),(a21,a22)),
yy= mat((y1),(y2));
2*det xx - 3*w;
2*(A11*A22 - A12*A21 - 5443200)
zz:= xx**(-1)*yy;
[ - A12*Y2 + A22*Y1 ]
[--------------------]
[ A11*A22 - A12*A21 ]
ZZ := [ ]
[ A11*Y2 - A21*Y1 ]
[------------------- ]
[ A11*A22 - A12*A21 ]
1/xx**2;
2
A12*A21 + A22
MAT((-------------------------------------------,
2 2 2 2
A11 *A22 - 2*A11*A12*A21*A22 + A12 *A21
- A12*(A11 + A22)
-------------------------------------------),
2 2 2 2
A11 *A22 - 2*A11*A12*A21*A22 + A12 *A21
- A21*(A11 + A22)
(-------------------------------------------,
2 2 2 2
A11 *A22 - 2*A11*A12*A21*A22 + A12 *A21
2
A11 + A12*A21
-------------------------------------------))
2 2 2 2
A11 *A22 - 2*A11*A12*A21*A22 + A12 *A21
comment end of matrix examples;
comment a physics example;
on div;
comment this gives us output in same form as Bjorken and Drell;
mass ki= 0, kf= 0, p1= m, pf= m;
vector ei,ef;
mshell ki,kf,p1,pf;
let p1.ei= 0, p1.ef= 0, p1.pf= m**2+ki.kf, p1.ki= m*k,p1.kf=
m*kp, pf.ei= -kf.ei, pf.ef= ki.ef, pf.ki= m*kp, pf.kf=
m*k, ki.ei= 0, ki.kf= m*(k-kp), kf.ef= 0, ei.ei= -1, ef.ef=
-1;
operator gp;
for all p let gp(p)= g(l,p)+m;
comment this is just to save us a lot of writing;
gp(pf)*(g(l,ef,ei,ki)/(2*ki.p1) + g(l,ei,ef,kf)/(2*kf.p1))
* gp(p1)*(g(l,ki,ei,ef)/(2*ki.p1) + g(l,kf,ef,ei)/(2*kf.p1))$
write "The Compton cross-section is ",ws;
2 1 -1 1 -1
The Compton cross-section is 2*EI.EF + ---*K*KP + ---*K *KP - 1
2 2
comment end of first physics example;
off div;
comment another physics example;
index ix,iy,iz;
mass p1=mm,p2=mm,p3= mm,p4= mm,k1=0;
mshell p1,p2,p3,p4,k1;
vector q1,q2;
factor mm,p1.p3;
operator ga,gb;
for all p let ga(p)=g(la,p)+mm, gb(p)= g(lb,p)+mm;
ga(-p2)*g(la,ix)*ga(-p4)*g(la,iy)* (gb(p3)*g(lb,ix)*gb(q1)
*g(lb,iz)*gb(p1)*g(lb,iy)*gb(q2)*g(lb,iz) + gb(p3)
*g(lb,iz)*gb(q2)*g(lb,ix)*gb(p1)*g(lb,iz)*gb(q1)*g(lb,iy))$
let q1=p1-k1, q2=p3+k1;
comment it is usually faster to make such substitutions after all the
trace algebra is done;
write "CXN =",ws;
4 4 2 2
CXN =32*MM *P1.P3 + 8*MM *(P1.K1 - P3.K1) - 16*MM *P1.P3
2
+ 16*MM *P1.P3*( - P1.K1 - P2.P4 + P3.K1)
2
+ 8*MM *( - P1.K1*P2.P4 + P2.P4*P3.K1 - 2*P2.K1*P4.K1) + 8
*P1.P3*(2*P1.P2*P3.P4 + P1.P2*P4.K1 + 2*P1.P4*P2.P3
+ P1.P4*P2.K1 - P2.P3*P4.K1 - P2.K1*P3.P4) + 8*(
- 2*P1.P2*P1.P4*P3.K1 + P1.P2*P1.K1*P3.P4
- P1.P2*P3.P4*P3.K1 + P1.P4*P1.K1*P2.P3 - P1.P4*P2.P3*P3.K1
+ 2*P1.K1*P2.P3*P3.P4)
comment end of second physics example;
showtime;
Time: 7344 ms
end;
4: 4:
Quitting
Sat Jun 29 13:37:36 PDT 1991