reduce-historical/r38/packages/excalc/excalc.rlg

Tue Feb 10 12:27:59 2004 run on Linux

*** ^ redefined 
% Problem: Calculate the PDE's for the isovector of the heat equation.
% --------
%         (c.f. B.K. Harrison, f.B. Estabrook, "Geometric Approach...",
%          J. Math. Phys. 12, 653, 1971)

% The heat equation @   psi = @  psi is equivalent to the set of exterior
%                    xx        t

% equations (with u=@ psi, y=@ psi):
%                    T        x


pform {psi,u,x,y,t}=0,a=1,{da,b}=2;



a := d psi - u*d t - y*d x;


a := d psi - d t*u - d x*y


da := - d u^d t - d y^d x;


da := d t^d u + d x^d y


b := u*d x^d t - d y^d t;


b :=  - d t^d x*u + d t^d y



% Now calculate the PDE's for the isovector.

tvector v;



pform {vpsi,vt,vu,vx,vy}=0;


fdomain vpsi=vpsi(psi,t,u,x,y),vt=vt(psi,t,u,x,y),vu=vu(psi,t,u,x,y),
                               vx=vx(psi,t,u,x,y),vy=vy(psi,t,u,x,y);



v := vpsi*@ psi + vt*@ t + vu*@ u + vx*@ x + vy*@ y;


v := @   *vpsi + @ *vt + @ *vu + @ *vx + @ *vy
      psi         t       u       x       y



factor d;


on rat;



i1 := v |_ a - l*a;


i1 := d psi*(@   vpsi - @   vt*u - @   vx*y - l)
              psi        psi        psi

       + d t*(@ vpsi - @ vt*u - @ vx*y + l*u - vu)
               t        t        t

       + d u*(@ vpsi - @ vt*u - @ vx*y)
               u        u        u

       + d x*(@ vpsi - @ vt*u - @ vx*y + l*y - vy)
               x        x        x

       + d y*(@ vpsi - @ vt*u - @ vx*y)
               y        y        y


pform o=1;



o := ot*d t + ox*d x + ou*d u + oy*d y;


o := d t*ot + d u*ou + d x*ox + d y*oy


fdomain f=f(psi,t,u,x,y);



i11 := v _| d a - l*a + d f;


i11 :=  - d psi*l + d t*(l*u - vu) + d u*vt + d x*(l*y - vy) + d y*vx


let vx=-@(f,y),vt=-@(f,u),vu=@(f,t)+u*@(f,psi),vy=@(f,x)+y*@(f,psi),
    vpsi=f-u*@(f,u)-y*@(f,y);



factor ^;



i2 := v |_ b - xi*b - o^a + zeta*da;


i2 := d psi^d t*( - @       f*y - @     f - @     f*u + ot) + d psi^d u*ou
                     psi psi       psi x     psi y

       + d psi^d x*(@     f*u + ox) + d psi^d y*( - @     f + oy)
                     psi u                           psi u

       + d t^d u*(@     f*y + @   f + @   f*u - ou*u + zeta) + d t^d x*(
                   psi u       u x     u y

         @     f*y - @   f*u + @   f*u - @ f + @   f + @   f*u + ot*y - ox*u
          psi x       psi       t u       t     x x     x y

          + u*xi)

       + d t^d y*(@     f*y + @   f - @   f + @   f + @   f*u - oy*u - xi)
                   psi y       psi     t u     x y     y y

       + d u^d x*(@   f*u + ou*y) - d u^d y*@   f
                   u u                       u u

       + d x^d y*( - @   f - @   f*u - oy*y + zeta)
                      u x     u y


let ou=0,oy=@(f,u,psi),ox=-u*@(f,u,psi),
    ot=@(f,x,psi)+u*@(f,y,psi)+y*@(f,psi,psi);



i2;


                                                                   2
d t^d u*(@     f*y + @   f + @   f*u + zeta) + d t^d x*(@       f*y
          psi u       u x     u y                        psi psi

               2
    + @     f*u  + 2*@     f*y + @     f*u*y - @   f*u + @   f*u - @ f + @   f
       psi u          psi x       psi y         psi       t u       t     x x

    + @   f*u + u*xi)
       x y

 + d t^d y*( - @     f*u + @     f*y + @   f - @   f + @   f + @   f*u - xi)
                psi u       psi y       psi     t u     x y     y y

 + d u^d x*@   f*u - d u^d y*@   f
            u u               u u

 + d x^d y*( - @     f*y - @   f - @   f*u + zeta)
                psi u       u x     u y


let zeta=-@(f,u,x)-@(f,u,y)*u-@(f,u,psi)*y;



i2;


                    2            2
d t^d x*(@       f*y  + @     f*u  + 2*@     f*y + @     f*u*y - @   f*u
          psi psi        psi u          psi x       psi y         psi

          + @   f*u - @ f + @   f + @   f*u + u*xi)
             t u       t     x x     x y

 + d t^d y*( - @     f*u + @     f*y + @   f - @   f + @   f + @   f*u - xi)
                psi u       psi y       psi     t u     x y     y y

 + d u^d x*@   f*u - d u^d y*@   f - 2*d x^d y*(@     f*y + @   f + @   f*u)
            u u               u u                psi u       u x     u y


let xi=-@(f,t,u)-u*@(f,u,psi)+@(f,x,y)+u*@(f,y,y)+y*@(f,y,psi)+@(f,psi);



i2;


                    2
d t^d x*(@       f*y  + 2*@     f*y + 2*@     f*u*y - @ f + @   f + 2*@   f*u
          psi psi          psi x         psi y         t     x x       x y

                   2
          + @   f*u ) + d u^d x*@   f*u - d u^d y*@   f
             y y                 u u               u u

 - 2*d x^d y*(@     f*y + @   f + @   f*u)
               psi u       u x     u y


let @(f,u,u)=0;



i2;


                    2
d t^d x*(@       f*y  + 2*@     f*y + 2*@     f*u*y - @ f + @   f + 2*@   f*u
          psi psi          psi x         psi y         t     x x       x y

                   2
          + @   f*u ) - 2*d x^d y*(@     f*y + @   f + @   f*u)
             y y                    psi u       u x     u y
      % These PDE's have to be solved.


clear a,da,b,v,i1,i11,o,i2,xi,t;


remfdomain f,vpsi,vt,vu,vx,vy;


clear @(f,u,u);




% Problem:
% --------
% Calculate the integrability conditions for the system of PDE's:
% (c.f. B.F. Schutz, "Geometrical Methods of Mathematical Physics"
% Cambridge University Press, 1984, p. 156)


% @ z /@ x + a1*z  + b1*z  = c1
%    1           1       2

% @ z /@ y + a2*z  + b2*z  = c2
%    1           1       2

% @ z /@ x + f1*z  + g1*z  = h1
%    2           1       2

% @ z /@ y + f2*z  + g2*z  = h2
%    2           1       2      ;


pform w(k)=1,integ(k)=4,{z(k),x,y}=0,{a,b,c,f,g,h}=1,
      {a1,a2,b1,b2,c1,c2,f1,f2,g1,g2,h1,h2}=0;



fdomain  a1=a1(x,y),a2=a2(x,y),b1=b1(x,y),b2=b2(x,y),
         c1=c1(x,y),c2=c2(x,y),f1=f1(x,y),f2=f2(x,y),
         g1=g1(x,y),g2=g2(x,y),h1=h1(x,y),h2=h2(x,y);




a:=a1*d x+a2*d y$


b:=b1*d x+b2*d y$


c:=c1*d x+c2*d y$


f:=f1*d x+f2*d y$


g:=g1*d x+g2*d y$


h:=h1*d x+h2*d y$



% The equivalent exterior system:
factor d;


w(1) := d z(-1) + z(-1)*a + z(-2)*b - c;


 1
w  := d z  + d x*(z *a1 + z *b1 - c1) + d y*(z *a2 + z *b2 - c2)
         1         1       2                  1       2

