reduce-historical/r34/xlog/reduce.log

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