1 (* ========================================================================== *)
2 (* FLYSPECK - CODE INFORMAL *)
4 (* Chapter: Nonlinear *)
5 (* Author: Thomas C. Hales *)
7 (* ========================================================================== *)
11 Define the functions mdtau and mdtau2.
14 mdtau = Sqrt[delta] * D[taumar,y1]
15 However, we need to rework the formula in such a way that it extends to delta=0.
17 mdtau2 = D[taumar,{y1,2}].
19 mdtau2uf = mdtau2 * uf, where uf is the positive factor defined below.
21 In /interval_code, taylorInterval.cc, wide.cc, there is a version mdtau_y_LC, mdtau2_y_LC
22 The "LC" signifies "locally constant", meaning no derivatives are to be used in interval computations.
24 ("mdtau_y_LC","Lib::mdtau_y_LC");
25 ("mdtau2uf_y_LC","Lib::mdtau2uf_y_LC");
26 ("mdtau_y","Lib::mdtau_y_LC");
27 ("mdtau2uf_y","Lib::mdtau2uf_y_LC");
33 let dua = new_definition `dua a b c = &2 * (b + c - a)`;;
35 (* corrected June 25, 2012, I had left out the sqrt! *)
37 let safesqrt = new_definition `safesqrt x = if (x >= &0) then sqrt x else (&0)`;;
39 let mdtau_y = new_definition `mdtau_y y1 y2 y3 y4 y5 y6 =
46 let chain0 = &2 * y1 in
48 let chain2 = &4 * x1 in
49 let u135 = ups_x x1 x3 x5 in
50 let u126 = ups_x x1 x2 x6 in
51 let u234 = ups_x x2 x3 x4 in
52 let uf = &4 * u135 * u126 * u234 * y1 * y2 * y3 in
53 let du135 = dua x1 x3 x5 in
54 let du126 = dua x1 x2 x6 in
55 let Luf = (du135 / u135 + du126/u126 ) * chain0 + &1 / y1 in
56 let n4 = x2*x3 + x1*x4 - x2*x5 - x3*x6 + x5*x6 -
57 x1*(-- x1 + x2 + x3 - x4 + x5 + x6) in
59 let n5 = x1*x3 - x1*x4 + x2*x5 - x3*x6 + x4*x6 -
60 x2*(x1 - x2 + x3 + x4 - x5 + x6) in // - del5
61 let n6 = x1*x2 - x1*x4 - x2*x5 + x4*x5 -
62 x3*(x1 + x2 - x3 + x4 + x5 - x6) + x3*x6 in // - del6
63 let Dn4 = &2* x1 - x2 - x3 + &2 *x4 - x5 - x6 in
64 let del = delta_x x1 x2 x3 x4 x5 x6 in
65 let del1 = -- (x1*x4) + x2*x5 - x3*x5 - x2*x6 + x3*x6 +
66 x4*(-- x1 + x2 + x3 - x4 + x5 + x6) in
67 let del2 = x1*x4 - x3*x4 - x2*x5 - x1*x6 +
68 x3*x6 + x5*(x1 - x2 + x3 + x4 - x5 + x6) in
69 let del3 = x1*x4 - x2*x4 - x1*x5 + x2*x5 -
70 x3*x6 + (x1 + x2 - x3 + x4 + x5 - x6)*x6 in
71 let Pdel = del1 * chain0 in
72 let Ldel = Pdel/del in
73 let sd4 = &4 *x1*del in
74 let sd5 = &4 *x2*del in
75 let sd6 = &4 *x3*del in
76 let Dsd4 = &4 *del + &4 *x1*del1 in
77 let m4diff = &2 *Dn4*sd4 - n4* Dsd4 in
78 let m4 = m4diff*chain0*u234*y2*y3 in
79 let m5 = -- &4 *x2*u234*del3* &2 *x1*u135*y3 in
80 let m6 = -- &4 *x3*u234*del2* &2 *x1*u126*y2 in
81 let const1 = (sol_y (&2) (&2) (&2) (&2) (&2) (&2)) /pi in
82 let rhoy1 = rho(y1) in
83 let rhoy2 = rho(y2) in
84 let rhoy3 = rho(y3) in
85 let Prhoy1 = const1 /(#0.52) in
86 let rr = rhoy1 * m4 + rhoy2 * m5 + rhoy3 * m6 in
87 let term1 = Prhoy1 * pi * safesqrt(del) in
88 let t = sqrt(&4 * x1)/del4 in
90 let term2a = del * t * matan(t2 *del) in
91 let term2 = term2a * Prhoy1 in
96 let mdtau2_y = new_definition `mdtau2_y y1 y2 y3 y4 y5 y6 =
103 let chain0 = &2 * y1 in
105 let chain2 = &4 * x1 in
106 let u135 = ups_x x1 x3 x5 in
107 let u126 = ups_x x1 x2 x6 in
108 let u234 = ups_x x2 x3 x4 in
109 let uf = &4 * u135 * u126 * u234 * y1 * y2 * y3 in
110 let du135 = dua x1 x3 x5 in
111 let du126 = dua x1 x2 x6 in
112 let Luf = (du135 / u135 + du126/u126 ) * chain0 + &1 / y1 in
113 let n4 = x2*x3 + x1*x4 - x2*x5 - x3*x6 + x5*x6 -
114 x1*(-- x1 + x2 + x3 - x4 + x5 + x6) in
116 let n5 = x1*x3 - x1*x4 + x2*x5 - x3*x6 + x4*x6 -
117 x2*(x1 - x2 + x3 + x4 - x5 + x6) in // - del5
118 let n6 = x1*x2 - x1*x4 - x2*x5 + x4*x5 -
119 x3*(x1 + x2 - x3 + x4 + x5 - x6) + x3*x6 in // - del6
120 let Dn4 = &2* x1 - x2 - x3 + &2 *x4 - x5 - x6 in
121 let del = delta_x x1 x2 x3 x4 x5 x6 in
122 let del1 = -- (x1*x4) + x2*x5 - x3*x5 - x2*x6 + x3*x6 +
123 x4*(-- x1 + x2 + x3 - x4 + x5 + x6) in
124 let del2 = x1*x4 - x3*x4 - x2*x5 - x1*x6 +
125 x3*x6 + x5*(x1 - x2 + x3 + x4 - x5 + x6) in
126 let del3 = x1*x4 - x2*x4 - x1*x5 + x2*x5 -
127 x3*x6 + (x1 + x2 - x3 + x4 + x5 - x6)*x6 in
128 let Pdel = del1 * chain0 in
129 let Ldel = Pdel/del in
130 let sd4 = &4 *x1*del in
131 let sd5 = &4 *x2*del in
132 let sd6 = &4 *x3*del in
133 let Dsd4 = &4 *del + &4 *x1*del1 in
134 let m4diff = &2 *Dn4*sd4 - n4* Dsd4 in
135 let m4 = m4diff*chain0*u234*y2*y3 in
136 let m5 = -- &4 *x2*u234*del3* &2 *x1*u135*y3 in
137 let m6 = -- &4 *x3*u234*del2* &2 *x1*u126*y2 in
138 let const1 = (sol_y (&2) (&2) (&2) (&2) (&2) (&2)) /pi in
139 let rhoy1 = rho(y1) in
140 let rhoy2 = rho(y2) in
141 let rhoy3 = rho(y3) in
142 let Prhoy1 = const1 /(#0.52) in
143 let rr = rhoy1 * m4 + rhoy2 * m5 + rhoy3 * m6 in
145 let D2sd4 = -- &8*x1*x4 + &8*(--(x1*x4) + x2*x5 - x3*x5 - x2*x6 + x3*x6 +
146 x4*(-- x1 + x2 + x3 - x4 + x5 + x6)) in
147 let Dm4diff = &2 * D2n4 * sd4 + Dn4 * Dsd4 - n4 *D2sd4 in
148 let Pm4 = (Dm4diff * chain2 + m4diff * Pchain ) * u234 * y2 * y3 in
149 let Ddel3 = x4 - x5 + x6 in
150 let Ddel2 = x4 + x5 - x6 in
151 let Pm5 = (Ddel3 * x1 * u135 + del3 * &1 * u135 + del3 * x1 * du135) *
152 chain0 * (-- &4 * x2 * u234 * &2 * y3) in
153 let Pm6 = (Ddel2 * x1 * u126 + del2 * &1 * u126 + del2 * x1 * du126) *
154 chain0 * (-- &4 * x3 * u234 * &2 * y2) in
155 let PrrC = &2 * Prhoy1 * m4 + rhoy1 * Pm4 + rhoy2 * Pm5 + rhoy3 * Pm6 in
156 let P2tauNum = (PrrC) + (-- Luf - #0.5 * Ldel) * rr in
157 let P2tau = P2tauNum/ (uf * safesqrt(del)) in
163 mdtau2uf = mdtau2 * uf, where uf = &4 * u135 * u126 * u234 * y1 * y2 * y3
166 let mdtau2uf_y = new_definition `mdtau2uf_y y1 y2 y3 y4 y5 y6 =
173 let chain0 = &2 * y1 in
175 let chain2 = &4 * x1 in
176 let u135 = ups_x x1 x3 x5 in
177 let u126 = ups_x x1 x2 x6 in
178 let u234 = ups_x x2 x3 x4 in
179 let uf = &4 * u135 * u126 * u234 * y1 * y2 * y3 in
180 let du135 = dua x1 x3 x5 in
181 let du126 = dua x1 x2 x6 in
182 let Luf = (du135 / u135 + du126/u126 ) * chain0 + &1 / y1 in
183 let n4 = x2*x3 + x1*x4 - x2*x5 - x3*x6 + x5*x6 -
184 x1*(-- x1 + x2 + x3 - x4 + x5 + x6) in
186 let n5 = x1*x3 - x1*x4 + x2*x5 - x3*x6 + x4*x6 -
187 x2*(x1 - x2 + x3 + x4 - x5 + x6) in // - del5
188 let n6 = x1*x2 - x1*x4 - x2*x5 + x4*x5 -
189 x3*(x1 + x2 - x3 + x4 + x5 - x6) + x3*x6 in // - del6
190 let Dn4 = &2* x1 - x2 - x3 + &2 *x4 - x5 - x6 in
191 let del = delta_x x1 x2 x3 x4 x5 x6 in
192 let del1 = -- (x1*x4) + x2*x5 - x3*x5 - x2*x6 + x3*x6 +
193 x4*(-- x1 + x2 + x3 - x4 + x5 + x6) in
194 let del2 = x1*x4 - x3*x4 - x2*x5 - x1*x6 +
195 x3*x6 + x5*(x1 - x2 + x3 + x4 - x5 + x6) in
196 let del3 = x1*x4 - x2*x4 - x1*x5 + x2*x5 -
197 x3*x6 + (x1 + x2 - x3 + x4 + x5 - x6)*x6 in
198 let Pdel = del1 * chain0 in
199 let Ldel = Pdel/del in
200 let sd4 = &4 *x1*del in
201 let sd5 = &4 *x2*del in
202 let sd6 = &4 *x3*del in
203 let Dsd4 = &4 *del + &4 *x1*del1 in
204 let m4diff = &2 *Dn4*sd4 - n4* Dsd4 in
205 let m4 = m4diff*chain0*u234*y2*y3 in
206 let m5 = -- &4 *x2*u234*del3* &2 *x1*u135*y3 in
207 let m6 = -- &4 *x3*u234*del2* &2 *x1*u126*y2 in
208 let const1 = (sol_y (&2) (&2) (&2) (&2) (&2) (&2)) /pi in
209 let rhoy1 = rho(y1) in
210 let rhoy2 = rho(y2) in
211 let rhoy3 = rho(y3) in
212 let Prhoy1 = const1 /(#0.52) in
213 let rr = rhoy1 * m4 + rhoy2 * m5 + rhoy3 * m6 in
215 let D2sd4 = -- &8*x1*x4 + &8*(--(x1*x4) + x2*x5 - x3*x5 - x2*x6 + x3*x6 +
216 x4*(-- x1 + x2 + x3 - x4 + x5 + x6)) in
217 let Dm4diff = &2 * D2n4 * sd4 + Dn4 * Dsd4 - n4 *D2sd4 in
218 let Pm4 = (Dm4diff * chain2 + m4diff * Pchain ) * u234 * y2 * y3 in
219 let Ddel3 = x4 - x5 + x6 in
220 let Ddel2 = x4 + x5 - x6 in
221 let Pm5 = (Ddel3 * x1 * u135 + del3 * &1 * u135 + del3 * x1 * du135) *
222 chain0 * (-- &4 * x2 * u234 * &2 * y3) in
223 let Pm6 = (Ddel2 * x1 * u126 + del2 * &1 * u126 + del2 * x1 * du126) *
224 chain0 * (-- &4 * x3 * u234 * &2 * y2) in
225 let PrrC = &2 * Prhoy1 * m4 + rhoy1 * Pm4 + rhoy2 * Pm5 + rhoy3 * Pm6 in
226 let P2tauNum = (PrrC) + (-- Luf - #0.5 * Ldel) * rr in
227 let P2tau_uf = P2tauNum/ (safesqrt(del)) in
231 let mdtau_y_LC = new_definition `mdtau_y_LC = mdtau_y`;;
233 (* let mdtau2_y_LC = new_definition `mdtau2_y_LC = mdtau2_y`;; *)
235 let mdtau2uf_y_LC = new_definition `mdtau2uf_y_LC = mdtau2uf_y`;;