1% Problem: Calculate the PDE's for the isovector of the heat equation. 2% -------- 3% (c.f. B.K. Harrison, f.B. Estabrook, "Geometric Approach...", 4% J. Math. Phys. 12, 653, 1971) 5 6% The heat equation @ psi = @ psi is equivalent to the set of exterior 7% xx t 8 9% equations (with u=@ psi, y=@ psi): 10% T x 11 12 13pform {psi,u,x,y,t}=0,a=1,{da,b}=2; 14 15 16 17a := d psi - u*d t - y*d x; 18 19 20a := d psi - d t*u - d x*y 21 22 23da := - d u^d t - d y^d x; 24 25 26da := d t^d u + d x^d y 27 28 29b := u*d x^d t - d y^d t; 30 31 32b := - d t^d x*u + d t^d y 33 34 35 36% Now calculate the PDE's for the isovector. 37 38tvector v; 39 40 41 42pform {vpsi,vt,vu,vx,vy}=0; 43 44 45fdomain vpsi=vpsi(psi,t,u,x,y),vt=vt(psi,t,u,x,y),vu=vu(psi,t,u,x,y), 46 vx=vx(psi,t,u,x,y),vy=vy(psi,t,u,x,y); 47 48 49 50v := vpsi*@ psi + vt*@ t + vu*@ u + vx*@ x + vy*@ y; 51 52 53v := @ *vpsi + @ *vt + @ *vu + @ *vx + @ *vy 54 psi t u x y 55 56 57 58factor d; 59 60 61on rat; 62 63 64 65i1 := v |_ a - l*a; 66 67 68i1 := d psi*(@ vpsi - @ vt*u - @ vx*y - l) 69 psi psi psi 70 71 + d t*(@ vpsi - @ vt*u - @ vx*y + l*u - vu) 72 t t t 73 74 + d u*(@ vpsi - @ vt*u - @ vx*y) 75 u u u 76 77 + d x*(@ vpsi - @ vt*u - @ vx*y + l*y - vy) 78 x x x 79 80 + d y*(@ vpsi - @ vt*u - @ vx*y) 81 y y y 82 83 84pform o=1; 85 86 87 88o := ot*d t + ox*d x + ou*d u + oy*d y; 89 90 91o := d t*ot + d u*ou + d x*ox + d y*oy 92 93 94fdomain f=f(psi,t,u,x,y); 95 96 97 98i11 := v _| d a - l*a + d f; 99 100 101i11 := - d psi*l + d t*(l*u - vu) + d u*vt + d x*(l*y - vy) + d y*vx 102 103 104let vx=-@(f,y),vt=-@(f,u),vu=@(f,t)+u*@(f,psi),vy=@(f,x)+y*@(f,psi), 105 vpsi=f-u*@(f,u)-y*@(f,y); 106 107 108 109factor ^; 110 111 112 113i2 := v |_ b - xi*b - o^a + zeta*da; 114 115 116i2 := d psi^d t*( - @ f*y - @ f - @ f*u + ot) + d psi^d u*ou 117 psi psi psi x psi y 118 119 + d psi^d x*(@ f*u + ox) + d psi^d y*( - @ f + oy) 120 psi u psi u 121 122 + d t^d u*(@ f*y + @ f + @ f*u - ou*u + zeta) + d t^d x*( 123 psi u u x u y 124 125 @ f*y - @ f*u + @ f*u - @ f + @ f + @ f*u + ot*y - ox*u 126 psi x psi t u t x x x y 127 128 + u*xi) 129 130 + d t^d y*(@ f*y + @ f - @ f + @ f + @ f*u - oy*u - xi) 131 psi y psi t u x y y y 132 133 + d u^d x*(@ f*u + ou*y) - d u^d y*@ f 134 u u u u 135 136 + d x^d y*( - @ f - @ f*u - oy*y + zeta) 137 u x u y 138 139 140let ou=0,oy=@(f,u,psi),ox=-u*@(f,u,psi), 141 ot=@(f,x,psi)+u*@(f,y,psi)+y*@(f,psi,psi); 142 143 144 145i2; 146 147 148 2 149d t^d u*(@ f*y + @ f + @ f*u + zeta) + d t^d x*(@ f*y 150 psi u u x u y psi psi 151 152 2 153 + @ f*u + 2*@ f*y + @ f*u*y - @ f*u + @ f*u - @ f + @ f 154 psi u psi x psi y psi t u t x x 155 156 + @ f*u + u*xi) 157 x y 158 159 + d t^d y*( - @ f*u + @ f*y + @ f - @ f + @ f + @ f*u - xi) 160 psi u psi y psi t u x y y y 161 162 + d u^d x*@ f*u - d u^d y*@ f 163 u u u u 164 165 + d x^d y*( - @ f*y - @ f - @ f*u + zeta) 166 psi u u x u y 167 168 169let zeta=-@(f,u,x)-@(f,u,y)*u-@(f,u,psi)*y; 170 171 172 173i2; 174 175 176 2 2 177d t^d x*(@ f*y + @ f*u + 2*@ f*y + @ f*u*y - @ f*u 178 psi psi psi u psi x psi y psi 179 180 + @ f*u - @ f + @ f + @ f*u + u*xi) 181 t u t x x x y 182 183 + d t^d y*( - @ f*u + @ f*y + @ f - @ f + @ f + @ f*u - xi) 184 psi u psi y psi t u x y y y 185 186 + d u^d x*@ f*u - d u^d y*@ f - 2*d x^d y*(@ f*y + @ f + @ f*u) 187 u u u u psi u u x u y 188 189 190let xi=-@(f,t,u)-u*@(f,u,psi)+@(f,x,y)+u*@(f,y,y)+y*@(f,y,psi)+@(f,psi); 191 192 193 194i2; 195 196 197 2 198d t^d x*(@ f*y + 2*@ f*y + 2*@ f*u*y - @ f + @ f + 2*@ f*u 199 psi psi psi x psi y t x x x y 200 201 2 202 + @ f*u ) + d u^d x*@ f*u - d u^d y*@ f 203 y y u u u u 204 205 - 2*d x^d y*(@ f*y + @ f + @ f*u) 206 psi u u x u y 207 208 209let @(f,u,u)=0; 210 211 212 213i2; 214 215 216 2 217d t^d x*(@ f*y + 2*@ f*y + 2*@ f*u*y - @ f + @ f + 2*@ f*u 218 psi psi psi x psi y t x x x y 219 220 2 221 + @ f*u ) - 2*d x^d y*(@ f*y + @ f + @ f*u) 222 y y psi u u x u y 223 % These PDE's have to be solved. 224 225 226clear a,da,b,v,i1,i11,o,i2,xi,t; 227 228 229remfdomain f,vpsi,vt,vu,vx,vy; 230 231 232clear @(f,u,u); 233 234 235 236 237% Problem: 238% -------- 239% Calculate the integrability conditions for the system of PDE's: 240% (c.f. B.F. Schutz, "Geometrical Methods of Mathematical Physics" 241% Cambridge University Press, 1984, p. 156) 242 243 244% @ z /@ x + a1*z + b1*z = c1 245% 1 1 2 246 247% @ z /@ y + a2*z + b2*z = c2 248% 1 1 2 249 250% @ z /@ x + f1*z + g1*z = h1 251% 2 1 2 252 253% @ z /@ y + f2*z + g2*z = h2 254% 2 1 2 ; 255 256 257pform w(k)=1,integ(k)=4,{z(k),x,y}=0,{a,b,c,f,g,h}=1, 258 {a1,a2,b1,b2,c1,c2,f1,f2,g1,g2,h1,h2}=0; 259 260 261 262fdomain a1=a1(x,y),a2=a2(x,y),b1=b1(x,y),b2=b2(x,y), 263 c1=c1(x,y),c2=c2(x,y),f1=f1(x,y),f2=f2(x,y), 264 g1=g1(x,y),g2=g2(x,y),h1=h1(x,y),h2=h2(x,y); 265 266 267 268 269a:=a1*d x+a2*d y$ 270 271 272b:=b1*d x+b2*d y$ 273 274 275c:=c1*d x+c2*d y$ 276 277 278f:=f1*d x+f2*d y$ 279 280 281g:=g1*d x+g2*d y$ 282 283 284h:=h1*d x+h2*d y$ 285 286 287 288% The equivalent exterior system: 289factor d; 290 291 292w(1) := d z(-1) + z(-1)*a + z(-2)*b - c; 293 294 295 1 296w := d z + d x*(z *a1 + z *b1 - c1) + d y*(z *a2 + z *b2 - c2) 297 1 1 2 1 2 298 299w(2) := d z(-2) + z(-1)*f + z(-2)*g - h; 300 301 302 2 303w := d z + d x*(z *f1 + z *g1 - h1) + d y*(z *f2 + z *g2 - h2) 304 2 1 2 1 2 305 306indexrange 1,2; 307 308 309factor z; 310 311 312% The integrability conditions: 313 314integ(k) := d w(k) ^ w(1) ^ w(2); 315 316 1 317integ := d z ^d z ^d x^d y*z *( - @ a1 + @ a2 + b1*f2 - b2*f1) + 318 1 2 1 y x 319 320 d z ^d z ^d x^d y*z *( - @ b1 + @ b2 + a1*b2 - a2*b1 + b1*g2 - b2*g1) 321 1 2 2 y x 322 323 + d z ^d z ^d x^d y*(@ c1 - @ c2 - a1*c2 + a2*c1 - b1*h2 + b2*h1) 324 1 2 y x 325 326 2 327integ := d z ^d z ^d x^d y*z *( - @ f1 + @ f2 - a1*f2 + a2*f1 - f1*g2 + f2*g1) 328 1 2 1 y x 329 330 + d z ^d z ^d x^d y*z *( - @ g1 + @ g2 - b1*f2 + b2*f1) 331 1 2 2 y x 332 333 + d z ^d z ^d x^d y*(@ h1 - @ h2 + c1*f2 - c2*f1 - g1*h2 + g2*h1) 334 1 2 y x 335 336 337 338clear a,b,c,f,g,h,x,y,w(k),integ(k),z(k); 339 340 341remfdomain a1,a2,b1,c1,c2,f1,f2,g1,g2,h1,h2; 342 343 344 345% Problem: 346% -------- 347% Calculate the PDE's for the generators of the d-theta symmetries of 348% the Lagrangian system of the planar Kepler problem. 349% c.f. W.Sarlet, F.Cantrijn, Siam Review 23, 467, 1981 350% Verify that time translation is a d-theta symmetry and calculate the 351% corresponding integral. 352 353pform {t,q(k),v(k),lam(k),tau,xi(k),eta(k)}=0,theta=1,f=0, 354 {l,glq(k),glv(k),glt}=0; 355 356 357 358tvector gam,y; 359 360 361 362indexrange 1,2; 363 364 365 366fdomain tau=tau(t,q(k),v(k)),xi=xi(t,q(k),v(k)),f=f(t,q(k),v(k)); 367 368 369 370l := 1/2*(v(1)**2 + v(2)**2) + m/r$ 371 372 % The Lagrangian. 373 374pform r=0; 375 376 377fdomain r=r(q(k)); 378 379 380let @(r,q 1)=q(1)/r,@(r,q 2)=q(2)/r,q(1)**2+q(2)**2=r**2; 381 382 383 384lam(k) := -m*q(k)/r; 385 386 1 387 1 - q *m 388lam := --------- 389 r 390 391 2 392 2 - q *m 393lam := --------- 394 r 395 396 % The force. 397 398gam := @ t + v(k)*@(q(k)) + lam(k)*@(v(k))$ 399 400 401 402eta(k) := gam _| d xi(k) - v(k)*gam _| d tau$ 403 404 405 406y := tau*@ t + xi(k)*@(q(k)) + eta(k)*@(v(k))$ 407 408 % Symmetry generator. 409 410theta := l*d t + @(l,v(k))*(d q(k) - v(k)*d t)$ 411 412 413 414factor @; 415 416 417 418s := y |_ theta - d f$ 419 420 421 422glq(k) := @(q k) _| s; 423 424 1 1 1 2 425 - @ (xi )*q *m - @ (xi )*q *m 426 1 2 427 1 1 1 1 2 v v 428glq := 2*@ (xi )*v + @ (xi )*v + ------------------ + ------------------ 429 1 2 r r 430 q q 431 432 1 2 2 433 + @ (xi ) + @ (xi )*v - @ f 434 t 1 1 435 q q 436 437 1 2 2 2 438 @ tau*( - 3*(v ) *r - (v ) *r + 2*m) 439 1 440 q 1 2 441 + --------------------------------------- - @ tau*v *v 442 2*r 2 443 q 444 445 1 1 2 1 446 @ tau*q *v *m @ tau*q *v *m 447 1 2 448 v v 1 449 + ---------------- + ---------------- - @ tau*v 450 r r t 451 452 2 1 453 - @ (xi )*q *m 454 1 455 2 1 1 2 1 2 2 v 456glq := @ (xi )*v + @ (xi )*v + 2*@ (xi )*v + ------------------ 457 2 1 2 r 458 q q q 459 460 2 2 461 - @ (xi )*q *m 462 2 463 v 2 1 2 464 + ------------------ + @ (xi ) - @ f - @ tau*v *v 465 r t 2 1 466 q q 467 468 1 2 2 2 1 2 469 @ tau*( - (v ) *r - 3*(v ) *r + 2*m) @ tau*q *v *m 470 2 1 471 q v 472 + --------------------------------------- + ---------------- 473 2*r r 474 475 2 2 476 @ tau*q *v *m 477 2 478 v 2 479 + ---------------- - @ tau*v 480 r t 481 482 483glv(k) := @(v k) _| s; 484 485 1 2 2 2 486 @ tau*( - (v ) *r - (v ) *r + 2*m) 487 1 488 1 1 1 2 2 v 489glv := @ (xi )*v + @ (xi )*v - @ f + ------------------------------------- 490 1 1 1 2*r 491 v v v 492 493 1 2 2 2 494 @ tau*( - (v ) *r - (v ) *r + 2*m) 495 2 496 2 1 1 2 2 v 497glv := @ (xi )*v + @ (xi )*v - @ f + ------------------------------------- 498 2 2 2 2*r 499 v v v 500 501 502glt := @(t) _| s; 503 504 505 1 1 1 506 @ (xi )*q *v *m 507 1 508 1 1 2 1 1 2 v 509glt := - @ (xi )*(v ) - @ (xi )*v *v + ------------------ 510 1 2 r 511 q q 512 513 1 2 1 514 @ (xi )*q *v *m 515 2 516 v 2 1 2 2 2 2 517 + ------------------ - @ (xi )*v *v - @ (xi )*(v ) 518 r 1 2 519 q q 520 521 2 1 2 2 2 2 522 @ (xi )*q *v *m @ (xi )*q *v *m 523 1 2 524 v v 525 + ------------------ + ------------------ - @ f 526 r r t 527 528 1 1 2 2 2 2 1 2 2 2 529 + @ tau*v *((v ) + (v ) ) + @ tau*v *((v ) + (v ) ) 530 1 2 531 q q 532 533 1 1 2 2 2 2 1 2 2 2 534 @ tau*q *m*((v ) + (v ) ) @ tau*q *m*((v ) + (v ) ) 535 1 2 536 v v 537 - ----------------------------- - ----------------------------- 538 r r 539 540 1 2 2 2 541 @ tau*((v ) *r + (v ) *r + 2*m) 1 1 2 2 542 t m*(q *xi + q *xi ) 543 + --------------------------------- - --------------------- 544 2*r 3 545 r 546 547 548% Translation in time must generate a symmetry. 549xi(k) := 0; 550 551 552 k 553xi := 0 554 555tau := 1; 556 557 558tau := 1 559 560 561glq k := glq k; 562 563 1 564glq := - @ f 565 1 566 q 567 568 2 569glq := - @ f 570 2 571 q 572 573 574glv k := glv k; 575 576 1 577glv := - @ f 578 1 579 v 580 581 2 582glv := - @ f 583 2 584 v 585 586 587glt; 588 589 590 - @ f 591 t 592 593 594% The corresponding integral is of course the energy. 595integ := - y _| theta; 596 597 598 1 2 2 2 599 (v ) *r + (v ) *r - 2*m 600integ := ------------------------- 601 2*r 602 603 604 605clear l,lam k,gam,eta k,y,theta,s,glq k,glv k,glt,t,q k,v k,tau,xi k; 606 607 608remfdomain r,f,tau,xi; 609 610 611 612% Problem: 613% -------- 614% Calculate the "gradient" and "Laplacian" of a function and the "curl" 615% and "divergence" of a one-form in elliptic coordinates. 616 617 618coframe e u = sqrt(cosh(v)**2 - sin(u)**2)*d u, 619 e v = sqrt(cosh(v)**2 - sin(u)**2)*d v, 620 e phi = cos u*sinh v*d phi; 621 622 623 624pform f=0; 625 626 627 628fdomain f=f(u,v,phi); 629 630 631 632factor e,^; 633 634 635on rat,gcd; 636 637 638order cosh v, sin u; 639 640 641% The gradient: 642d f; 643 644 645 phi u v 646 e *@ f e *@ f e *@ f 647 phi u v 648---------------- + -------------------------- + -------------------------- 649 cos(u)*sinh(v) 2 2 2 2 650 sqrt(cosh(v) - sin(u) ) sqrt(cosh(v) - sin(u) ) 651 652 653factor @; 654 655 656% The Laplacian: 657# d # d f; 658 659 660 @ f @ f @ f*sin(u) 661 phi phi u u u 662------------------ + -------------------- - ----------------------------- 663 2 2 2 2 2 2 664 cos(u) *sinh(v) cosh(v) - sin(u) cos(u)*(cosh(v) - sin(u) ) 665 666 @ f @ f*cosh(v) 667 v v v 668 + -------------------- + ------------------------------ 669 2 2 2 2 670 cosh(v) - sin(u) sinh(v)*(cosh(v) - sin(u) ) 671 672 673% Another way of calculating the Laplacian: 674-#vardf(1/2*d f^#d f,f); 675 676 677 @ f @ f @ f*sin(u) 678 phi phi u u u 679------------------ + -------------------- - ----------------------------- 680 2 2 2 2 2 2 681 cos(u) *sinh(v) cosh(v) - sin(u) cos(u)*(cosh(v) - sin(u) ) 682 683 @ f @ f*cosh(v) 684 v v v 685 + -------------------- + ------------------------------ 686 2 2 2 2 687 cosh(v) - sin(u) sinh(v)*(cosh(v) - sin(u) ) 688 689 690remfac @; 691 692 693 694% Now calculate the "curl" and the "divergence" of a one-form. 695 696pform w=1,a(k)=0; 697 698 699 700fdomain a=a(u,v,phi); 701 702 703 704w := a(-k)*e k; 705 706 707 phi u v 708w := e *a + e *a + e *a 709 phi u v 710 711% The curl: 712x := # d w; 713 714 715 phi 2 2 716x := (e *( - cosh(v) *@ (a ) + cosh(v) *@ (a ) - cosh(v)*a *sinh(v) 717 v u u v u 718 719 2 2 720 + sin(u) *@ (a ) - sin(u) *@ (a ) - sin(u)*a *cos(u)))/( 721 v u u v v 722 723 2 2 2 2 u 724 sqrt(cosh(v) - sin(u) )*(cosh(v) - sin(u) )) + (e *( 725 726 cosh(v)*a *cos(u) + cos(u)*@ (a )*sinh(v) 727 phi v phi 728 729 2 2 2 2 730 - sqrt(cosh(v) - sin(u) )*@ (a )))/(sqrt(cosh(v) - sin(u) ) 731 phi v 732 733 v 734 *cos(u)*sinh(v)) + (e *(sin(u)*a *sinh(v) - cos(u)*@ (a )*sinh(v) 735 phi u phi 736 737 2 2 2 2 738 + sqrt(cosh(v) - sin(u) )*@ (a )))/(sqrt(cosh(v) - sin(u) ) 739 phi u 740 741 *cos(u)*sinh(v)) 742 743 744factor @; 745 746 747% The divergence: 748y := # d # w; 749 750 751 @ (a ) @ (a ) @ (a ) 752 phi phi u u v v 753y := ---------------- + -------------------------- + -------------------------- 754 cos(u)*sinh(v) 2 2 2 2 755 sqrt(cosh(v) - sin(u) ) sqrt(cosh(v) - sin(u) ) 756 757 3 2 758 + (cosh(v) *a *cos(u) - cosh(v) *sin(u)*a *sinh(v) 759 v u 760 761 2 2 762 - cosh(v)*sin(u) *a *cos(u) + cosh(v)*a *cos(u)*sinh(v) 763 v v 764 765 3 2 766 + sin(u) *a *sinh(v) - sin(u)*a *cos(u) *sinh(v))/( 767 u u 768 769 2 2 2 2 770 sqrt(cosh(v) - sin(u) )*cos(u)*sinh(v)*(cosh(v) - sin(u) )) 771 772 773 774remfac @; 775 776 777clear x,y,w,u,v,phi,e k,a k; 778 779 780remfdomain a,f; 781 782 783 784 785% Problem: 786% -------- 787% Calculate in a spherical coordinate system the Navier Stokes equations. 788 789coframe e r=d r, e theta =r*d theta, e phi = r*sin theta *d phi; 790 791 792frame x; 793 794 795 796fdomain v=v(t,r,theta,phi),p=p(r,theta,phi); 797 798 799 800pform v(k)=0,p=0,w=1; 801 802 803 804% We first calculate the convective derivative. 805 806w := v(-k)*e(k)$ 807 808 809 810factor e; 811 812 on rat; 813 814 815 816cdv := @(w,t) + (v(k)*x(-k)) |_ w - 1/2*d(v(k)*v(-k)); 817 818 819 phi 820cdv := (e *(cos(theta)*v *v + @ (v )*v 821 phi theta phi phi phi 822 823 + @ (v )*sin(theta)*v *r + @ (v )*sin(theta)*r 824 r phi r t phi 825 826 + @ (v )*sin(theta)*v + sin(theta)*v *v ))/( 827 theta phi theta phi r 828 829 r 830 sin(theta)*r) + (e *(@ (v )*v + @ (v )*sin(theta)*v *r 831 phi r phi r r r 832 833 + @ (v )*sin(theta)*r + @ (v )*sin(theta)*v 834 t r theta r theta 835 836 2 2 837 - sin(theta)*(v ) - sin(theta)*(v ) ))/(sin(theta)*r) + ( 838 phi theta 839 840 theta 2 841 e *( - cos(theta)*(v ) + @ (v )*v 842 phi phi theta phi 843 844 + @ (v )*sin(theta)*v *r + @ (v )*sin(theta)*r 845 r theta r t theta 846 847 + @ (v )*sin(theta)*v + sin(theta)*v *v ))/( 848 theta theta theta r theta 849 850 sin(theta)*r) 851 852 853%next we calculate the viscous terms; 854 855visc := nu*(d#d# w - #d#d w) + mu*d#d# w; 856 857 858 phi 2 859visc := (e *( - cos(theta) *v *nu + cos(theta)*@ (v )*sin(theta)*nu 860 phi theta phi 861 862 + cos(theta)*@ (v )*mu + 2*cos(theta)*@ (v )*nu 863 phi theta phi theta 864 865 + @ (v )*mu + @ (v )*nu 866 phi phi phi phi phi phi 867 868 2 2 2 869 + @ (v )*sin(theta) *nu*r + 2*@ (v )*sin(theta) *nu*r 870 r r phi r phi 871 872 2 873 + @ (v )*sin(theta) *nu + @ (v )*sin(theta)*mu*r 874 theta theta phi phi r r 875 876 + 2*@ (v )*sin(theta)*mu + 2*@ (v )*sin(theta)*nu 877 phi r phi r 878 879 2 880 + @ (v )*sin(theta)*mu - sin(theta) *v *nu))/( 881 phi theta theta phi 882 883 2 2 r 884 sin(theta) *r ) + (e *(cos(theta)*@ (v )*sin(theta)*nu 885 theta r 886 887 + cos(theta)*@ (v )*sin(theta)*mu*r 888 r theta 889 890 - cos(theta)*sin(theta)*v *mu 891 theta 892 893 - 2*cos(theta)*sin(theta)*v *nu 894 theta 895 896 + @ (v )*sin(theta)*mu*r - @ (v )*sin(theta)*mu 897 phi r phi phi phi 898 899 - 2*@ (v )*sin(theta)*nu + @ (v )*nu 900 phi phi phi phi r 901 902 2 2 2 2 903 + @ (v )*sin(theta) *mu*r + @ (v )*sin(theta) *nu*r 904 r r r r r r 905 906 2 2 907 + 2*@ (v )*sin(theta) *mu*r + 2*@ (v )*sin(theta) *nu*r 908 r r r r 909 910 2 911 + @ (v )*sin(theta) *nu 912 theta theta r 913 914 2 915 + @ (v )*sin(theta) *mu*r 916 r theta theta 917 918 2 2 919 - @ (v )*sin(theta) *mu - 2*@ (v )*sin(theta) *nu 920 theta theta theta theta 921 922 2 2 2 2 923 - 2*sin(theta) *v *mu - 2*sin(theta) *v *nu))/(sin(theta) *r ) + 924 r r 925 926 theta 2 2 927 (e *( - cos(theta) *v *mu - cos(theta) *v *nu 928 theta theta 929 930 - cos(theta)*@ (v )*mu - 2*cos(theta)*@ (v )*nu 931 phi phi phi phi 932 933 + cos(theta)*@ (v )*sin(theta)*mu 934 theta theta 935 936 + cos(theta)*@ (v )*sin(theta)*nu 937 theta theta 938 939 + @ (v )*sin(theta)*mu 940 phi theta phi 941 942 2 2 943 + @ (v )*sin(theta) *mu*r + 2*@ (v )*sin(theta) *mu 944 r theta r theta r 945 946 2 947 + 2*@ (v )*sin(theta) *nu + @ (v )*nu 948 theta r phi phi theta 949 950 2 2 951 + @ (v )*sin(theta) *nu*r 952 r r theta 953 954 2 955 + 2*@ (v )*sin(theta) *nu*r 956 r theta 957 958 2 959 + @ (v )*sin(theta) *mu 960 theta theta theta 961 962 2 2 963 + @ (v )*sin(theta) *nu - sin(theta) *v *mu 964 theta theta theta theta 965 966 2 2 2 967 - sin(theta) *v *nu))/(sin(theta) *r ) 968 theta 969 970 971% Finally we add the pressure term and print the components of the 972% whole equation. 973 974pform nasteq=1,nast(k)=0; 975 976 977 978nasteq := cdv - visc + 1/rho*d p$ 979 980 981 982factor @; 983 984 985 986nast(-k) := x(-k) _| nasteq; 987 988 - @ (v )*mu @ (v )*(mu + 2*nu) - @ (v )*nu 989 phi r phi phi phi phi phi r 990nast := -------------------- + ------------------------ + -------------------- 991 r sin(theta)*r 2 2 2 992 sin(theta)*r sin(theta) *r 993 994 @ (v )*v @ (v )*(v *r - 2*mu - 2*nu) 995 phi r phi r r r 996 + --------------- - @ (v )*(mu + nu) + ----------------------------- 997 sin(theta)*r r r r r 998 999 - @ (v )*nu 1000 theta theta r 1001 + @ (v ) + ------------------------ 1002 t r 2 1003 r 1004 1005 @ (v )*( - cos(theta)*nu + sin(theta)*v *r) 1006 theta r theta 1007 + ----------------------------------------------------- 1008 2 1009 sin(theta)*r 1010 1011 - @ (v )*mu - @ (v )*cos(theta)*mu 1012 r theta theta r theta 1013 + ------------------------ + ----------------------------- 1014 r sin(theta)*r 1015 1016 @ (v )*(mu + 2*nu) @ p 1017 theta theta r 1018 + ---------------------------- + ----- + (cos(theta)*v *mu 1019 2 rho theta 1020 r 1021 1022 2 1023 + 2*cos(theta)*v *nu - sin(theta)*(v ) *r 1024 theta phi 1025 1026 2 1027 + 2*sin(theta)*v *mu + 2*sin(theta)*v *nu - sin(theta)*(v ) *r) 1028 r r theta 1029 1030 2 1031 /(sin(theta)*r ) 1032 1033 - @ (v )*mu @ (v )*cos(theta)*(mu + 2*nu) 1034 phi theta phi phi phi 1035nast := ------------------------ + ----------------------------------- 1036 theta 2 2 2 1037 sin(theta)*r sin(theta) *r 1038 1039 - @ (v )*mu 2*@ (v )*(mu + nu) 1040 r theta r theta r 1041 + -------------------- - ------------------------ 1042 r 2 1043 r 1044 1045 - @ (v )*nu @ (v )*v 1046 phi phi theta phi theta phi 1047 + ------------------------ + ------------------- - @ (v )*nu 1048 2 2 sin(theta)*r r r theta 1049 sin(theta) *r 1050 1051 @ (v )*(v *r - 2*nu) 1052 r theta r 1053 + -------------------------- + @ (v ) 1054 r t theta 1055 1056 @ (v )*(mu + nu) 1057 theta theta theta 1058 - -------------------------------- + (@ (v ) 1059 2 theta theta 1060 r 1061 1062 *( - cos(theta)*mu - cos(theta)*nu + sin(theta)*v *r))/( 1063 theta 1064 1065 @ p 1066 2 theta 2 1067 sin(theta)*r ) + --------- + (cos(theta) *v *mu 1068 r*rho theta 1069 1070 2 2 1071 + cos(theta) *v *nu - cos(theta)*sin(theta)*(v ) *r 1072 theta phi 1073 1074 2 2 1075 + sin(theta) *v *v *r + sin(theta) *v *mu 1076 r theta theta 1077 1078 2 2 2 1079 + sin(theta) *v *nu)/(sin(theta) *r ) 1080 theta 1081 1082 @ (v )*(mu + nu) @ (v )*v 1083 phi phi phi phi phi phi 1084nast := - -------------------------- + ----------------- - @ (v )*nu 1085 phi 2 2 sin(theta)*r r r phi 1086 sin(theta) *r 1087 1088 @ (v )*(v *r - 2*nu) - @ (v )*nu 1089 r phi r theta theta phi 1090 + ------------------------ + @ (v ) + -------------------------- 1091 r t phi 2 1092 r 1093 1094 @ (v )*( - cos(theta)*nu + sin(theta)*v *r) 1095 theta phi theta 1096 + ------------------------------------------------------- 1097 2 1098 sin(theta)*r 1099 1100 - @ (v )*mu 2*@ (v )*(mu + nu) 1101 phi r r phi r 1102 + ------------------ - ---------------------- 1103 sin(theta)*r 2 1104 sin(theta)*r 1105 1106 - @ (v )*mu 1107 phi theta theta 1108 + -------------------------- 1109 2 1110 sin(theta)*r 1111 1112 @ (v )*cos(theta)*( - mu - 2*nu) @ p 1113 phi theta phi 1114 + ---------------------------------------- + ------------------ + ( 1115 2 2 sin(theta)*r*rho 1116 sin(theta) *r 1117 1118 2 1119 v *(cos(theta) *nu + cos(theta)*sin(theta)*v *r 1120 phi theta 1121 1122 2 2 2 2 1123 + sin(theta) *v *r + sin(theta) *nu))/(sin(theta) *r ) 1124 r 1125 1126 1127 1128remfac @,e; 1129 1130 1131 1132clear v k,x k,nast k,cdv,visc,p,w,nasteq,e k; 1133 1134 1135remfdomain p,v; 1136 1137 1138 1139 1140% Problem: 1141% -------- 1142% Calculate from the Lagrangian of a vibrating rod the equation of 1143% motion and show that the invariance under time translation leads 1144% to a conserved current. 1145 1146pform {y,x,t,q,j}=0,lagr=2; 1147 1148 1149 1150fdomain y=y(x,t),q=q(x),j=j(x); 1151 1152 1153 1154factor ^; 1155 1156 1157 1158lagr := 1/2*(rho*q*@(y,t)**2 - e*j*@(y,x,x)**2)*d x^d t; 1159 1160 1161 2 2 1162 d t^d x*( - @ y *q*rho + @ y *e*j) 1163 t x x 1164lagr := -------------------------------------- 1165 2 1166 1167 1168vardf(lagr,y); 1169 1170 1171d t^d x*(@ j*@ y*e + 2*@ j*@ y*e + @ y*q*rho + @ y*e*j) 1172 x x x x x x x x t t x x x x 1173 1174 1175% The Lagrangian does not explicitly depend on time; therefore the 1176% vector field @ t generates a symmetry. The conserved current is 1177 1178pform c=1; 1179 1180 1181factor d; 1182 1183 1184 1185c := noether(lagr,y,@ t); 1186 1187 1188c := d t*e*(@ j*@ y*@ y - @ y*@ y*j + @ y*@ y*j) 1189 x t x x t x x x t x x x 1190 1191 2 2 1192 d x*(@ y *q*rho + @ y *e*j) 1193 t x x 1194 - ------------------------------- 1195 2 1196 1197 1198% The exterior derivative of this must be zero or a multiple of the 1199% equation of motion (weak conservation law) to be a conserved current. 1200 1201remfac d; 1202 1203 1204 1205d c; 1206 1207 1208d t^d x*@ y*( - @ j*@ y*e - 2*@ j*@ y*e - @ y*q*rho - @ y*e*j) 1209 t x x x x x x x x t t x x x x 1210 1211 1212% i.e. it is a multiple of the equation of motion. 1213 1214clear lagr,c,j,y,q; 1215 1216 1217remfdomain y,q,j; 1218 1219 1220 1221% Problem: 1222% -------- 1223% Show that the metric structure given by Eguchi and Hanson induces a 1224% self-dual curvature. 1225% c.f. T. Eguchi, P.B. Gilkey, A.J. Hanson, "Gravitation, Gauge Theories 1226% and Differential Geometry", Physics Reports 66, 213, 1980 1227 1228for all x let cos(x)**2=1-sin(x)**2; 1229 1230 1231 1232pform f=0,g=0; 1233 1234 1235fdomain f=f(r), g=g(r); 1236 1237 1238 1239coframe o(r) = f*d r, 1240 o(theta) = (r/2)*(sin(psi)*d theta - sin(theta)*cos(psi)*d phi), 1241 o(phi) = (r/2)*(-cos(psi)*d theta - sin(theta)*sin(psi)*d phi), 1242 o(psi) = (r/2)*g*(d psi + cos(theta)*d phi); 1243 1244 1245 1246frame e; 1247 1248 1249 1250 1251pform gamma(a,b)=1,curv2(a,b)=2; 1252 1253 1254index_symmetries gamma(a,b),curv2(a,b): antisymmetric; 1255 1256 1257 1258factor o; 1259 1260 1261 1262gamma(-a,-b) := -(1/2)*( e(-a) _| (e(-c) _| (d o(-b))) 1263 -e(-b) _| (e(-a) _| (d o(-c))) 1264 +e(-c) _| (e(-b) _| (d o(-a))) )*o(c)$ 1265 1266 1267 1268 1269curv2(-a,b) := d gamma(-a,b) + gamma(-c,b)^gamma(-a,c)$ 1270 1271 1272 1273let f=1/g,g=sqrt(1-(a/r)**4); 1274 1275 1276 1277pform chck(k,l)=2; 1278 1279 1280index_symmetries chck(k,l): antisymmetric; 1281 1282 1283 1284% The following has to be zero for a self-dual curvature. 1285 1286chck(k,l) := 1/2*eps(k,l,m,n)*curv2(-m,-n) + curv2(k,l); 1287 1288 r theta 1289chck := 0 1290 1291 r phi 1292chck := 0 1293 1294 theta phi 1295chck := 0 1296 1297 r psi 1298chck := 0 1299 1300 theta psi 1301chck := 0 1302 1303 phi psi 1304chck := 0 1305 1306 1307 1308clear gamma(a,b),curv2(a,b),f,g,chck(a,b),o(k),e(k),r,phi,psi; 1309 1310 1311remfdomain f,g; 1312 1313 1314 1315% Example: 6-dimensional FRW model with quadratic curvature terms in 1316% ------- 1317% the Lagrangian (Lanczos and Gauss-Bonnet terms). 1318% cf. Henriques, Nuclear Physics, B277, 621 (1986) 1319 1320for all x let cos(x)**2+sin(x)**2=1; 1321 1322 1323 1324pform {r,s}=0; 1325 1326 1327fdomain r=r(t),s=s(t); 1328 1329 1330 1331coframe o(t) = d t, 1332 o(1) = r*d u/(1 + k*(u**2)/4), 1333 o(2) = r*u*d theta/(1 + k*(u**2)/4), 1334 o(3) = r*u*sin(theta)*d phi/(1 + k*(u**2)/4), 1335 o(4) = s*d v1, 1336 o(5) = s*sin(v1)*d v2 1337 with metric g =-o(t)*o(t)+o(1)*o(1)+o(2)*o(2)+o(3)*o(3) 1338 +o(4)*o(4)+o(5)*o(5); 1339 1340 1341 1342frame e; 1343 1344 1345 1346on nero; 1347 1348 factor o,^; 1349 1350 1351 1352riemannconx om; 1353 1354 1355 1356pform curv(k,l)=2,{riemann(a,b,c,d),ricci(a,b),riccisc}=0; 1357 1358 1359 1360index_symmetries curv(k,l): antisymmetric, 1361 riemann(k,l,m,n): antisymmetric in {k,l},{m,n} 1362 symmetric in {{k,l},{m,n}}, 1363 ricci(k,l): symmetric; 1364 1365 1366 1367curv(k,l) := d om(k,l) + om(k,-m)^om(m,l); 1368 1369 t 1 1370 o ^o *@ r 1371 t 1 t t 1372curv := ------------- 1373 r 1374 1375 t 2 1376 o ^o *@ r 1377 t 2 t t 1378curv := ------------- 1379 r 1380 1381 1 2 2 1382 o ^o *(@ r + k) 1383 1 2 t 1384curv := ------------------ 1385 2 1386 r 1387 1388 t 3 1389 o ^o *@ r 1390 t 3 t t 1391curv := ------------- 1392 r 1393 1394 1 3 2 1395 o ^o *(@ r + k) 1396 1 3 t 1397curv := ------------------ 1398 2 1399 r 1400 1401 2 3 2 1402 o ^o *(@ r + k) 1403 2 3 t 1404curv := ------------------ 1405 2 1406 r 1407 1408 t 4 1409 o ^o *@ s 1410 t 4 t t 1411curv := ------------- 1412 s 1413 1414 1 4 1415 o ^o *@ r*@ s 1416 1 4 t t 1417curv := --------------- 1418 r*s 1419 1420 2 4 1421 o ^o *@ r*@ s 1422 2 4 t t 1423curv := --------------- 1424 r*s 1425 1426 3 4 1427 o ^o *@ r*@ s 1428 3 4 t t 1429curv := --------------- 1430 r*s 1431 1432 t 5 1433 o ^o *@ s 1434 t 5 t t 1435curv := ------------- 1436 s 1437 1438 1 5 1439 o ^o *@ r*@ s 1440 1 5 t t 1441curv := --------------- 1442 r*s 1443 1444 2 5 1445 o ^o *@ r*@ s 1446 2 5 t t 1447curv := --------------- 1448 r*s 1449 1450 3 5 1451 o ^o *@ r*@ s 1452 3 5 t t 1453curv := --------------- 1454 r*s 1455 1456 4 5 2 1457 o ^o *(@ s + 1) 1458 4 5 t 1459curv := ------------------ 1460 2 1461 s 1462 1463 1464 1465riemann(a,b,c,d) := e(d) _| (e (c) _| curv(a,b)); 1466 1467 - @ r 1468 t 1 t 1 t t 1469riemann := ---------- 1470 r 1471 1472 - @ r 1473 t 2 t 2 t t 1474riemann := ---------- 1475 r 1476 1477 2 1478 @ r + k 1479 1 2 1 2 t 1480riemann := ---------- 1481 2 1482 r 1483 1484 - @ r 1485 t 3 t 3 t t 1486riemann := ---------- 1487 r 1488 1489 2 1490 @ r + k 1491 1 3 1 3 t 1492riemann := ---------- 1493 2 1494 r 1495 1496 2 1497 @ r + k 1498 2 3 2 3 t 1499riemann := ---------- 1500 2 1501 r 1502 1503 - @ s 1504 t 4 t 4 t t 1505riemann := ---------- 1506 s 1507 1508 @ r*@ s 1509 1 4 1 4 t t 1510riemann := --------- 1511 r*s 1512 1513 @ r*@ s 1514 2 4 2 4 t t 1515riemann := --------- 1516 r*s 1517 1518 @ r*@ s 1519 3 4 3 4 t t 1520riemann := --------- 1521 r*s 1522 1523 - @ s 1524 t 5 t 5 t t 1525riemann := ---------- 1526 s 1527 1528 @ r*@ s 1529 1 5 1 5 t t 1530riemann := --------- 1531 r*s 1532 1533 @ r*@ s 1534 2 5 2 5 t t 1535riemann := --------- 1536 r*s 1537 1538 @ r*@ s 1539 3 5 3 5 t t 1540riemann := --------- 1541 r*s 1542 1543 2 1544 @ s + 1 1545 4 5 4 5 t 1546riemann := ---------- 1547 2 1548 s 1549 1550 1551 1552% The rest is done in the Ricci calculus language, 1553 1554ricci(-a,-b) := riemann(c,-a,-d,-b)*g(-c,d); 1555 1556 - 3*@ r*s - 2*@ s*r 1557 t t t t 1558ricci := -------------------------- 1559 t t r*s 1560 1561 2 1562 @ r*r*s + 2*@ r *s + 2*@ r*@ s*r + 2*k*s 1563 t t t t t 1564ricci := -------------------------------------------- 1565 1 1 2 1566 r *s 1567 1568 2 1569 @ r*r*s + 2*@ r *s + 2*@ r*@ s*r + 2*k*s 1570 t t t t t 1571ricci := -------------------------------------------- 1572 2 2 2 1573 r *s 1574 1575 2 1576 @ r*r*s + 2*@ r *s + 2*@ r*@ s*r + 2*k*s 1577 t t t t t 1578ricci := -------------------------------------------- 1579 3 3 2 1580 r *s 1581 1582 2 1583 3*@ r*@ s*s + @ s*r*s + @ s *r + r 1584 t t t t t 1585ricci := -------------------------------------- 1586 4 4 2 1587 r*s 1588 1589 2 1590 3*@ r*@ s*s + @ s*r*s + @ s *r + r 1591 t t t t t 1592ricci := -------------------------------------- 1593 5 5 2 1594 r*s 1595 1596 1597 1598riccisc := ricci(-a,-b)*g(a,b); 1599 1600 1601 2 2 2 2 2 2 1602riccisc := (2*(3*@ r*r*s + 3*@ r *s + 6*@ r*@ s*r*s + 2*@ s*r *s + @ s *r 1603 t t t t t t t t 1604 1605 2 2 2 2 1606 + 3*k*s + r ))/(r *s ) 1607 1608 1609pform {laglanc,inv1,inv2} = 0; 1610 1611 1612 1613index_symmetries riemc3(k,l),riemri(k,l), 1614 hlang(k,l),einst(k,l): symmetric; 1615 1616 1617 1618pform {riemc3(i,j),riemri(i,j)}=0; 1619 1620 1621 1622riemc3(-i,-j) := riemann(-i,-k,-l,-m)*riemann(-j,k,l,m)$ 1623 1624 1625inv1 := riemc3(-i,-j)*g(i,j); 1626 1627 1628 2 2 4 4 4 2 2 2 2 2 4 1629inv1 := (4*(3*@ r *r *s + 3*@ r *s + 6*@ r *@ s *r *s + 6*@ r *k*s 1630 t t t t t t 1631 1632 2 4 2 4 4 2 4 2 4 4 4 4 1633 + 2*@ s *r *s + @ s *r + 2*@ s *r + 3*k *s + r ))/(r *s ) 1634 t t t t 1635 1636riemri(-i,-j) := 2*riemann(-i,-k,-j,-l)*ricci(k,l)$ 1637 1638 1639inv2 := ricci(-a,-b)*ricci(a,b); 1640 1641 1642 2 2 4 2 4 2 3 1643inv2 := (2*(6*@ r *r *s + 6*@ r*@ r *r*s + 6*@ r*@ r*@ s*r *s 1644 t t t t t t t t t 1645 1646 3 3 4 4 4 1647 + 6*@ r*@ s*r *s + 6*@ r*k*r*s + 6*@ r *s 1648 t t t t t t t 1649 1650 3 3 2 2 2 2 2 4 1651 + 12*@ r *@ s*r*s + 15*@ r *@ s *r *s + 12*@ r *k*s 1652 t t t t t 1653 1654 3 2 3 3 3 1655 + 6*@ r*@ s*@ s*r *s + 6*@ r*@ s *r *s + 12*@ r*@ s*k*r*s 1656 t t t t t t t t 1657 1658 3 2 4 2 2 4 1659 + 6*@ r*@ s*r *s + 3*@ s *r *s + 2*@ s*@ s *r *s 1660 t t t t t t t 1661 1662 4 4 4 2 4 2 4 4 4 4 1663 + 2*@ s*r *s + @ s *r + 2*@ s *r + 6*k *s + r ))/(r *s ) 1664 t t t t 1665 1666laglanc := (1/2)*(inv1 - 4*inv2 + riccisc**2); 1667 1668 1669 2 2 2 2 2 1670laglanc := (12*(@ r*@ r *s + 4*@ r*@ r*@ s*r*s + @ r*@ s *r + @ r*k*s 1671 t t t t t t t t t t t t 1672 1673 2 3 2 2 2 1674 + @ r*r + 2*@ r *@ s*s + 2*@ r *@ s*r*s + 3*@ r *@ s *r 1675 t t t t t t t t t 1676 1677 2 2 1678 + @ r *r + 2*@ r*@ s*@ s*r + 2*@ r*@ s*k*s + 2*@ s*k*r*s 1679 t t t t t t t t t 1680 1681 2 3 2 1682 + @ s *k*r + k*r))/(r *s ) 1683 t 1684 1685 1686 1687pform {einst(a,b),hlang(a,b)}=0; 1688 1689 1690 1691hlang(-i,-j) := 2*(riemc3(-i,-j) - riemri(-i,-j) - 1692 2*ricci(-i,-k)*ricci(-j,K) + 1693 riccisc*ricci(-i,-j) - (1/2)*laglanc*g(-i,-j)); 1694 1695hlang := 1696 t t 1697 1698 3 2 2 2 2 1699 12*(2*@ r *@ s*s + 3*@ r *@ s *r + @ r *r + 2*@ r*@ s*k*s + @ s *k*r + k*r) 1700 t t t t t t t t 1701----------------------------------------------------------------------------- 1702 3 2 1703 r *s 1704 1705 2 1706hlang := (4*( - 4*@ r*@ r*@ s*s - 2*@ r*@ s *r - 2*@ r*r 1707 1 1 t t t t t t t t t 1708 1709 2 2 2 2 1710 - 2*@ r *@ s*s - 3*@ r *@ s - @ r - 4*@ r*@ s*@ s*r 1711 t t t t t t t t t t 1712 1713 2 2 2 1714 - 2*@ s*k*s - @ s *k - k))/(r *s ) 1715 t t t 1716 1717 2 1718hlang := (4*( - 4*@ r*@ r*@ s*s - 2*@ r*@ s *r - 2*@ r*r 1719 2 2 t t t t t t t t t 1720 1721 2 2 2 2 1722 - 2*@ r *@ s*s - 3*@ r *@ s - @ r - 4*@ r*@ s*@ s*r 1723 t t t t t t t t t t 1724 1725 2 2 2 1726 - 2*@ s*k*s - @ s *k - k))/(r *s ) 1727 t t t 1728 1729 2 1730hlang := (4*( - 4*@ r*@ r*@ s*s - 2*@ r*@ s *r - 2*@ r*r 1731 3 3 t t t t t t t t t 1732 1733 2 2 2 2 1734 - 2*@ r *@ s*s - 3*@ r *@ s - @ r - 4*@ r*@ s*@ s*r 1735 t t t t t t t t t t 1736 1737 2 2 2 1738 - 2*@ s*k*s - @ s *k - k))/(r *s ) 1739 t t t 1740 1741 2 3 1742hlang := (12*( - @ r*@ r *s - 2*@ r*@ r*@ s*r - @ r*k*s - @ r *@ s 1743 4 4 t t t t t t t t t t t 1744 1745 2 3 1746 - @ r *@ s*r - @ r*@ s*k - @ s*k*r))/(r *s) 1747 t t t t t t t 1748 1749 2 3 1750hlang := (12*( - @ r*@ r *s - 2*@ r*@ r*@ s*r - @ r*k*s - @ r *@ s 1751 5 5 t t t t t t t t t t t 1752 1753 2 3 1754 - @ r *@ s*r - @ r*@ s*k - @ s*k*r))/(r *s) 1755 t t t t t t t 1756 1757 1758 1759% The complete Einstein tensor: 1760 1761einst(-i,-j) := (ricci(-i,-j) - (1/2)*riccisc*g(-i,-j))*alp1 + 1762 hlang(-i,-j)*alp2$ 1763 1764 1765 1766alp1 := 1$ 1767 1768 1769factor alp2; 1770 1771 1772 1773einst(-i,-j) := einst(-i,-j); 1774 1775 3 2 2 2 1776einst := (12*alp2*(2*@ r *@ s*s + 3*@ r *@ s *r + @ r *r + 2*@ r*@ s*k*s 1777 t t t t t t t t t 1778 1779 2 3 2 1780 + @ s *k*r + k*r))/(r *s ) 1781 t 1782 1783 2 2 2 2 2 2 1784 3*@ r *s + 6*@ r*@ s*r*s + @ s *r + 3*k*s + r 1785 t t t t 1786 + --------------------------------------------------- 1787 2 2 1788 r *s 1789 1790 2 1791einst := (4*alp2*( - 4*@ r*@ r*@ s*s - 2*@ r*@ s *r - 2*@ r*r 1792 1 1 t t t t t t t t t 1793 1794 2 2 2 2 1795 - 2*@ r *@ s*s - 3*@ r *@ s - @ r - 4*@ r*@ s*@ s*r 1796 t t t t t t t t t t 1797 1798 2 2 2 2 1799 - 2*@ s*k*s - @ s *k - k))/(r *s ) + ( - 2*@ r*r*s 1800 t t t t t 1801 1802 2 2 2 2 2 2 2 1803 - @ r *s - 4*@ r*@ s*r*s - 2*@ s*r *s - @ s *r - k*s - r )/ 1804 t t t t t t 1805 1806 2 2 1807 (r *s ) 1808 1809 2 1810einst := (4*alp2*( - 4*@ r*@ r*@ s*s - 2*@ r*@ s *r - 2*@ r*r 1811 2 2 t t t t t t t t t 1812 1813 2 2 2 2 1814 - 2*@ r *@ s*s - 3*@ r *@ s - @ r - 4*@ r*@ s*@ s*r 1815 t t t t t t t t t t 1816 1817 2 2 2 2 1818 - 2*@ s*k*s - @ s *k - k))/(r *s ) + ( - 2*@ r*r*s 1819 t t t t t 1820 1821 2 2 2 2 2 2 2 1822 - @ r *s - 4*@ r*@ s*r*s - 2*@ s*r *s - @ s *r - k*s - r )/ 1823 t t t t t t 1824 1825 2 2 1826 (r *s ) 1827 1828 2 1829einst := (4*alp2*( - 4*@ r*@ r*@ s*s - 2*@ r*@ s *r - 2*@ r*r 1830 3 3 t t t t t t t t t 1831 1832 2 2 2 2 1833 - 2*@ r *@ s*s - 3*@ r *@ s - @ r - 4*@ r*@ s*@ s*r 1834 t t t t t t t t t t 1835 1836 2 2 2 2 1837 - 2*@ s*k*s - @ s *k - k))/(r *s ) + ( - 2*@ r*r*s 1838 t t t t t 1839 1840 2 2 2 2 2 2 2 1841 - @ r *s - 4*@ r*@ s*r*s - 2*@ s*r *s - @ s *r - k*s - r )/ 1842 t t t t t t 1843 1844 2 2 1845 (r *s ) 1846 1847 2 3 1848einst := (12*alp2*( - @ r*@ r *s - 2*@ r*@ r*@ s*r - @ r*k*s - @ r *@ s 1849 4 4 t t t t t t t t t t t 1850 1851 2 3 1852 - @ r *@ s*r - @ r*@ s*k - @ s*k*r))/(r *s) 1853 t t t t t t t 1854 1855 2 2 1856 - 3*@ r*r*s - 3*@ r *s - 3*@ r*@ s*r - @ s*r - 3*k*s 1857 t t t t t t t 1858 + ------------------------------------------------------------ 1859 2 1860 r *s 1861 1862 2 3 1863einst := (12*alp2*( - @ r*@ r *s - 2*@ r*@ r*@ s*r - @ r*k*s - @ r *@ s 1864 5 5 t t t t t t t t t t t 1865 1866 2 3 1867 - @ r *@ s*r - @ r*@ s*k - @ s*k*r))/(r *s) 1868 t t t t t t t 1869 1870 2 2 1871 - 3*@ r*r*s - 3*@ r *s - 3*@ r*@ s*r - @ s*r - 3*k*s 1872 t t t t t t t 1873 + ------------------------------------------------------------ 1874 2 1875 r *s 1876 1877 1878 1879clear o(k),e(k),riemc3(i,j),riemri(i,j),curv(k,l),riemann(a,b,c,d), 1880 ricci(a,b),riccisc,t,u,v1,v2,theta,phi,r,om(k,l),einst(a,b), 1881 hlang(a,b); 1882 1883 1884 1885remfdomain r,s; 1886 1887 1888 1889% Problem: 1890% -------- 1891% Calculate for a given coframe and given torsion the Riemannian part and 1892% the torsion induced part of the connection. Calculate the curvature. 1893 1894% For a more elaborate example see E.Schruefer, F.W. Hehl, J.D. McCrea, 1895% "Application of the REDUCE package EXCALC to the Poincare gauge field 1896% theory of gravity", GRG Journal, vol. 19, (1988) 197--218 1897 1898pform {ff, gg}=0; 1899 1900 1901 1902fdomain ff=ff(r), gg=gg(r); 1903 1904 1905 1906coframe o(4) = d u + 2*b0*cos(theta)*d phi, 1907 o(1) = ff*(d u + 2*b0*cos(theta)*d phi) + d r, 1908 o(2) = gg*d theta, 1909 o(3) = gg*sin(theta)*d phi 1910 with metric g = -o(4)*o(1)-o(4)*o(1)+o(2)*o(2)+o(3)*o(3); 1911 1912 1913 1914frame e; 1915 1916 1917 1918pform {tor(a),gwt(a)}=2,gamma(a,b)=1, 1919 {u1,u3,u5}=0; 1920 1921 1922 1923index_symmetries gamma(a,b): antisymmetric; 1924 1925 1926 1927fdomain u1=u1(r),u3=u3(r),u5=u5(r); 1928 1929 1930 1931tor(4) := 0$ 1932 1933 1934 1935tor(1) := -u5*o(4)^o(1) - 2*u3*o(2)^o(3)$ 1936 1937 1938 1939tor(2) := u1*o(4)^o(2) + u3*o(4)^o(3)$ 1940 1941 1942 1943tor(3) := u1*o(4)^o(3) - u3*o(4)^o(2)$ 1944 1945 1946 1947gwt(-a) := d o(-a) - tor(-a)$ 1948 1949 1950 1951% The following is the combined connection. 1952% The Riemannian part could have equally well been calculated by the 1953% RIEMANNCONX statement. 1954 1955gamma(-a,-b) := (1/2)*( e(-b) _| (e(-c) _| gwt(-a)) 1956 +e(-c) _| (e(-a) _| gwt(-b)) 1957 -e(-a) _| (e(-b) _| gwt(-c)) )*o(c); 1958 1959 4 1960gamma := o *(@ ff - u5) 1961 4 1 r 1962 1963 2 1964 o *(@ gg*ff + gg*u1) 3 2 1965 r o *( - b0*ff + gg *u3) 1966gamma := - ---------------------- + ------------------------ 1967 4 2 gg 2 1968 gg 1969 1970 2 1971 o *@ gg 3 1972 r - o *b0 1973gamma := --------- + ---------- 1974 1 2 gg 2 1975 gg 1976 1977 3 1978 2 2 o *(@ gg*ff + gg*u1) 1979 o *(b0*ff - gg *u3) r 1980gamma := --------------------- - ---------------------- 1981 4 3 2 gg 1982 gg 1983 1984 3 1985 2 o *@ gg 1986 o *b0 r 1987gamma := ------- + --------- 1988 1 3 2 gg 1989 gg 1990 1991 1 3 4 2 1992 o *b0 o *cos(theta) o *(b0*ff - 2*gg *u3) 1993gamma := ------- + --------------- + ----------------------- 1994 2 3 2 sin(theta)*gg 2 1995 gg gg 1996 1997 1998 1999pform curv(a,b)=2; 2000 2001 2002index_symmetries curv(a,b): antisymmetric; 2003 2004 2005factor ^; 2006 2007 2008 2009curv(-a,b) := d gamma(-a,b) + gamma(-c,b)^gamma(-a,c); 2010 2011 4 2 3 2012curv := (2*o ^o 2013 4 2014 2015 2 2016 *(@ ff*b0*gg - 2*@ gg*b0*ff + @ gg*gg *u3 - b0*gg*u1 - b0*gg*u5))/ 2017 r r r 2018 2019 3 4 1 2020 gg + o ^o *(@ ff - @ u5) 2021 r r r 2022 2023 1 2 3 2 2024 o ^o *(@ gg*gg - b0 ) 2025 4 r r 4 2 3 3 2026curv := -------------------------- + (o ^o *( - @ ff*@ gg*gg - @ gg*ff*gg 2027 2 4 r r r r 2028 gg 2029 2030 3 2 2 4 2031 + @ gg*gg *u5 - b0 *ff + 2*b0*gg *u3))/gg 2032 r 2033 2034 4 3 2 2035 o ^o *(@ ff*b0*gg - 2*@ gg*b0*ff + 2*@ gg*gg *u3 - b0*gg*u5) 2036 r r r 2037 + -------------------------------------------------------------- 2038 3 2039 gg 2040 2041 1 3 3 2 2042 o ^o *(@ gg*gg - b0 ) 2043 4 r r 2044curv := -------------------------- 2045 3 4 2046 gg 2047 2048 4 2 2 2049 o ^o *( - @ ff*b0*gg + 2*@ gg*b0*ff - 2*@ gg*gg *u3 + b0*gg*u5) 2050 r r r 2051 + ----------------------------------------------------------------- 2052 3 2053 gg 2054 2055 4 3 3 3 3 2 2056 + (o ^o *( - @ ff*@ gg*gg - @ gg*ff*gg + @ gg*gg *u5 - b0 *ff 2057 r r r r r 2058 2059 2 4 2060 + 2*b0*gg *u3))/gg 2061 2062 1 2 3 2063curv := (2*o ^o 2064 1 2065 2066 2 2067 *( - @ ff*b0*gg + 2*@ gg*b0*ff - @ gg*gg *u3 + b0*gg*u1 + b0*gg*u5)) 2068 r r r 2069 2070 3 4 1 2071 /gg + o ^o *( - @ ff + @ u5) 2072 r r r 2073 2074 1 1 2 3 3 3 4 2075curv := (o ^o *( - @ ff*@ gg*gg - @ gg*ff*gg - @ gg*gg *u1 - @ u1*gg 2076 2 r r r r r r 2077 2078 2 2 4 1 3 2079 - b0 *ff + b0*gg *u3))/gg + (o ^o *( - @ ff*b0*gg 2080 r 2081 2082 2 3 3 2083 + 2*@ gg*b0*ff + @ gg*gg *u3 + @ u3*gg + b0*gg*u1))/gg + ( 2084 r r r 2085 2086 4 2 4 2 3 3 2087 o ^o *( - @ ff*gg *u1 + @ gg*ff *gg + @ gg*ff*gg *u1 2088 r r r r 2089 2090 3 4 2 2 2 2091 + @ gg*ff*gg *u5 + @ u1*ff*gg - b0 *ff + 3*b0*ff*gg *u3 2092 r r 2093 2094 4 4 2 4 4 3 2 2095 + gg *u1*u5 - 2*gg *u3 ))/gg + (o ^o *(@ ff*gg *u3 2096 r 2097 2098 2 2099 - 3*@ gg*ff*gg*u3 - @ u3*ff*gg + b0*ff*u1 + b0*ff*u5 2100 r r 2101 2102 2 2 2 2103 - 2*gg *u1*u3 - gg *u3*u5))/gg 2104 2105 1 1 2 2106curv := (o ^o 2107 3 2108 2109 2 3 2110 *(@ ff*b0*gg - 2*@ gg*b0*ff - @ gg*gg *u3 - @ u3*gg - b0*gg*u1))/ 2111 r r r r 2112 2113 3 1 3 3 3 3 2114 gg + (o ^o *( - @ ff*@ gg*gg - @ gg*ff*gg - @ gg*gg *u1 2115 r r r r r 2116 2117 4 2 2 4 4 2 2 2118 - @ u1*gg - b0 *ff + b0*gg *u3))/gg + (o ^o *( - @ ff*gg *u3 2119 r r 2120 2121 2 2122 + 3*@ gg*ff*gg*u3 + @ u3*ff*gg - b0*ff*u1 - b0*ff*u5 2123 r r 2124 2125 2 2 2 4 3 4 2126 + 2*gg *u1*u3 + gg *u3*u5))/gg + (o ^o *( - @ ff*gg *u1 2127 r 2128 2129 2 3 3 3 2130 + @ gg*ff *gg + @ gg*ff*gg *u1 + @ gg*ff*gg *u5 2131 r r r r 2132 2133 4 2 2 2 4 2134 + @ u1*ff*gg - b0 *ff + 3*b0*ff*gg *u3 + gg *u1*u5 2135 r 2136 2137 4 2 4 2138 - 2*gg *u3 ))/gg 2139 2140 2 1 2 3 3 3 4 2141curv := (o ^o *( - @ ff*@ gg*gg - @ gg*ff*gg - @ gg*gg *u1 - @ u1*gg 2142 4 r r r r r r 2143 2144 2 2 4 1 3 2145 - b0 *ff + b0*gg *u3))/gg + (o ^o *( - @ ff*b0*gg 2146 r 2147 2148 2 3 3 2149 + 2*@ gg*b0*ff + @ gg*gg *u3 + @ u3*gg + b0*gg*u1))/gg + ( 2150 r r r 2151 2152 4 2 4 2 3 3 2153 o ^o *( - @ ff*gg *u1 + @ gg*ff *gg + @ gg*ff*gg *u1 2154 r r r r 2155 2156 3 4 2 2 2 2157 + @ gg*ff*gg *u5 + @ u1*ff*gg - b0 *ff + 3*b0*ff*gg *u3 2158 r r 2159 2160 4 4 2 4 4 3 2 2161 + gg *u1*u5 - 2*gg *u3 ))/gg + (o ^o *(@ ff*gg *u3 2162 r 2163 2164 2 2165 - 3*@ gg*ff*gg*u3 - @ u3*ff*gg + b0*ff*u1 + b0*ff*u5 2166 r r 2167 2168 2 2 2 2169 - 2*gg *u1*u3 - gg *u3*u5))/gg 2170 2171 1 2 3 2 2172 o ^o *(@ gg*gg - b0 ) 2173 2 r r 4 2 3 3 2174curv := -------------------------- + (o ^o *( - @ ff*@ gg*gg - @ gg*ff*gg 2175 1 4 r r r r 2176 gg 2177 2178 3 2 2 4 2179 + @ gg*gg *u5 - b0 *ff + 2*b0*gg *u3))/gg 2180 r 2181 2182 4 3 2 2183 o ^o *(@ ff*b0*gg - 2*@ gg*b0*ff + 2*@ gg*gg *u3 - b0*gg*u5) 2184 r r r 2185 + -------------------------------------------------------------- 2186 3 2187 gg 2188 2189 2 2 3 2190curv := (o ^o 2191 3 2192 2193 2 2 3 2 2 2 2194 *( - 2*@ gg *ff*gg - 2*@ gg*gg *u1 + 6*b0 *ff - 6*b0*gg *u3 + gg )) 2195 r r 2196 2197 4 1 3 2198 2*o ^o *(@ ff*b0*gg - 2*@ gg*b0*ff - @ u3*gg ) 2199 4 r r r 2200 /gg + ------------------------------------------------ 2201 3 2202 gg 2203 2204 3 1 2 2205curv := (o ^o 2206 4 2207 2208 2 3 2209 *(@ ff*b0*gg - 2*@ gg*b0*ff - @ gg*gg *u3 - @ u3*gg - b0*gg*u1))/ 2210 r r r r 2211 2212 3 1 3 3 3 3 2213 gg + (o ^o *( - @ ff*@ gg*gg - @ gg*ff*gg - @ gg*gg *u1 2214 r r r r r 2215 2216 4 2 2 4 4 2 2 2217 - @ u1*gg - b0 *ff + b0*gg *u3))/gg + (o ^o *( - @ ff*gg *u3 2218 r r 2219 2220 2 2221 + 3*@ gg*ff*gg*u3 + @ u3*ff*gg - b0*ff*u1 - b0*ff*u5 2222 r r 2223 2224 2 2 2 4 3 4 2225 + 2*gg *u1*u3 + gg *u3*u5))/gg + (o ^o *( - @ ff*gg *u1 2226 r 2227 2228 2 3 3 3 2229 + @ gg*ff *gg + @ gg*ff*gg *u1 + @ gg*ff*gg *u5 2230 r r r r 2231 2232 4 2 2 2 4 2233 + @ u1*ff*gg - b0 *ff + 3*b0*ff*gg *u3 + gg *u1*u5 2234 r 2235 2236 4 2 4 2237 - 2*gg *u3 ))/gg 2238 2239 1 3 3 2 2240 o ^o *(@ gg*gg - b0 ) 2241 3 r r 2242curv := -------------------------- 2243 1 4 2244 gg 2245 2246 4 2 2 2247 o ^o *( - @ ff*b0*gg + 2*@ gg*b0*ff - 2*@ gg*gg *u3 + b0*gg*u5) 2248 r r r 2249 + ----------------------------------------------------------------- 2250 3 2251 gg 2252 2253 4 3 3 3 3 2 2254 + (o ^o *( - @ ff*@ gg*gg - @ gg*ff*gg + @ gg*gg *u5 - b0 *ff 2255 r r r r r 2256 2257 2 4 2258 + 2*b0*gg *u3))/gg 2259 2260 3 2 3 2261curv := (o ^o 2262 2 2263 2264 2 2 3 2 2 2 2265 *(2*@ gg *ff*gg + 2*@ gg*gg *u1 - 6*b0 *ff + 6*b0*gg *u3 - gg ))/ 2266 r r 2267 2268 4 1 3 2269 2*o ^o *( - @ ff*b0*gg + 2*@ gg*b0*ff + @ u3*gg ) 2270 4 r r r 2271 gg + --------------------------------------------------- 2272 3 2273 gg 2274 2275 2276 2277clear o(k),e(k),curv(a,b),gamma(a,b),theta,phi,x,y,z,r,s,t,u,v,p,q,c,cs; 2278 2279 2280remfdomain u1,u3,u5,ff,gg; 2281 2282 2283 2284showtime; 2285 2286 2287 2288end; 2289 2290Tested on x86_64-pc-windows CSL 2291Time (counter 1): 203 ms plus GC time: 16 ms 2292 2293End of Lisp run after 0.20+0.07 seconds 2294real 0.42 2295user 0.01 2296sys 0.04 2297