w(2) := d z(-2) + z(-1)*f + z(-2)*g - h;


 2
w  := d z  + d x*(z *f1 + z *g1 - h1) + d y*(z *f2 + z *g2 - h2)
         2         1       2                  1       2

indexrange 1,2;


factor z;


% The integrability conditions:

integ(k) := d w(k) ^ w(1) ^ w(2);

     1
integ  := d z ^d z ^d x^d y*z *( - @ a1 + @ a2 + b1*f2 - b2*f1) + 
             1    2          1      y      x

          d z ^d z ^d x^d y*z *( - @ b1 + @ b2 + a1*b2 - a2*b1 + b1*g2 - b2*g1)
             1    2          2      y      x

           + d z ^d z ^d x^d y*(@ c1 - @ c2 - a1*c2 + a2*c1 - b1*h2 + b2*h1)
                1    2           y      x

     2
integ  := d z ^d z ^d x^d y*z *( - @ f1 + @ f2 - a1*f2 + a2*f1 - f1*g2 + f2*g1)
             1    2          1      y      x

           + d z ^d z ^d x^d y*z *( - @ g1 + @ g2 - b1*f2 + b2*f1)
                1    2          2      y      x

           + d z ^d z ^d x^d y*(@ h1 - @ h2 + c1*f2 - c2*f1 - g1*h2 + g2*h1)
                1    2           y      x



clear a,b,c,f,g,h,x,y,w(k),integ(k),z(k);


remfdomain a1,a2,b1,c1,c2,f1,f2,g1,g2,h1,h2;



% Problem:
% --------
% Calculate the PDE's for the generators of the d-theta symmetries of
% the Lagrangian system of the planar Kepler problem.
% c.f. W.Sarlet, F.Cantrijn, Siam Review 23, 467, 1981
% Verify that time translation is a d-theta symmetry and calculate the
% corresponding integral.

pform {t,q(k),v(k),lam(k),tau,xi(k),eta(k)}=0,theta=1,f=0,
      {l,glq(k),glv(k),glt}=0;



tvector gam,y;



indexrange 1,2;



fdomain tau=tau(t,q(k),v(k)),xi=xi(t,q(k),v(k)),f=f(t,q(k),v(k));



l := 1/2*(v(1)**2 + v(2)**2) + m/r$

      % The Lagrangian.

pform r=0;


fdomain r=r(q(k));


let @(r,q 1)=q(1)/r,@(r,q 2)=q(2)/r,q(1)**2+q(2)**2=r**2;



lam(k) := -m*q(k)/r;

             1
   1      - q *m
lam  := ---------
            r

             2
   2      - q *m
lam  := ---------
            r

                                % The force.

gam := @ t + v(k)*@(q(k)) + lam(k)*@(v(k))$



eta(k) := gam _| d xi(k) - v(k)*gam _| d tau$



y  := tau*@ t + xi(k)*@(q(k)) + eta(k)*@(v(k))$

     % Symmetry generator.

theta := l*d t + @(l,v(k))*(d q(k) - v(k)*d t)$



factor @;



s := y |_ theta - d f$



glq(k) := @(q k) _| s;

                                                1   1                1   2
                                        - @  (xi )*q *m      - @  (xi )*q *m
                                            1                    2
   1            1   1         1   2        v                    v
glq  := 2*@  (xi )*v  + @  (xi )*v  + ------------------ + ------------------
            1             2                   r                    r
           q             q

                1          2   2
         + @ (xi ) + @  (xi )*v  - @  f
            t          1             1
                      q             q

                           1 2       2 2
            @  tau*( - 3*(v ) *r - (v ) *r + 2*m)
              1
             q                                               1  2
         + --------------------------------------- - @  tau*v *v
                             2*r                       2
                                                      q

                    1  1               2  1
            @  tau*q *v *m     @  tau*q *v *m
              1                  2
             v                  v                       1
         + ---------------- + ---------------- - @ tau*v
                  r                  r            t

                                                              2   1
                                                      - @  (xi )*q *m
                                                          1
   2          1   1         2   1           2   2        v
glq  := @  (xi )*v  + @  (xi )*v  + 2*@  (xi )*v  + ------------------
          2             1               2                   r
         q             q               q

                     2   2
             - @  (xi )*q *m
                 2
                v                    2                   1  2
         + ------------------ + @ (xi ) - @  f - @  tau*v *v
                   r             t          2      1
                                           q      q

                         1 2         2 2                      1  2
            @  tau*( - (v ) *r - 3*(v ) *r + 2*m)     @  tau*q *v *m
              2                                         1
             q                                         v
         + --------------------------------------- + ----------------
                             2*r                            r

                    2  2
            @  tau*q *v *m
              2
             v                       2
         + ---------------- - @ tau*v
                  r            t


glv(k) := @(v k) _| s;

                                                         1 2       2 2
                                            @  tau*( - (v ) *r - (v ) *r + 2*m)
                                              1
   1          1   1         2   2            v
glv  := @  (xi )*v  + @  (xi )*v  - @  f + -------------------------------------
          1             1             1                     2*r
         v             v             v

                                                         1 2       2 2
                                            @  tau*( - (v ) *r - (v ) *r + 2*m)
                                              2
   2          1   1         2   2            v
glv  := @  (xi )*v  + @  (xi )*v  - @  f + -------------------------------------
          2             2             2                     2*r
         v             v             v


glt := @(t) _| s;


                                                   1   1  1
                                             @  (xi )*q *v *m
                                               1
                1    1 2         1   1  2     v
glt :=  - @  (xi )*(v )  - @  (xi )*v *v  + ------------------
            1                2                      r
           q                q

                 1   2  1
           @  (xi )*q *v *m
             2
            v                        2   1  2         2    2 2
        + ------------------ - @  (xi )*v *v  - @  (xi )*(v )
                  r              1                2
                                q                q

                 2   1  2             2   2  2
           @  (xi )*q *v *m     @  (xi )*q *v *m
             1                    2
            v                    v
        + ------------------ + ------------------ - @ f
                  r                    r             t

                  1    1 2     2 2            2    1 2     2 2
        + @  tau*v *((v )  + (v ) ) + @  tau*v *((v )  + (v ) )
            1                           2
           q                           q

                   1      1 2     2 2              2      1 2     2 2
           @  tau*q *m*((v )  + (v ) )     @  tau*q *m*((v )  + (v ) )
             1                               2
            v                               v
        - ----------------------------- - -----------------------------
                        r                               r

                    1 2       2 2
           @ tau*((v ) *r + (v ) *r + 2*m)         1   1    2   2
            t                                  m*(q *xi  + q *xi )
        + --------------------------------- - ---------------------
                         2*r                            3
                                                       r


% Translation in time must generate a symmetry.
xi(k) := 0;


  k
xi  := 0

tau := 1;


tau := 1


glq k := glq k;

   1
glq  :=  - @  f
             1
            q

   2
glq  :=  - @  f
             2
            q


glv k := glv k;

   1
glv  :=  - @  f
             1
            v

   2
glv  :=  - @  f
             2
            v


glt;


 - @ f
    t


% The corresponding integral is of course the energy.
integ := - y _| theta;


            1 2       2 2
          (v ) *r + (v ) *r - 2*m
integ := -------------------------
                    2*r



clear l,lam k,gam,eta k,y,theta,s,glq k,glv k,glt,t,q k,v k,tau,xi k;


remfdomain r,f,tau,xi;



% Problem:
% --------
% Calculate the "gradient" and "Laplacian" of a function and the "curl"
% and "divergence" of a one-form in elliptic coordinates.


coframe e u = sqrt(cosh(v)**2 - sin(u)**2)*d u,
        e v = sqrt(cosh(v)**2 - sin(u)**2)*d v,
        e phi = cos u*sinh v*d phi;



pform f=0;



fdomain f=f(u,v,phi);



