Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with
or
.
Download ZIP
Newer
Older
100644 204 lines (172 sloc) 8.338 kb
3990de3 @ept ML program to calculate derivatives and to do symbolic algebra.
authored
1 infix 6 power;
2 infix 5 times;
3 infix 4 cross;
4 infix 3 plus;
5
6 datatype expr =
7 V of string |
8 N of real |
9 times of expr*expr |
10 cross of expr*expr |
11 plus of expr*expr |
12 power of expr*real |
13 ddt of expr |
14 dual of expr |
15 trans of expr;
16
17
18 fun equals (V x) (V y) = (x = y)
19 | equals (N x) (N y) = (x = y)
4cb7127 @ept Lots of changes and fixes to algebra program, but in retrospect I'm n…
authored
20 | equals (a times b) (x times y) = (equals a x) andalso (equals b y)
3990de3 @ept ML program to calculate derivatives and to do symbolic algebra.
authored
21 | equals (a cross b) (x cross y) = (equals a x) andalso (equals b y)
22 | equals (a plus b) (x plus y) = ((equals a x) andalso (equals b y)) orelse
23 ((equals a y) andalso (equals a y))
24 | equals (a power b) (x power y) = (equals a x) andalso b=y
25 | equals (ddt x ) (ddt y ) = (equals x y)
26 | equals (dual x ) (dual y ) = (equals x y)
27 | equals (trans x ) (trans y ) = (equals x y)
28 | equals _ _ = false;
29
30
31 fun summands (a plus b) = (summands a) @ (summands b)
32 | summands a = [a];
33
34 fun countInSum(x, y::ys) = if equals x y then 1+(countInSum(x,ys)) else countInSum(x,ys)
35 | countInSum _ = 0;
36
37 fun combineSummands(prev, x::xs) =
38 if countInSum(x,prev) > 0 then combineSummands(prev, xs) else
39 let val cnt = countInSum(x,xs)
40 in let val rest = combineSummands(x::prev, xs) and
41 this = (if cnt = 0 then x else ((N (real(1+cnt))) times x))
42 in case rest of (N 0.0) => this | _ => (this plus rest) end end
43 | combineSummands(_, []) = N 0.0;
44
45
4cb7127 @ept Lots of changes and fixes to algebra program, but in retrospect I'm n…
authored
46 fun multExpr (N x, N y) = N(x*y)
47 | multExpr (N x, e ) = if x = 0.0 then N 0.0 else
48 if x = 1.0 then e else (N x) times e
49 | multExpr (e, N x) = if x = 0.0 then N 0.0 else
50 if x = 1.0 then e else (N x) times e
51 | multExpr (e, f ) = e times f;
52
53
54 fun simplifyFactors ((N x) times e) =
55 let val (ep, c) = simplifyFactors e in (ep, multExpr(N x, c)) end
56 | simplifyFactors (e times (N x)) =
57 let val (ep, c) = simplifyFactors e in (ep, multExpr(N x, c)) end
58 | simplifyFactors ((a power x) times e) =
59 let val (ap, b) = simplifyFactors a and
60 (ep, c) = simplifyFactors e in (ep, multExpr(ap power x, c)) end
61 | simplifyFactors (e times (a power x)) =
62 let val (ap, b) = simplifyFactors a and
63 (ep, c) = simplifyFactors e in (ep, multExpr(ap power x, c)) end
64 | simplifyFactors (e times f) =
65 let val (ep, c) = simplifyFactors e and
66 (fp, d) = simplifyFactors f in (ep times fp, multExpr(c, d)) end
67 | simplifyFactors (e cross f) =
68 let val (ep, c) = simplifyFactors e and
69 (fp, d) = simplifyFactors f in (ep cross fp, multExpr(c, d)) end
70 | simplifyFactors (e plus f) =
71 let val (ep, c) = simplifyFactors e and
72 (fp, d) = simplifyFactors f in
73 ((multExpr(ep, c)) plus (multExpr(fp, d)), N 1.0) end
74 | simplifyFactors (a power b) =
75 let val (ap, c) = simplifyFactors a in
76 ((multExpr(ap, c)) power b, N 1.0) end
77 | simplifyFactors (ddt a) =
78 let val (ap, c) = simplifyFactors a in (ddt ap, c) end
79 | simplifyFactors (dual a) =
80 let val (ap, c) = simplifyFactors a in (dual ap, c) end
81 | simplifyFactors (trans a) =
82 let val (ap, c) = simplifyFactors a in (trans a, c) end
83 | simplifyFactors e = (e, N 1.0);
84
85
86 fun simplify expr =
87 let fun simp (a times b) = (simp a) times (simp b)
88 | simp (a cross b) = (simp a) cross (simp b)
89 | simp (a plus b) =
90 let val sa = simp a and sb = simp b
91 in case (sa,sb) of (N x, N y) => N (x+y)
92 | (N x, e ) => if x = 0.0 then e else (N x) plus e
93 | (e, N y) => if y = 0.0 then e else (N y) plus e
94 | _ => combineSummands([], summands(sa plus sb)) end
95 | simp (a power b) = (simp a) power b
96 | simp (ddt e) = ddt (simp e)
97 | simp (dual e) = dual (simp e)
98 | simp (trans(dual e)) = (N ~1.0) times (dual(simp e))
99 | simp (trans(trans e)) = simp e
100 | simp (trans(e times f)) = (simp(trans f)) times (simp(trans e))
101 | simp (trans(N x)) = N x
102 | simp (trans(a power b)) = (simp a) power b
103 | simp (trans(e)) = trans(simp e)
104 | simp e = e in
105 let val (ep, c) = simplifyFactors(simp expr)
106 in simp(multExpr(ep, c)) end end;
3990de3 @ept ML program to calculate derivatives and to do symbolic algebra.
authored
107
108
109 fun diff ((dn,de)::ds) (V n) = if dn = n then de else (diff ds (V n))
110 | diff _ (N n) = N 0.0
111 | diff ds (a times b) = (diff ds a) times b plus a times (diff ds b)
112 | diff ds (a cross b) = (diff ds a) cross b plus a cross (diff ds b)
113 | diff ds (a plus b) = (diff ds a) plus (diff ds b)
114 | diff ds (a power b) = (N b) times (a power (b-1.0)) times (diff ds a)
4cb7127 @ept Lots of changes and fixes to algebra program, but in retrospect I'm n…
authored
115 | diff ds (dual a) = dual (diff ds a)
116 | diff ds (trans a) = trans (diff ds a)
3990de3 @ept ML program to calculate derivatives and to do symbolic algebra.
authored
117 | diff _ e = ddt(e);
118
119
120 fun combine f _ [] = N 0.0
121 | combine f [] _ = N 0.0
122 | combine f (x::[]) (y::[]) = f(x,y)
123 | combine f (x::xs) (y::[]) = (f(x,y)) plus (combine f xs [y])
124 | combine f xs (y::ys) = (combine f xs [y]) plus (combine f xs ys);
125
126 fun sumOfProducts (a times b) = combine (fn (a,b) => a times b)
127 (summands (sumOfProducts a)) (summands (sumOfProducts b))
128 | sumOfProducts (a cross b) = combine (fn (a,b) => a cross b)
129 (summands (sumOfProducts a)) (summands (sumOfProducts b))
130 | sumOfProducts (a plus b) = (sumOfProducts a) plus (sumOfProducts b)
4cb7127 @ept Lots of changes and fixes to algebra program, but in retrospect I'm n…
authored
131 | sumOfProducts (trans a) = foldr (fn (x,y) => x plus y) (N 0.0)
132 (map (fn x => trans x) (summands(sumOfProducts a)))
3990de3 @ept ML program to calculate derivatives and to do symbolic algebra.
authored
133 | sumOfProducts e = e;
134
135
136 fun exists f [] = false
137 | exists f (x::xn) = (f x) orelse (exists f xn);
138
139 fun exprInProduct vars e = (exists (equals e) vars) orelse (
140 case e of
141 (a times b) => (exprInProduct vars a) orelse (exprInProduct vars b)
142 | (a cross b) => (exprInProduct vars a) orelse (exprInProduct vars b)
143 | (a plus b) => (exprInProduct vars a) orelse (exprInProduct vars b)
144 | (ddt a) => exprInProduct vars a
4cb7127 @ept Lots of changes and fixes to algebra program, but in retrospect I'm n…
authored
145 | (trans a) => exprInProduct vars a
3990de3 @ept ML program to calculate derivatives and to do symbolic algebra.
authored
146 | _ => false);
147
148 fun varToRight vars (a times b) =
149 if (exprInProduct vars a) then
4cb7127 @ept Lots of changes and fixes to algebra program, but in retrospect I'm n…
authored
150 trans((trans (varToRight vars b)) times (trans(varToRight vars a)))
151 else (varToRight vars a) times (varToRight vars b)
3990de3 @ept ML program to calculate derivatives and to do symbolic algebra.
authored
152 | varToRight vars (a cross b) =
153 if (exprInProduct vars a) then
154 (N ~1.0) times ((varToRight vars b) cross (varToRight vars a))
155 else ((varToRight vars a) cross (varToRight vars b))
156 | varToRight _ e = e;
157
158
159 fun toMatrices (a times b) = (toMatrices a) times (toMatrices b)
160 | toMatrices (a cross b) = (dual(toMatrices a)) times (toMatrices b)
161 | toMatrices (a plus b) = (toMatrices a) plus (toMatrices b)
162 | toMatrices (a power b) = (toMatrices a) power b
163 | toMatrices (ddt e) = ddt (toMatrices e)
164 | toMatrices (dual e) = dual (toMatrices e)
165 | toMatrices (trans e) = trans(toMatrices e)
166 | toMatrices e = e;
167
168
4cb7127 @ept Lots of changes and fixes to algebra program, but in retrospect I'm n…
authored
169 fun massageVars vars e =
3990de3 @ept ML program to calculate derivatives and to do symbolic algebra.
authored
170 foldr (fn (x,y) => x plus y) (N 0.0)
4cb7127 @ept Lots of changes and fixes to algebra program, but in retrospect I'm n…
authored
171 (map (varToRight vars) (summands(sumOfProducts(simplify e))));
172
173 fun massage e = simplify(
174 massageVars [V "w", ddt(V "w"), V "p", ddt(V "p")] e);
3990de3 @ept ML program to calculate derivatives and to do symbolic algebra.
authored
175
176 fun print (V n, _) = n
177 | print (N n, _) = makestring n
178 | print (a times b, level) =
179 (if (level > 3) andalso (level <> 5) then "(" else "") ^
180 (print(a, 5)) ^ " " ^ (print(b, 5)) ^
181 (if (level > 3) andalso (level <> 5) then ")" else "")
182 | print (a cross b, level) =
183 (if level > 3 then "(" else "") ^
184 (print(a, 4)) ^ " x " ^ (print(b,4)) ^
185 (if level > 3 then ")" else "")
4cb7127 @ept Lots of changes and fixes to algebra program, but in retrospect I'm n…
authored
186 | print (a plus b, level) =
187 (if level > 3 then "(" else "") ^
188 (print(a,3)) ^ " + " ^ (print(b,3)) ^
189 (if level > 3 then ")" else "")
3990de3 @ept ML program to calculate derivatives and to do symbolic algebra.
authored
190 | print (a power b, _) = (print(a,6)) ^ "^" ^ (makestring b)
191 | print (ddt a, _) = (print(a,6)) ^ "'"
192 | print (dual a, _) = (print(a,6)) ^ "*"
193 | print (trans a, _) = (print(a,6)) ^ "T";
194
195 fun prettyprint e = print(e,3);
196
197
198 val derivs = [("e", (V "w") cross (V "e")), ("f", (V "p") cross (V "f"))];
199 val n = (V "e") cross (V "f");
4cb7127 @ept Lots of changes and fixes to algebra program, but in retrospect I'm n…
authored
200 val nm = (trans n) times n;
201 val nmd = diff derivs nm;
3990de3 @ept ML program to calculate derivatives and to do symbolic algebra.
authored
202 val nhat = n times (((trans n) times n) power ~0.5);
4cb7127 @ept Lots of changes and fixes to algebra program, but in retrospect I'm n…
authored
203 prettyprint(massage nmd);
Something went wrong with that request. Please try again.