legacy/oldnonlinear/mdtau.hl

   1 (* ========================================================================== *)
   2 (* FLYSPECK - CODE INFORMAL                                                   *)
   3 (*                                                                            *)
   4 (* Chapter: Nonlinear                                                         *)
   5 (* Author: Thomas C. Hales                                                    *)
   6 (* Date: 2012-04-15                                                           *)
   7 (* ========================================================================== *)
   8
   9
  10 (*
  11 Define the functions mdtau and mdtau2.
  12
  13 In basic terms,
  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.
  16
  17 mdtau2 = D[taumar,{y1,2}].
  18
  19 mdtau2uf = mdtau2 * uf, where uf is the positive factor defined below.
  20
  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.
  23
  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");
  28
  29 *)
  30
  31 module Mdtau = struct
  32
  33 let dua = new_definition `dua a b c = &2 * (b + c - a)`;;
  34
  35 (* corrected June 25, 2012, I had left out the sqrt! *)
  36
  37 let safesqrt = new_definition `safesqrt x = if (x >= &0) then sqrt x else (&0)`;;
  38
  39 let mdtau_y = new_definition `mdtau_y y1 y2 y3 y4 y5 y6 =
  40   let x1 = y1 * y1 in
  41   let x2 = y2 * y2 in
  42   let x3 = y3 * y3 in
  43   let x4 = y4 * y4 in
  44   let x5 = y5 * y5 in
  45   let x6 = y6 * y6 in
  46   let chain0 = &2 * y1 in
  47   let Pchain = &2 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
  58   let del4 = -- n4 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
  89   let t2 = t*t in
  90   let term2a = del * t * matan(t2 *del) in
  91   let term2 = term2a * Prhoy1 in
  92   let term3 = rr/uf in
  93     term1+term2+term3`;;
  94
  95 (*
  96 let mdtau2_y = new_definition `mdtau2_y y1 y2 y3 y4 y5 y6 =
  97   let x1 = y1 * y1 in
  98   let x2 = y2 * y2 in
  99   let x3 = y3 * y3 in
 100   let x4 = y4 * y4 in
 101   let x5 = y5 * y5 in
 102   let x6 = y6 * y6 in
 103   let chain0 = &2 * y1 in
 104   let Pchain = &2 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
 115   let del4 = -- n4 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
 144   let D2n4 = &2  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
 158     P2tau`;;
 159 *)
 160
 161
 162 (*
 163 mdtau2uf = mdtau2 * uf, where uf = &4 * u135 * u126 * u234 * y1 * y2 * y3
 164 *)
 165
 166 let mdtau2uf_y = new_definition `mdtau2uf_y y1 y2 y3 y4 y5 y6 =
 167   let x1 = y1 * y1 in
 168   let x2 = y2 * y2 in
 169   let x3 = y3 * y3 in
 170   let x4 = y4 * y4 in
 171   let x5 = y5 * y5 in
 172   let x6 = y6 * y6 in
 173   let chain0 = &2 * y1 in
 174   let Pchain = &2 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
 185   let del4 = -- n4 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
 214   let D2n4 = &2  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
 228     P2tau_uf`;;
 229
 230
 231   let mdtau_y_LC = new_definition `mdtau_y_LC = mdtau_y`;;
 232
 233 (*  let mdtau2_y_LC = new_definition `mdtau2_y_LC = mdtau2_y`;; *)
 234
 235   let mdtau2uf_y_LC = new_definition `mdtau2uf_y_LC = mdtau2uf_y`;;
 236
 237  end;;