factor e,^;


on rat,gcd;


order cosh v, sin u;


% The gradient:
d f;


    phi                       u                            v
   e   *@   f                e *@ f                       e *@ f
         phi                     u                            v
---------------- + -------------------------- + --------------------------
 cos(u)*sinh(v)                 2         2                  2         2
                    sqrt(cosh(v)  - sin(u) )     sqrt(cosh(v)  - sin(u) )


factor @;


% The Laplacian:
# d # d f;


    @       f               @   f                    @ f*sin(u)
     phi phi                 u u                      u
------------------ + -------------------- - -----------------------------
       2        2            2         2                    2         2
 cos(u) *sinh(v)      cosh(v)  - sin(u)      cos(u)*(cosh(v)  - sin(u) )

          @   f                    @ f*cosh(v)
           v v                      v
 + -------------------- + ------------------------------
           2         2                     2         2
    cosh(v)  - sin(u)      sinh(v)*(cosh(v)  - sin(u) )


% Another way of calculating the Laplacian:
-#vardf(1/2*d f^#d f,f);


    @       f               @   f                    @ f*sin(u)
     phi phi                 u u                      u
------------------ + -------------------- - -----------------------------
       2        2            2         2                    2         2
 cos(u) *sinh(v)      cosh(v)  - sin(u)      cos(u)*(cosh(v)  - sin(u) )

          @   f                    @ f*cosh(v)
           v v                      v
 + -------------------- + ------------------------------
           2         2                     2         2
    cosh(v)  - sin(u)      sinh(v)*(cosh(v)  - sin(u) )


remfac @;



% Now calculate the "curl" and the "divergence" of a one-form.

pform w=1,a(k)=0;



fdomain a=a(u,v,phi);



w := a(-k)*e k;


      phi         u       v
w := e   *a    + e *a  + e *a
           phi       u       v

% The curl:
x := # d w;


       phi            2                 2
x := (e   *( - cosh(v) *@ (a ) + cosh(v) *@ (a ) - cosh(v)*a *sinh(v)
                         v  u              u  v             u

                     2                2
             + sin(u) *@ (a ) - sin(u) *@ (a ) - sin(u)*a *cos(u)))/(
                        v  u             u  v            v

                    2         2          2         2       u
        sqrt(cosh(v)  - sin(u) )*(cosh(v)  - sin(u) )) + (e *(

           cosh(v)*a   *cos(u) + cos(u)*@ (a   )*sinh(v)
                    phi                  v  phi

                          2         2                          2         2
            - sqrt(cosh(v)  - sin(u) )*@   (a )))/(sqrt(cosh(v)  - sin(u) )
                                        phi  v

                             v
        *cos(u)*sinh(v)) + (e *(sin(u)*a   *sinh(v) - cos(u)*@ (a   )*sinh(v)
                                        phi                   u  phi

                          2         2                          2         2
            + sqrt(cosh(v)  - sin(u) )*@   (a )))/(sqrt(cosh(v)  - sin(u) )
                                        phi  u

        *cos(u)*sinh(v))


factor @;


% The divergence:
y := # d # w;


        @   (a   )                @ (a )                       @ (a )
         phi  phi                  u  u                         v  v
y := ---------------- + -------------------------- + -------------------------- 
      cos(u)*sinh(v)                 2         2                  2         2
                         sqrt(cosh(v)  - sin(u) )     sqrt(cosh(v)  - sin(u) )

               3                    2
     + (cosh(v) *a *cos(u) - cosh(v) *sin(u)*a *sinh(v)
                  v                           u

                         2                                      2
         - cosh(v)*sin(u) *a *cos(u) + cosh(v)*a *cos(u)*sinh(v)
                            v                   v

                 3                              2
         + sin(u) *a *sinh(v) - sin(u)*a *cos(u) *sinh(v))/(
                    u                   u

                    2         2                         2         2
        sqrt(cosh(v)  - sin(u) )*cos(u)*sinh(v)*(cosh(v)  - sin(u) ))



remfac @;


clear x,y,w,u,v,phi,e k,a k;


remfdomain a,f;




% Problem:
% --------
% Calculate in a spherical coordinate system the Navier Stokes equations.

coframe e r=d r, e theta =r*d theta, e phi = r*sin theta *d phi;


frame x;



fdomain v=v(t,r,theta,phi),p=p(r,theta,phi);



pform v(k)=0,p=0,w=1;



% We first calculate the convective derivative.

w := v(-k)*e(k)$



factor e;

 on rat;



cdv := @(w,t) + (v(k)*x(-k)) |_ w - 1/2*d(v(k)*v(-k));


         phi
cdv := (e   *(cos(theta)*v   *v      + @   (v   )*v
                          phi  theta    phi  phi   phi

               + @ (v   )*sin(theta)*v *r + @ (v   )*sin(theta)*r
                  r  phi              r      t  phi

               + @     (v   )*sin(theta)*v      + sin(theta)*v   *v ))/(
                  theta  phi              theta               phi  r

                            r
          sin(theta)*r) + (e *(@   (v )*v    + @ (v )*sin(theta)*v *r
                                phi  r   phi    r  r              r

              + @ (v )*sin(theta)*r + @     (v )*sin(theta)*v
                 t  r                  theta  r              theta

                                 2                      2
              - sin(theta)*(v   )  - sin(theta)*(v     ) ))/(sin(theta)*r) + (
                             phi                  theta

           theta                      2
          e     *( - cos(theta)*(v   )  + @   (v     )*v
                                  phi      phi  theta   phi

              + @ (v     )*sin(theta)*v *r + @ (v     )*sin(theta)*r
                 r  theta              r      t  theta

              + @     (v     )*sin(theta)*v      + sin(theta)*v *v     ))/(
                 theta  theta              theta               r  theta

          sin(theta)*r)


%next we calculate the viscous terms;

