-
Notifications
You must be signed in to change notification settings - Fork 3
/
riemannian.rkt
166 lines (148 loc) · 8.27 KB
/
riemannian.rkt
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
#lang racket
(require (only-in "fundamental.rkt" list-take list-reverse map-n))
(require "simplify.rkt"
"tensor.rkt"
"linear-algebra.rkt")
(provide (all-defined-out))
;;;
(define (identity-mat-as-tensor index-lst dim)
(cond ((not (= 2 (length index-lst))) (error "Length of index-lst doesn't match" index-lst))
((eq? (car (car index-lst)) (car (cadr index-lst))) (error "Indices should be one sub one super" index-lst))
(else (make-tensor index-lst (identity-mat dim)))))
;(identity-mat-as-tensor '((_ a) (^ b)) 4) ;works
(define (tensor-order tnsr) (length (get-index tnsr)))
(define (tensor-dimension tnsr) (length (get-matrix tnsr)))
;For tensors in tensor.rkt, the dimensions belong to different indices need not be the same.
;In Riemannian geometry, the dimension (of space) should be the same for all orders.
;In here, we don't check whether it is true or not.
(define (einstein-summation tnsr)
(define (index-same-back ele lst)
(define (list-flip new-lst old-lst)
(cond ((null? old-lst) false)
((eq? (cadr ele) (cadr (car old-lst))) (list-reverse new-lst (append (cdr old-lst) (list (car old-lst)))))
(else (list-flip (cons (car old-lst) new-lst) (cdr old-lst)))))
(list-flip '() lst))
;(index-same-back '(^ a) '((_ b) (_ a) (^ c))) ;'((_ b) (^ c) (_ a))
(define (list-same-right lst)
(define (list-same-right-recur lst)
(if (null? lst)
false
(let ([remain (index-same-back (car lst) (cdr lst))])
(if remain
(append remain (list (car lst)))
(let ([recur (list-same-right-recur (cdr lst))])
(if recur
(cons (car lst) recur)
false))))))
(list-same-right-recur lst))
(list-same-right '((_ a) (_ b) (^ a) (_ d) (^ b))) ;'((_ b) (_ d) (^ b) (^ a) (_ a))
(list-same-right '((_ a) (_ b) (^ c) (_ d))) ;#f
(define (tensor-trace tnsr)
(let ([len (- (length (get-index tnsr)) 2)])
(make-tensor
(list-take (get-index tnsr) len)
(map-n len mat-trace (get-matrix-without-tag tnsr)))))
(if (= (length (get-index tnsr)) 2)
(if (eq? (cadr (car (get-index tnsr))) (cadr (cadr (get-index tnsr))))
(make-scalar (simplify (mat-trace (get-matrix-without-tag tnsr))))
(simplify-generic tnsr))
(let ([new-index (list-same-right (get-index tnsr))])
(if new-index
(simplify-generic (einstein-summation (tensor-trace (switch-index new-index tnsr))))
(simplify-generic tnsr)))))
;(define t (make-tensor '((^ a) (_ a)) '((a b) (c d))))
;(einstein-summation t) ;'(scalar + a d)
;(define s (make-tensor '((^ a) (_ b)) '((a b) (c d))))
;(einstein-summation s)
;(define ts (make-tensor '((_ b) (^ a) (^ b)) (list (list (list 1 2) (list 3 4)) (list (list 5 6) (list 7 8)))))
;(einstein-summation ts) ;'(tensor (a) (scalar . 7) (scalar . 11))
;(define tss (make-tensor '((_ a) (^ b) (^ c)) (list (list (list 1 2) (list 3 4)) (list (list 5 6) (list 7 8)))))
;(einstein-summation tss) ;same as the original
(define (metric upper-lower-lst tnsr)
(cond ((not (= 2 (tensor-order tnsr))) (error "Not a metric" tnsr))
((not (= 2 (length upper-lower-lst))) (error "Upper-lower-lst doesn't match" upper-lower-lst))
((or (equal? upper-lower-lst '(_ ^)) (equal? upper-lower-lst '(^ _)))
(identity-mat-as-tensor (list (list (car upper-lower-lst) (cadr (car (get-index tnsr))))
(list (cadr upper-lower-lst) (cadr (cadr (get-index tnsr)))))
(tensor-dimension tnsr)))
((and (eq? (car upper-lower-lst) (car (car (get-index tnsr))))
(eq? (cadr upper-lower-lst) (car (cadr (get-index tnsr))))) tnsr)
(else
(make-tensor (list (list (car upper-lower-lst) (cadr (car (get-index tnsr))))
(list (cadr upper-lower-lst) (cadr (cadr (get-index tnsr)))))
(mat-inverse (get-matrix-without-tag tnsr))))))
;(define gg (make-tensor '((_ a) (_ b)) '((1 -3 1) (-3 -1 0) (1 0 2))))
;(define gg (make-tensor '((_ a) (_ b)) '((1 -3 a) (-3 -1 0) (a 0 2))))
;The second one also works, but it doesn't know how to do simplification right now.
;(metric '(_ ^) gg) ;works
;(metric '(_ _) gg) ;works
;(metric '(^ ^) gg) ;works
;(einstein-summation (mul (change-index '((^ a) (^ b)) (metric '(^ ^) gg))
; (change-index '((_ b) (_ c)) (metric '(_ _) gg)))) ;'(tensor ((^ a) (_ c)) identity)
;successfully got the simple diagonal matrix :-)
;'(tensor ((^ a) (_ c)) ((scalar . 1) (scalar . 0) (scalar . 0)) ((scalar . 0) (scalar . 1) (scalar . 0)) ((scalar . 0) (scalar . 0) (scalar . 1)))
(define (christoffel index-lst g-tensor coordinate-lst)
(if (nand (= 3 (length index-lst))
(eq? (car (car index-lst)) '^)
(eq? (car (cadr index-lst)) '_)
(eq? (car (caddr index-lst)) '_))
(error "Christoffel symbol needs 3 indices, has super-sub-sub in order. You give" index-lst)
(let ([first-term (partial-deriv (change-index '((_ dummy) (_ j)) g-tensor)
(make-tensor '((_ k)) coordinate-lst))])
(change-index
index-lst
(einstein-summation
(mul
(change-index '((^ i) (^ dummy)) (metric '(^ ^) g-tensor))
(scalar-mul
(/ 1 2)
(add (add first-term
(switch-index '((_ dummy) (_ j) (_ k)) (change-index '((_ dummy) (_ k) (_ j)) first-term)))
(scalar-mul -1 (switch-index '((_ dummy) (_ j) (_ k)) (change-index '((_ k) (_ j) (_ dummy)) first-term)))))))))))
(define (riemann-tensor index-lst christoffel-gamma coordinate-lst)
(if (nand (= 4 (length index-lst))
(eq? (car (car index-lst)) '^)
(eq? (car (cadr index-lst)) '_)
(eq? (car (caddr index-lst)) '_)
(eq? (car (cadddr index-lst)) '_))
(error "Riemann curvature tensor needs 4 indices, has super-sub-sub-sub in order. You give" index-lst)
(let* ([gamma (change-index '((^ i) (_ j) (_ k)) christoffel-gamma)]
[partial-gamma (partial-deriv gamma (make-tensor '((_ l)) coordinate-lst))])
(change-index
index-lst
(add
(add (switch-index '((^ i) (_ k) (_ l) (_ j)) partial-gamma)
(scalar-mul -1 (switch-index '((^ i) (_ k) (_ l) (_ j))
(change-index '((^ i) (_ l) (_ k) (_ j)) partial-gamma))))
(add
(switch-index '((^ i) (_ k) (_ l) (_ j))
(einstein-summation
(mul (change-index '((^ i) (_ l) (_ dummy)) gamma)
(change-index '((^ dummy) (_ j) (_ k)) gamma))))
(switch-index '((^ i) (_ k) (_ l) (_ j))
(scalar-mul -1 (einstein-summation
(mul (change-index '((^ i) (_ j) (_ dummy)) gamma)
(change-index '((^ dummy) (_ l) (_ k)) gamma)))))))))))
(define (ricci-curvature-tensor index-lst riemann-tnsr)
(if (nand (= 2 (length index-lst))
(eq? (car (car index-lst)) '_)
(eq? (car (cadr index-lst)) '_))
(error "Ricci curvature tensor needs 2 indices, has sub-sub in order. You give" index-lst)
(change-index
index-lst
(einstein-summation
(change-index '((^ k) (_ i) (_ k) (_ j)) riemann-tnsr)))))
(define (ricci-scalar g-tnsr ricci-tnsr)
(einstein-summation
(mul (change-index '((^ i) (^ j)) (metric '(^ ^) g-tnsr))
(change-index '((_ j) (_ i)) ricci-tnsr))))
;(define g (make-tensor '((_ a) (_ b)) '((x1 0) (0 x2))))
;(define gamma (christoffel '((^ i) (_ j) (_ k)) g '(x1 x2))) ;seems work, no simplification right now.
;'(tensor
; ((^ i) (_ j) (_ k))
; (((scalar * (1/2)(* x2 (** (* x2 x1) -1))) (scalar . 0)) ((scalar . 0) (scalar . 0)))
; (((scalar . 0) (scalar . 0)) ((scalar . 0) (scalar * (1/2)(* x1 (** (* x2 x1) -1))))))
;;;(christoffel '((^ i) (_ j) (^ k)) g '(x1 x2)) ;index error
;(define r_abcd (riemann-tensor '((^ a) (_ b) (_ c) (_ d)) gamma '(x1 x2)))
;(define r_ab (ricci-curvature-tensor '((_ a) (_ b)) r_abcd)) ;It is symmetric right now.
;(ricci-scalar g r_ab) ;works. no simplification. so currently can't check whether right or now.