visc := nu*(d#d# w - #d#d w) + mu*d#d# w;


          phi               2
visc := (e   *( - cos(theta) *v   *nu + cos(theta)*@     (v   )*sin(theta)*nu
                               phi                  theta  phi

                + cos(theta)*@   (v     )*mu + 2*cos(theta)*@   (v     )*nu
                              phi  theta                     phi  theta

                + @       (v   )*mu + @       (v   )*nu
                   phi phi  phi        phi phi  phi

                                       2     2                        2
                + @   (v   )*sin(theta) *nu*r  + 2*@ (v   )*sin(theta) *nu*r
                   r r  phi                         r  phi

                                               2
                + @           (v   )*sin(theta) *nu + @     (v )*sin(theta)*mu*r
                   theta theta  phi                    phi r  r

                + 2*@   (v )*sin(theta)*mu + 2*@   (v )*sin(theta)*nu
                     phi  r                     phi  r

                                                               2
                + @         (v     )*sin(theta)*mu - sin(theta) *v   *nu))/(
                   phi theta  theta                               phi

                     2  2      r
           sin(theta) *r ) + (e *(cos(theta)*@     (v )*sin(theta)*nu
                                              theta  r

               + cos(theta)*@ (v     )*sin(theta)*mu*r
                             r  theta

               - cos(theta)*sin(theta)*v     *mu
                                        theta

               - 2*cos(theta)*sin(theta)*v     *nu
                                          theta

               + @     (v   )*sin(theta)*mu*r - @   (v   )*sin(theta)*mu
                  phi r  phi                     phi  phi

               - 2*@   (v   )*sin(theta)*nu + @       (v )*nu
                    phi  phi                   phi phi  r

                                    2     2                      2     2
               + @   (v )*sin(theta) *mu*r  + @   (v )*sin(theta) *nu*r
                  r r  r                       r r  r

                                    2                           2
               + 2*@ (v )*sin(theta) *mu*r + 2*@ (v )*sin(theta) *nu*r
                    r  r                        r  r

                                            2
               + @           (v )*sin(theta) *nu
                  theta theta  r

                                            2
               + @       (v     )*sin(theta) *mu*r
                  r theta  theta

                                          2                                 2
               - @     (v     )*sin(theta) *mu - 2*@     (v     )*sin(theta) *nu
                  theta  theta                      theta  theta

                             2                     2                    2  2
               - 2*sin(theta) *v *mu - 2*sin(theta) *v *nu))/(sin(theta) *r ) + 
                                r                     r

          theta               2                       2
        (e     *( - cos(theta) *v     *mu - cos(theta) *v     *nu
                                 theta                   theta

                  - cos(theta)*@   (v   )*mu - 2*cos(theta)*@   (v   )*nu
                                phi  phi                     phi  phi

                  + cos(theta)*@     (v     )*sin(theta)*mu
                                theta  theta

                  + cos(theta)*@     (v     )*sin(theta)*nu
                                theta  theta

                  + @         (v   )*sin(theta)*mu
                     phi theta  phi

                                           2                               2
                  + @       (v )*sin(theta) *mu*r + 2*@     (v )*sin(theta) *mu
                     r theta  r                        theta  r

                                           2
                  + 2*@     (v )*sin(theta) *nu + @       (v     )*nu
                       theta  r                    phi phi  theta

                                           2     2
                  + @   (v     )*sin(theta) *nu*r
                     r r  theta

                                           2
                  + 2*@ (v     )*sin(theta) *nu*r
                       r  theta

                                                   2
                  + @           (v     )*sin(theta) *mu
                     theta theta  theta

                                                   2                2
                  + @           (v     )*sin(theta) *nu - sin(theta) *v     *mu
                     theta theta  theta                                theta

                              2                        2  2
                  - sin(theta) *v     *nu))/(sin(theta) *r )
                                 theta


% Finally we add the pressure term and print the components of the
% whole equation.

pform nasteq=1,nast(k)=0;



nasteq := cdv - visc + 1/rho*d p$



factor @;



nast(-k) := x(-k) _| nasteq;

           - @     (v   )*mu     @   (v   )*(mu + 2*nu)      - @       (v )*nu
              phi r  phi          phi  phi                      phi phi  r
nast  := -------------------- + ------------------------ + --------------------
    r        sin(theta)*r                        2                      2  2
                                     sin(theta)*r             sin(theta) *r

             @   (v )*v                             @ (v )*(v *r - 2*mu - 2*nu)
              phi  r   phi                           r  r    r
          + --------------- - @   (v )*(mu + nu) + -----------------------------
             sin(theta)*r      r r  r                            r

                       - @           (v )*nu
                          theta theta  r
          + @ (v ) + ------------------------
             t  r                2
                                r

             @     (v )*( - cos(theta)*nu + sin(theta)*v     *r)
              theta  r                                  theta
          + -----------------------------------------------------
                                            2
                                sin(theta)*r

              - @       (v     )*mu      - @ (v     )*cos(theta)*mu
                 r theta  theta             r  theta
          + ------------------------ + -----------------------------
                       r                       sin(theta)*r

             @     (v     )*(mu + 2*nu)     @ p
              theta  theta                   r
          + ---------------------------- + ----- + (cos(theta)*v     *mu
                          2                 rho                 theta
                         r

                                                         2
             + 2*cos(theta)*v     *nu - sin(theta)*(v   ) *r
                             theta                   phi

                                                                            2
             + 2*sin(theta)*v *mu + 2*sin(theta)*v *nu - sin(theta)*(v     ) *r)
                             r                    r                   theta

                       2
         /(sin(theta)*r )

               - @         (v   )*mu     @   (v   )*cos(theta)*(mu + 2*nu)
                  phi theta  phi          phi  phi
nast      := ------------------------ + -----------------------------------
    theta                     2                             2  2
                  sin(theta)*r                    sin(theta) *r

                  - @       (v )*mu     2*@     (v )*(mu + nu)
                     r theta  r            theta  r
              + -------------------- - ------------------------
                         r                         2
                                                  r

                  - @       (v     )*nu     @   (v     )*v
                     phi phi  theta          phi  theta   phi
              + ------------------------ + ------------------- - @   (v     )*nu
                               2  2           sin(theta)*r        r r  theta
                     sin(theta) *r

                 @ (v     )*(v *r - 2*nu)
                  r  theta    r
              + -------------------------- + @ (v     )
                            r                 t  theta

                 @           (v     )*(mu + nu)
                  theta theta  theta
              - -------------------------------- + (@     (v     )
                                2                    theta  theta
                               r

                *( - cos(theta)*mu - cos(theta)*nu + sin(theta)*v     *r))/(
                                                                 theta

                                  @     p
                            2      theta                2
                sin(theta)*r ) + --------- + (cos(theta) *v     *mu
                                   r*rho                   theta

                             2                                         2
                 + cos(theta) *v     *nu - cos(theta)*sin(theta)*(v   ) *r
                                theta                              phi

                             2                         2
                 + sin(theta) *v *v     *r + sin(theta) *v     *mu
                                r  theta                  theta

                             2                       2  2
                 + sin(theta) *v     *nu)/(sin(theta) *r )
                                theta

               @       (v   )*(mu + nu)     @   (v   )*v
                phi phi  phi                 phi  phi   phi
nast    :=  - -------------------------- + ----------------- - @   (v   )*nu
    phi                       2  2           sin(theta)*r       r r  phi
                    sin(theta) *r

               @ (v   )*(v *r - 2*nu)                 - @           (v   )*nu
                r  phi    r                              theta theta  phi
            + ------------------------ + @ (v   ) + --------------------------
                         r                t  phi                 2
                                                                r

               @     (v   )*( - cos(theta)*nu + sin(theta)*v     *r)
                theta  phi                                  theta
            + -------------------------------------------------------
                                               2
                                   sin(theta)*r

                - @     (v )*mu     2*@   (v )*(mu + nu)
                   phi r  r            phi  r
            + ------------------ - ----------------------
                 sin(theta)*r                      2
                                       sin(theta)*r

                - @         (v     )*mu
                   phi theta  theta
            + --------------------------
                                2
                    sin(theta)*r

               @   (v     )*cos(theta)*( - mu - 2*nu)          @   p
                phi  theta                                      phi
            + ---------------------------------------- + ------------------ + (
                                     2  2                 sin(theta)*r*rho
                           sin(theta) *r

                              2
              v   *(cos(theta) *nu + cos(theta)*sin(theta)*v     *r
               phi                                          theta

                              2                  2                 2  2
                  + sin(theta) *v *r + sin(theta) *nu))/(sin(theta) *r )
                                 r



remfac @,e;



clear v k,x k,nast k,cdv,visc,p,w,nasteq,e k;


remfdomain p,v;




% Problem:
% --------
% Calculate from the Lagrangian of a vibrating rod the equation of
% motion and show that the invariance under time translation leads
% to a conserved current.

pform {y,x,t,q,j}=0,lagr=2;



fdomain y=y(x,t),q=q(x),j=j(x);



factor ^;



lagr := 1/2*(rho*q*@(y,t)**2 - e*j*@(y,x,x)**2)*d x^d t;


                        2              2
         d t^d x*( - @ y *q*rho + @   y *e*j)
                      t            x x
lagr := --------------------------------------
                          2


vardf(lagr,y);


d t^d x*(@   j*@   y*e + 2*@ j*@     y*e + @   y*q*rho + @       y*e*j)
          x x   x x         x   x x x       t t           x x x x


% The Lagrangian does not explicitly depend on time; therefore the
% vector field @ t generates a symmetry. The conserved current is

pform c=1;


factor d;



c := noether(lagr,y,@ t);


c := d t*e*(@ j*@ y*@   y - @   y*@   y*j + @ y*@     y*j)
             x   t   x x     t x   x x       t   x x x

                 2              2
         d x*(@ y *q*rho + @   y *e*j)
               t            x x
      - -------------------------------
                       2


% The exterior derivative of this must be zero or a multiple of the
% equation of motion (weak conservation law) to be a conserved current.

remfac d;



d c;


d t^d x*@ y*( - @   j*@   y*e - 2*@ j*@     y*e - @   y*q*rho - @       y*e*j)
         t       x x   x x         x   x x x       t t           x x x x


% i.e. it is a multiple of the equation of motion.

clear lagr,c,j,y,q;


remfdomain y,q,j;



% Problem:
% --------
% Show that the metric structure given by Eguchi and Hanson induces a
% self-dual curvature.
% c.f. T. Eguchi, P.B. Gilkey, A.J. Hanson, "Gravitation, Gauge Theories
% and Differential Geometry", Physics Reports 66, 213, 1980

for all x let cos(x)**2=1-sin(x)**2;



pform f=0,g=0;


fdomain f=f(r), g=g(r);



coframe   o(r) = f*d r,
      o(theta) = (r/2)*(sin(psi)*d theta - sin(theta)*cos(psi)*d phi),
        o(phi) = (r/2)*(-cos(psi)*d theta - sin(theta)*sin(psi)*d phi),
        o(psi) = (r/2)*g*(d psi + cos(theta)*d phi);



frame e;




pform gamma(a,b)=1,curv2(a,b)=2;


index_symmetries gamma(a,b),curv2(a,b): antisymmetric;



factor o;



gamma(-a,-b) := -(1/2)*( e(-a) _| (e(-c) _| (d o(-b)))
		        -e(-b) _| (e(-a) _| (d o(-c)))
		        +e(-c) _| (e(-b) _| (d o(-a))) )*o(c)$




curv2(-a,b) := d gamma(-a,b) + gamma(-c,b)^gamma(-a,c)$



let f=1/g,g=sqrt(1-(a/r)**4);



pform chck(k,l)=2;


index_symmetries chck(k,l): antisymmetric;



% The following has to be zero for a self-dual curvature.

chck(k,l) := 1/2*eps(k,l,m,n)*curv2(-m,-n) + curv2(k,l);

    r theta
chck        := 0

    r phi
chck      := 0

    theta phi
chck          := 0

    r psi
chck      := 0

    theta psi
chck          := 0

    phi psi
chck        := 0



clear gamma(a,b),curv2(a,b),f,g,chck(a,b),o(k),e(k),r,phi,psi;


remfdomain f,g;



% Example: 6-dimensional FRW model with quadratic curvature terms in
% -------
% the Lagrangian (Lanczos and Gauss-Bonnet terms).
% cf. Henriques, Nuclear Physics, B277, 621 (1986)

for all x let cos(x)**2+sin(x)**2=1;



pform {r,s}=0;


fdomain r=r(t),s=s(t);



coframe o(t) = d t,
        o(1) = r*d u/(1 + k*(u**2)/4),
        o(2) = r*u*d theta/(1 + k*(u**2)/4),
        o(3) = r*u*sin(theta)*d phi/(1 + k*(u**2)/4),
        o(4) = s*d v1,
        o(5) = s*sin(v1)*d v2
 with metric g =-o(t)*o(t)+o(1)*o(1)+o(2)*o(2)+o(3)*o(3)
                +o(4)*o(4)+o(5)*o(5);



frame e;



on nero;

 factor o,^;



riemannconx om;



pform curv(k,l)=2,{riemann(a,b,c,d),ricci(a,b),riccisc}=0;



index_symmetries curv(k,l): antisymmetric,
                 riemann(k,l,m,n): antisymmetric in {k,l},{m,n}
                                   symmetric in {{k,l},{m,n}},
                 ricci(k,l): symmetric;



curv(k,l) := d om(k,l) + om(k,-m)^om(m,l);

             t  1
            o ^o *@   r
    t 1            t t
curv    := -------------
                 r

             t  2
            o ^o *@   r
    t 2            t t
curv    := -------------
                 r

             1  2     2
            o ^o *(@ r  + k)
    1 2             t
curv    := ------------------
                    2
                   r

             t  3
            o ^o *@   r
    t 3            t t
curv    := -------------
                 r

             1  3     2
            o ^o *(@ r  + k)
    1 3             t
curv    := ------------------
                    2
                   r

             2  3     2
            o ^o *(@ r  + k)
    2 3             t
curv    := ------------------
                    2
                   r

             t  4
            o ^o *@   s
    t 4            t t
curv    := -------------
                 s

             1  4
            o ^o *@ r*@ s
    1 4            t   t
curv    := ---------------
                 r*s

             2  4
            o ^o *@ r*@ s
    2 4            t   t
curv    := ---------------
                 r*s

             3  4
            o ^o *@ r*@ s
    3 4            t   t
curv    := ---------------
                 r*s

             t  5
            o ^o *@   s
    t 5            t t
curv    := -------------
                 s

             1  5
            o ^o *@ r*@ s
    1 5            t   t
curv    := ---------------
                 r*s

             2  5
            o ^o *@ r*@ s
    2 5            t   t
curv    := ---------------
                 r*s

             3  5
            o ^o *@ r*@ s
    3 5            t   t
curv    := ---------------
                 r*s

             4  5     2
            o ^o *(@ s  + 1)
    4 5             t
curv    := ------------------
                    2
                   s



riemann(a,b,c,d) := e(d) _| (e (c) _| curv(a,b));

                    - @   r
       t 1 t 1         t t
riemann        := ----------
                      r

                    - @   r
       t 2 t 2         t t
riemann        := ----------
                      r

                      2
                   @ r  + k
       1 2 1 2      t
riemann        := ----------
                       2
                      r

                    - @   r
       t 3 t 3         t t
riemann        := ----------
                      r

                      2
                   @ r  + k
       1 3 1 3      t
riemann        := ----------
                       2
                      r

                      2
                   @ r  + k
       2 3 2 3      t
riemann        := ----------
                       2
                      r

                    - @   s
       t 4 t 4         t t
riemann        := ----------
                      s

                   @ r*@ s
       1 4 1 4      t   t
riemann        := ---------
                     r*s

                   @ r*@ s
       2 4 2 4      t   t
riemann        := ---------
                     r*s

                   @ r*@ s
       3 4 3 4      t   t
riemann        := ---------
                     r*s

                    - @   s
       t 5 t 5         t t
riemann        := ----------
                      s

                   @ r*@ s
       1 5 1 5      t   t
riemann        := ---------
                     r*s

                   @ r*@ s
       2 5 2 5      t   t
riemann        := ---------
                     r*s

                   @ r*@ s
       3 5 3 5      t   t
riemann        := ---------
                     r*s

                      2
                   @ s  + 1
       4 5 4 5      t
riemann        := ----------
                       2
                      s



% The rest is done in the Ricci calculus language,

ricci(-a,-b) := riemann(c,-a,-d,-b)*g(-c,d);

              - 3*@   r*s - 2*@   s*r
                   t t         t t
ricci    := --------------------------
     t t               r*s

                              2
             @   r*r*s + 2*@ r *s + 2*@ r*@ s*r + 2*k*s
              t t           t          t   t
ricci    := --------------------------------------------
     1 1                         2
                                r *s

                              2
             @   r*r*s + 2*@ r *s + 2*@ r*@ s*r + 2*k*s
              t t           t          t   t
ricci    := --------------------------------------------
     2 2                         2
                                r *s

                              2
             @   r*r*s + 2*@ r *s + 2*@ r*@ s*r + 2*k*s
              t t           t          t   t
ricci    := --------------------------------------------
     3 3                         2
                                r *s

                                          2
             3*@ r*@ s*s + @   s*r*s + @ s *r + r
                t   t       t t         t
ricci    := --------------------------------------
     4 4                        2
                             r*s

                                          2
             3*@ r*@ s*s + @   s*r*s + @ s *r + r
                t   t       t t         t
ricci    := --------------------------------------
     5 5                        2
                             r*s



riccisc := ricci(-a,-b)*g(a,b);


                          2        2  2                            2        2  2
riccisc := (2*(3*@   r*r*s  + 3*@ r *s  + 6*@ r*@ s*r*s + 2*@   s*r *s + @ s *r
                  t t            t           t   t           t t          t

                       2    2     2  2
                + 3*k*s  + r ))/(r *s )


pform {laglanc,inv1,inv2} = 0;



index_symmetries riemc3(k,l),riemri(k,l),
                 hlang(k,l),einst(k,l): symmetric;



pform {riemc3(i,j),riemri(i,j)}=0;



riemc3(-i,-j) := riemann(-i,-k,-l,-m)*riemann(-j,k,l,m)$


inv1 := riemc3(-i,-j)*g(i,j);


                   2  2  4        4  4        2    2  2  2        2    4
inv1 := (4*(3*@   r *r *s  + 3*@ r *s  + 6*@ r *@ s *r *s  + 6*@ r *k*s
               t t              t           t    t              t

                      2  4  2      4  4        2  4      2  4    4     4  4
             + 2*@   s *r *s  + @ s *r  + 2*@ s *r  + 3*k *s  + r ))/(r *s )
                  t t            t           t

riemri(-i,-j) := 2*riemann(-i,-k,-j,-l)*ricci(k,l)$


inv2 := ricci(-a,-b)*ricci(a,b);


                   2  2  4              2    4                    2  3
inv2 := (2*(6*@   r *r *s  + 6*@   r*@ r *r*s  + 6*@   r*@ r*@ s*r *s
               t t              t t   t             t t   t   t

                              3  3                4        4  4
             + 6*@   r*@   s*r *s  + 6*@   r*k*r*s  + 6*@ r *s
                  t t   t t             t t              t

                     3        3         2    2  2  2         2    4
             + 12*@ r *@ s*r*s  + 15*@ r *@ s *r *s  + 12*@ r *k*s
                   t    t             t    t               t

                                3  2            3  3                     3
             + 6*@ r*@   s*@ s*r *s  + 6*@ r*@ s *r *s + 12*@ r*@ s*k*r*s
                  t   t t   t             t   t              t   t

                          3            2  4  2              2  4
             + 6*@ r*@ s*r *s + 3*@   s *r *s  + 2*@   s*@ s *r *s
                  t   t            t t              t t   t

                        4        4  4        2  4      2  4    4     4  4
             + 2*@   s*r *s + @ s *r  + 2*@ s *r  + 6*k *s  + r ))/(r *s )
                  t t          t           t

laglanc := (1/2)*(inv1 - 4*inv2 + riccisc**2);


                         2  2                                  2  2            2
laglanc := (12*(@   r*@ r *s  + 4*@   r*@ r*@ s*r*s + @   r*@ s *r  + @   r*k*s
                 t t   t           t t   t   t         t t   t         t t

                          2        3              2                  2    2
                 + @   r*r  + 2*@ r *@ s*s + 2*@ r *@   s*r*s + 3*@ r *@ s *r
                    t t          t    t         t    t t           t    t

                      2                      2
                 + @ r *r + 2*@ r*@   s*@ s*r  + 2*@ r*@ s*k*s + 2*@   s*k*r*s
                    t          t   t t   t          t   t           t t

                      2               3  2
                 + @ s *k*r + k*r))/(r *s )
                    t



pform {einst(a,b),hlang(a,b)}=0;



hlang(-i,-j) := 2*(riemc3(-i,-j) - riemri(-i,-j) -
		   2*ricci(-i,-k)*ricci(-j,K) +
		   riccisc*ricci(-i,-j) - (1/2)*laglanc*g(-i,-j));

hlang    := 
     t t

          3              2    2        2                        2
 12*(2*@ r *@ s*s + 3*@ r *@ s *r + @ r *r + 2*@ r*@ s*k*s + @ s *k*r + k*r)
        t    t         t    t        t          t   t         t
-----------------------------------------------------------------------------
                                     3  2
                                    r *s

                                                  2
hlang    := (4*( - 4*@   r*@ r*@ s*s - 2*@   r*@ s *r - 2*@   r*r
     1 1              t t   t   t         t t   t          t t

                        2                2    2      2
                 - 2*@ r *@   s*s - 3*@ r *@ s  - @ r  - 4*@ r*@   s*@ s*r
                      t    t t         t    t      t        t   t t   t

                                    2           2  2
                 - 2*@   s*k*s - @ s *k - k))/(r *s )
                      t t         t

                                                  2
hlang    := (4*( - 4*@   r*@ r*@ s*s - 2*@   r*@ s *r - 2*@   r*r
     2 2              t t   t   t         t t   t          t t

                        2                2    2      2
                 - 2*@ r *@   s*s - 3*@ r *@ s  - @ r  - 4*@ r*@   s*@ s*r
                      t    t t         t    t      t        t   t t   t

                                    2           2  2
                 - 2*@   s*k*s - @ s *k - k))/(r *s )
                      t t         t

                                                  2
hlang    := (4*( - 4*@   r*@ r*@ s*s - 2*@   r*@ s *r - 2*@   r*r
     3 3              t t   t   t         t t   t          t t

                        2                2    2      2
                 - 2*@ r *@   s*s - 3*@ r *@ s  - @ r  - 4*@ r*@   s*@ s*r
                      t    t t         t    t      t        t   t t   t

                                    2           2  2
                 - 2*@   s*k*s - @ s *k - k))/(r *s )
                      t t         t

                             2                                        3
hlang    := (12*( - @   r*@ r *s - 2*@   r*@ r*@ s*r - @   r*k*s - @ r *@ s
     4 4             t t   t          t t   t   t       t t         t    t

                       2                                     3
                  - @ r *@   s*r - @ r*@ s*k - @   s*k*r))/(r *s)
                     t    t t       t   t       t t

                             2                                        3
hlang    := (12*( - @   r*@ r *s - 2*@   r*@ r*@ s*r - @   r*k*s - @ r *@ s
     5 5             t t   t          t t   t   t       t t         t    t

                       2                                     3
                  - @ r *@   s*r - @ r*@ s*k - @   s*k*r))/(r *s)
                     t    t t       t   t       t t



% The complete Einstein tensor: 

einst(-i,-j) := (ricci(-i,-j) - (1/2)*riccisc*g(-i,-j))*alp1 +
		hlang(-i,-j)*alp2$



alp1 := 1$


factor alp2;



einst(-i,-j) := einst(-i,-j);

                           3              2    2        2
einst    := (12*alp2*(2*@ r *@ s*s + 3*@ r *@ s *r + @ r *r + 2*@ r*@ s*k*s
     t t                 t    t         t    t        t          t   t

                      2               3  2
                 + @ s *k*r + k*r))/(r *s )
                    t

                     2  2                      2  2        2    2
                3*@ r *s  + 6*@ r*@ s*r*s + @ s *r  + 3*k*s  + r
                   t           t   t         t
             + ---------------------------------------------------
                                       2  2
                                      r *s

                                                       2
einst    := (4*alp2*( - 4*@   r*@ r*@ s*s - 2*@   r*@ s *r - 2*@   r*r
     1 1                   t t   t   t         t t   t          t t

                        2                2    2      2
                 - 2*@ r *@   s*s - 3*@ r *@ s  - @ r  - 4*@ r*@   s*@ s*r
                      t    t t         t    t      t        t   t t   t

                                    2           2  2                   2
                 - 2*@   s*k*s - @ s *k - k))/(r *s ) + ( - 2*@   r*r*s
                      t t         t                            t t

                     2  2                            2        2  2      2    2
                - @ r *s  - 4*@ r*@ s*r*s - 2*@   s*r *s - @ s *r  - k*s  - r )/
                   t           t   t           t t          t

              2  2
            (r *s )

                                                       2
einst    := (4*alp2*( - 4*@   r*@ r*@ s*s - 2*@   r*@ s *r - 2*@   r*r
     2 2                   t t   t   t         t t   t          t t

                        2                2    2      2
                 - 2*@ r *@   s*s - 3*@ r *@ s  - @ r  - 4*@ r*@   s*@ s*r
                      t    t t         t    t      t        t   t t   t

                                    2           2  2                   2
                 - 2*@   s*k*s - @ s *k - k))/(r *s ) + ( - 2*@   r*r*s
                      t t         t                            t t

                     2  2                            2        2  2      2    2
                - @ r *s  - 4*@ r*@ s*r*s - 2*@   s*r *s - @ s *r  - k*s  - r )/
                   t           t   t           t t          t

              2  2
            (r *s )

                                                       2
einst    := (4*alp2*( - 4*@   r*@ r*@ s*s - 2*@   r*@ s *r - 2*@   r*r
     3 3                   t t   t   t         t t   t          t t

                        2                2    2      2
                 - 2*@ r *@   s*s - 3*@ r *@ s  - @ r  - 4*@ r*@   s*@ s*r
                      t    t t         t    t      t        t   t t   t

                                    2           2  2                   2
                 - 2*@   s*k*s - @ s *k - k))/(r *s ) + ( - 2*@   r*r*s
                      t t         t                            t t

                     2  2                            2        2  2      2    2
                - @ r *s  - 4*@ r*@ s*r*s - 2*@   s*r *s - @ s *r  - k*s  - r )/
                   t           t   t           t t          t

              2  2
            (r *s )

                                  2                                        3
einst    := (12*alp2*( - @   r*@ r *s - 2*@   r*@ r*@ s*r - @   r*k*s - @ r *@ s
     4 4                  t t   t          t t   t   t       t t         t    t

                      2                                     3
                 - @ r *@   s*r - @ r*@ s*k - @   s*k*r))/(r *s)
                    t    t t       t   t       t t

                                      2                          2
                 - 3*@   r*r*s - 3*@ r *s - 3*@ r*@ s*r - @   s*r  - 3*k*s
                      t t           t          t   t       t t
             + ------------------------------------------------------------
                                            2
                                           r *s

                                  2                                        3
einst    := (12*alp2*( - @   r*@ r *s - 2*@   r*@ r*@ s*r - @   r*k*s - @ r *@ s
     5 5                  t t   t          t t   t   t       t t         t    t

                      2                                     3
                 - @ r *@   s*r - @ r*@ s*k - @   s*k*r))/(r *s)
                    t    t t       t   t       t t

                                      2                          2
                 - 3*@   r*r*s - 3*@ r *s - 3*@ r*@ s*r - @   s*r  - 3*k*s
                      t t           t          t   t       t t
             + ------------------------------------------------------------
                                            2
                                           r *s



clear o(k),e(k),riemc3(i,j),riemri(i,j),curv(k,l),riemann(a,b,c,d),
      ricci(a,b),riccisc,t,u,v1,v2,theta,phi,r,om(k,l),einst(a,b),
      hlang(a,b);



remfdomain r,s;

 

% Problem:
% --------
% Calculate for a given coframe and given torsion the Riemannian part and
% the torsion induced part of the connection. Calculate the curvature.

% For a more elaborate example see E.Schruefer, F.W. Hehl, J.D. McCrea,
% "Application of the REDUCE package EXCALC to the Poincare gauge field
% theory of gravity", GRG Journal, vol. 19, (1988)  197--218

pform {ff, gg}=0;



fdomain ff=ff(r), gg=gg(r);



coframe o(4) = d u + 2*b0*cos(theta)*d phi,
        o(1) = ff*(d u + 2*b0*cos(theta)*d phi) + d r,
        o(2) = gg*d theta,
        o(3) = gg*sin(theta)*d phi
 with metric g = -o(4)*o(1)-o(4)*o(1)+o(2)*o(2)+o(3)*o(3);



frame e;



pform {tor(a),gwt(a)}=2,gamma(a,b)=1,
      {u1,u3,u5}=0;



index_symmetries gamma(a,b): antisymmetric;



fdomain u1=u1(r),u3=u3(r),u5=u5(r);



tor(4) := 0$



tor(1) := -u5*o(4)^o(1) - 2*u3*o(2)^o(3)$



tor(2) := u1*o(4)^o(2) + u3*o(4)^o(3)$



tor(3) := u1*o(4)^o(3) - u3*o(4)^o(2)$



gwt(-a) := d o(-a) - tor(-a)$



% The following is the combined connection.
% The Riemannian part could have equally well been calculated by the
% RIEMANNCONX statement.

gamma(-a,-b) := (1/2)*( e(-b) _| (e(-c) _| gwt(-a))
	               +e(-c) _| (e(-a) _| gwt(-b))
                       -e(-a) _| (e(-b) _| gwt(-c)) )*o(c);

             4
gamma    := o *(@ ff - u5)
     4 1         r

                 2
                o *(@ gg*ff + gg*u1)      3               2
                     r                   o *( - b0*ff + gg *u3)
gamma    :=  - ---------------------- + ------------------------
     4 2                 gg                         2
                                                  gg

              2
             o *@ gg         3
                 r        - o *b0
gamma    := --------- + ----------
     1 2       gg            2
                           gg

                                      3
              2            2         o *(@ gg*ff + gg*u1)
             o *(b0*ff - gg *u3)          r
gamma    := --------------------- - ----------------------
     4 3               2                      gg
                     gg

                        3
              2        o *@ gg
             o *b0         r
gamma    := ------- + ---------
     1 3        2        gg
              gg

              1         3                 4              2
             o *b0     o *cos(theta)     o *(b0*ff - 2*gg *u3)
gamma    := ------- + --------------- + -----------------------
     2 3        2      sin(theta)*gg                2
              gg                                  gg



pform curv(a,b)=2;


index_symmetries curv(a,b): antisymmetric;


factor ^;



curv(-a,b) := d gamma(-a,b) + gamma(-c,b)^gamma(-a,c);

      4        2  3
curv    := (2*o ^o
    4 

                                                 2
            *(@ ff*b0*gg - 2*@ gg*b0*ff + @ gg*gg *u3 - b0*gg*u1 - b0*gg*u5))/
               r              r            r

             3    4  1
           gg  + o ^o *(@   ff - @ u5)
                         r r      r

             1  2           3     2
            o ^o *(@   gg*gg  - b0 )
      4             r r                   4  2                 3               3
curv    := -------------------------- + (o ^o *( - @ ff*@ gg*gg  - @   gg*ff*gg
    2                   4                           r    r          r r
                      gg

                           3        2             2        4
                  + @ gg*gg *u5 - b0 *ff + 2*b0*gg *u3))/gg
                     r

                4  3                                       2
               o ^o *(@ ff*b0*gg - 2*@ gg*b0*ff + 2*@ gg*gg *u3 - b0*gg*u5)
                       r              r              r
            + --------------------------------------------------------------
                                             3
                                           gg

             1  3           3     2
            o ^o *(@   gg*gg  - b0 )
      4             r r
curv    := --------------------------
    3                   4
                      gg

                4  2                                          2
               o ^o *( - @ ff*b0*gg + 2*@ gg*b0*ff - 2*@ gg*gg *u3 + b0*gg*u5)
                          r              r              r
            + ----------------------------------------------------------------- 
                                               3
                                             gg

               4  3                 3               3          3        2
           + (o ^o *( - @ ff*@ gg*gg  - @   gg*ff*gg  + @ gg*gg *u5 - b0 *ff
                         r    r          r r             r

                           2        4
                  + 2*b0*gg *u3))/gg

      1        2  3
curv    := (2*o ^o
    1 

                                                    2
            *( - @ ff*b0*gg + 2*@ gg*b0*ff - @ gg*gg *u3 + b0*gg*u1 + b0*gg*u5))
                  r              r            r

              3    4  1
           /gg  + o ^o *( - @   ff + @ u5)
                             r r      r

      1      1  2                 3               3          3             4
curv    := (o ^o *( - @ ff*@ gg*gg  - @   gg*ff*gg  - @ gg*gg *u1 - @ u1*gg
    2                  r    r          r r             r             r

                    2           2        4     1  3
                - b0 *ff + b0*gg *u3))/gg  + (o ^o *( - @ ff*b0*gg
                                                         r

                                          2             3                3
                  + 2*@ gg*b0*ff + @ gg*gg *u3 + @ u3*gg  + b0*gg*u1))/gg  + (
                       r            r             r

               4  2            4               2   3             3
              o ^o *( - @ ff*gg *u1 + @   gg*ff *gg  + @ gg*ff*gg *u1
                         r             r r              r

                              3                4     2   2             2
                  + @ gg*ff*gg *u5 + @ u1*ff*gg  - b0 *ff  + 3*b0*ff*gg *u3
                     r                r

                      4             4   2     4     4  3         2
                  + gg *u1*u5 - 2*gg *u3 ))/gg  + (o ^o *(@ ff*gg *u3
                                                           r

                                                2
                  - 3*@ gg*ff*gg*u3 - @ u3*ff*gg  + b0*ff*u1 + b0*ff*u5
                       r               r

                        2           2           2
                  - 2*gg *u1*u3 - gg *u3*u5))/gg

      1      1  2
curv    := (o ^o
    3 

                                                 2             3
            *(@ ff*b0*gg - 2*@ gg*b0*ff - @ gg*gg *u3 - @ u3*gg  - b0*gg*u1))/
               r              r            r             r

             3     1  3                 3               3          3
           gg  + (o ^o *( - @ ff*@ gg*gg  - @   gg*ff*gg  - @ gg*gg *u1
                             r    r          r r             r

                           4     2           2        4     4  2            2
                  - @ u1*gg  - b0 *ff + b0*gg *u3))/gg  + (o ^o *( - @ ff*gg *u3
                     r                                                r

                                                2
                  + 3*@ gg*ff*gg*u3 + @ u3*ff*gg  - b0*ff*u1 - b0*ff*u5
                       r               r

                        2           2           2     4  3            4
                  + 2*gg *u1*u3 + gg *u3*u5))/gg  + (o ^o *( - @ ff*gg *u1
                                                                r

                             2   3             3                3
                  + @   gg*ff *gg  + @ gg*ff*gg *u1 + @ gg*ff*gg *u5
                     r r              r                r

                              4     2   2             2        4
                  + @ u1*ff*gg  - b0 *ff  + 3*b0*ff*gg *u3 + gg *u1*u5
                     r

                        4   2     4
                  - 2*gg *u3 ))/gg

      2      1  2                 3               3          3             4
curv    := (o ^o *( - @ ff*@ gg*gg  - @   gg*ff*gg  - @ gg*gg *u1 - @ u1*gg
    4                  r    r          r r             r             r

                    2           2        4     1  3
                - b0 *ff + b0*gg *u3))/gg  + (o ^o *( - @ ff*b0*gg
                                                         r

                                          2             3                3
                  + 2*@ gg*b0*ff + @ gg*gg *u3 + @ u3*gg  + b0*gg*u1))/gg  + (
                       r            r             r

               4  2            4               2   3             3
              o ^o *( - @ ff*gg *u1 + @   gg*ff *gg  + @ gg*ff*gg *u1
                         r             r r              r

                              3                4     2   2             2
                  + @ gg*ff*gg *u5 + @ u1*ff*gg  - b0 *ff  + 3*b0*ff*gg *u3
                     r                r

                      4             4   2     4     4  3         2
                  + gg *u1*u5 - 2*gg *u3 ))/gg  + (o ^o *(@ ff*gg *u3
                                                           r

                                                2
                  - 3*@ gg*ff*gg*u3 - @ u3*ff*gg  + b0*ff*u1 + b0*ff*u5
                       r               r

                        2           2           2
                  - 2*gg *u1*u3 - gg *u3*u5))/gg

             1  2           3     2
            o ^o *(@   gg*gg  - b0 )
      2             r r                   4  2                 3               3
curv    := -------------------------- + (o ^o *( - @ ff*@ gg*gg  - @   gg*ff*gg
    1                   4                           r    r          r r
                      gg

                           3        2             2        4
                  + @ gg*gg *u5 - b0 *ff + 2*b0*gg *u3))/gg
                     r

                4  3                                       2
               o ^o *(@ ff*b0*gg - 2*@ gg*b0*ff + 2*@ gg*gg *u3 - b0*gg*u5)
                       r              r              r
            + --------------------------------------------------------------
                                             3
                                           gg

      2      2  3
curv    := (o ^o
    3 

                       2      2            3          2             2        2
            *( - 2*@ gg *ff*gg  - 2*@ gg*gg *u1 + 6*b0 *ff - 6*b0*gg *u3 + gg ))
                    r                r

                      4  1                                     3
                   2*o ^o *(@ ff*b0*gg - 2*@ gg*b0*ff - @ u3*gg )
              4              r              r            r
           /gg  + ------------------------------------------------
                                          3
                                        gg

      3      1  2
curv    := (o ^o
    4 

                                                 2             3
            *(@ ff*b0*gg - 2*@ gg*b0*ff - @ gg*gg *u3 - @ u3*gg  - b0*gg*u1))/
               r              r            r             r

             3     1  3                 3               3          3
           gg  + (o ^o *( - @ ff*@ gg*gg  - @   gg*ff*gg  - @ gg*gg *u1
                             r    r          r r             r

                           4     2           2        4     4  2            2
                  - @ u1*gg  - b0 *ff + b0*gg *u3))/gg  + (o ^o *( - @ ff*gg *u3
                     r                                                r

                                                2
                  + 3*@ gg*ff*gg*u3 + @ u3*ff*gg  - b0*ff*u1 - b0*ff*u5
                       r               r

                        2           2           2     4  3            4
                  + 2*gg *u1*u3 + gg *u3*u5))/gg  + (o ^o *( - @ ff*gg *u1
                                                                r

                             2   3             3                3
                  + @   gg*ff *gg  + @ gg*ff*gg *u1 + @ gg*ff*gg *u5
                     r r              r                r

                              4     2   2             2        4
                  + @ u1*ff*gg  - b0 *ff  + 3*b0*ff*gg *u3 + gg *u1*u5
                     r

                        4   2     4
                  - 2*gg *u3 ))/gg

             1  3           3     2
            o ^o *(@   gg*gg  - b0 )
      3             r r
curv    := --------------------------
    1                   4
                      gg

                4  2                                          2
               o ^o *( - @ ff*b0*gg + 2*@ gg*b0*ff - 2*@ gg*gg *u3 + b0*gg*u5)
                          r              r              r
            + ----------------------------------------------------------------- 
                                               3
                                             gg

               4  3                 3               3          3        2
           + (o ^o *( - @ ff*@ gg*gg  - @   gg*ff*gg  + @ gg*gg *u5 - b0 *ff
                         r    r          r r             r

                           2        4
                  + 2*b0*gg *u3))/gg

      3      2  3
curv    := (o ^o
    2 

                    2      2            3          2             2        2
            *(2*@ gg *ff*gg  + 2*@ gg*gg *u1 - 6*b0 *ff + 6*b0*gg *u3 - gg ))/
                 r                r

                     4  1                                        3
                  2*o ^o *( - @ ff*b0*gg + 2*@ gg*b0*ff + @ u3*gg )
             4                 r              r            r
           gg  + ---------------------------------------------------
                                           3
                                         gg



clear o(k),e(k),curv(a,b),gamma(a,b),theta,phi,x,y,z,r,s,t,u,v,p,q,c,cs;


remfdomain u1,u3,u5,ff,gg;



showtime;


Time: 1110 ms  plus GC time: 20 ms

end;


Time for test: 1110 ms, plus GC time: 20 ms