(* ========================================================================== *)
(* FLYSPECK - CODE INFORMAL *)
(* *)
(* Chapter: Nonlinear *)
(* Author: Thomas C. Hales *)
(* Date: 2012-04-15 *)
(* ========================================================================== *)
(*
Define the functions mdtau and mdtau2.
In basic terms,
mdtau = Sqrt[delta] * D[taumar,y1]
However, we need to rework the formula in such a way that it extends to delta=0.
mdtau2 = D[taumar,{y1,2}].
mdtau2uf = mdtau2 * uf, where uf is the positive factor defined below.
In /interval_code, taylorInterval.cc, wide.cc, there is a version mdtau_y_LC, mdtau2_y_LC
The "LC" signifies "locally constant", meaning no derivatives are to be used in interval computations.
("mdtau_y_LC","Lib::mdtau_y_LC");
("mdtau2uf_y_LC","Lib::mdtau2uf_y_LC");
("mdtau_y","Lib::mdtau_y_LC");
("mdtau2uf_y","Lib::mdtau2uf_y_LC");
*)
module Mdtau = struct
(* corrected June 25, 2012, I had left out the sqrt! *)
let mdtau_y = new_definition `mdtau_y y1 y2 y3 y4 y5 y6 =
let x1 = y1 * y1 in
let x2 = y2 * y2 in
let x3 = y3 * y3 in
let x4 = y4 * y4 in
let x5 = y5 * y5 in
let x6 = y6 * y6 in
let chain0 = &2 * y1 in
let Pchain = &2 in
let chain2 = &4 * x1 in
let u135 = ups_x x1 x3 x5 in
let u126 = ups_x x1 x2 x6 in
let u234 = ups_x x2 x3 x4 in
let uf = &4 * u135 * u126 * u234 * y1 * y2 * y3 in
let du135 = dua x1 x3 x5 in
let du126 = dua x1 x2 x6 in
let Luf = (du135 / u135 + du126/u126 ) * chain0 + &1 / y1 in
let n4 = x2*x3 + x1*x4 - x2*x5 - x3*x6 + x5*x6 -
x1*(-- x1 + x2 + x3 - x4 + x5 + x6) in
let del4 = -- n4 in
let n5 = x1*x3 - x1*x4 + x2*x5 - x3*x6 + x4*x6 -
x2*(x1 - x2 + x3 + x4 - x5 + x6) in // - del5
let n6 = x1*x2 - x1*x4 - x2*x5 + x4*x5 -
x3*(x1 + x2 - x3 + x4 + x5 - x6) + x3*x6 in // - del6
let Dn4 = &2* x1 - x2 - x3 + &2 *x4 - x5 - x6 in
let del = delta_x x1 x2 x3 x4 x5 x6 in
let del1 = -- (x1*x4) + x2*x5 - x3*x5 - x2*x6 + x3*x6 +
x4*(-- x1 + x2 + x3 - x4 + x5 + x6) in
let del2 = x1*x4 - x3*x4 - x2*x5 - x1*x6 +
x3*x6 + x5*(x1 - x2 + x3 + x4 - x5 + x6) in
let del3 = x1*x4 - x2*x4 - x1*x5 + x2*x5 -
x3*x6 + (x1 + x2 - x3 + x4 + x5 - x6)*x6 in
let Pdel = del1 * chain0 in
let Ldel = Pdel/del in
let sd4 = &4 *x1*del in
let sd5 = &4 *x2*del in
let sd6 = &4 *x3*del in
let Dsd4 = &4 *del + &4 *x1*del1 in
let m4diff = &2 *Dn4*sd4 - n4* Dsd4 in
let m4 = m4diff*chain0*u234*y2*y3 in
let m5 = -- &4 *x2*u234*del3* &2 *x1*u135*y3 in
let m6 = -- &4 *x3*u234*del2* &2 *x1*u126*y2 in
let const1 = (sol_y (&2) (&2) (&2) (&2) (&2) (&2)) /pi in
let rhoy1 = rho(y1) in
let rhoy2 = rho(y2) in
let rhoy3 = rho(y3) in
let Prhoy1 = const1 /(#0.52) in
let rr = rhoy1 * m4 + rhoy2 * m5 + rhoy3 * m6 in
let term1 = Prhoy1 * pi * safesqrt(del) in
let t = sqrt(&4 * x1)/del4 in
let t2 = t*t in
let term2a = del * t * matan(t2 *del) in
let term2 = term2a * Prhoy1 in
let term3 = rr/uf in
term1+term2+term3`;;let mdtau2uf_y = new_definition `mdtau2uf_y y1 y2 y3 y4 y5 y6 =
let x1 = y1 * y1 in
let x2 = y2 * y2 in
let x3 = y3 * y3 in
let x4 = y4 * y4 in
let x5 = y5 * y5 in
let x6 = y6 * y6 in
let chain0 = &2 * y1 in
let Pchain = &2 in
let chain2 = &4 * x1 in
let u135 = ups_x x1 x3 x5 in
let u126 = ups_x x1 x2 x6 in
let u234 = ups_x x2 x3 x4 in
let uf = &4 * u135 * u126 * u234 * y1 * y2 * y3 in
let du135 = dua x1 x3 x5 in
let du126 = dua x1 x2 x6 in
let Luf = (du135 / u135 + du126/u126 ) * chain0 + &1 / y1 in
let n4 = x2*x3 + x1*x4 - x2*x5 - x3*x6 + x5*x6 -
x1*(-- x1 + x2 + x3 - x4 + x5 + x6) in
let del4 = -- n4 in
let n5 = x1*x3 - x1*x4 + x2*x5 - x3*x6 + x4*x6 -
x2*(x1 - x2 + x3 + x4 - x5 + x6) in // - del5
let n6 = x1*x2 - x1*x4 - x2*x5 + x4*x5 -
x3*(x1 + x2 - x3 + x4 + x5 - x6) + x3*x6 in // - del6
let Dn4 = &2* x1 - x2 - x3 + &2 *x4 - x5 - x6 in
let del = delta_x x1 x2 x3 x4 x5 x6 in
let del1 = -- (x1*x4) + x2*x5 - x3*x5 - x2*x6 + x3*x6 +
x4*(-- x1 + x2 + x3 - x4 + x5 + x6) in
let del2 = x1*x4 - x3*x4 - x2*x5 - x1*x6 +
x3*x6 + x5*(x1 - x2 + x3 + x4 - x5 + x6) in
let del3 = x1*x4 - x2*x4 - x1*x5 + x2*x5 -
x3*x6 + (x1 + x2 - x3 + x4 + x5 - x6)*x6 in
let Pdel = del1 * chain0 in
let Ldel = Pdel/del in
let sd4 = &4 *x1*del in
let sd5 = &4 *x2*del in
let sd6 = &4 *x3*del in
let Dsd4 = &4 *del + &4 *x1*del1 in
let m4diff = &2 *Dn4*sd4 - n4* Dsd4 in
let m4 = m4diff*chain0*u234*y2*y3 in
let m5 = -- &4 *x2*u234*del3* &2 *x1*u135*y3 in
let m6 = -- &4 *x3*u234*del2* &2 *x1*u126*y2 in
let const1 = (sol_y (&2) (&2) (&2) (&2) (&2) (&2)) /pi in
let rhoy1 = rho(y1) in
let rhoy2 = rho(y2) in
let rhoy3 = rho(y3) in
let Prhoy1 = const1 /(#0.52) in
let rr = rhoy1 * m4 + rhoy2 * m5 + rhoy3 * m6 in
let D2n4 = &2 in
let D2sd4 = -- &8*x1*x4 + &8*(--(x1*x4) + x2*x5 - x3*x5 - x2*x6 + x3*x6 +
x4*(-- x1 + x2 + x3 - x4 + x5 + x6)) in
let Dm4diff = &2 * D2n4 * sd4 + Dn4 * Dsd4 - n4 *D2sd4 in
let Pm4 = (Dm4diff * chain2 + m4diff * Pchain ) * u234 * y2 * y3 in
let Ddel3 = x4 - x5 + x6 in
let Ddel2 = x4 + x5 - x6 in
let Pm5 = (Ddel3 * x1 * u135 + del3 * &1 * u135 + del3 * x1 * du135) *
chain0 * (-- &4 * x2 * u234 * &2 * y3) in
let Pm6 = (Ddel2 * x1 * u126 + del2 * &1 * u126 + del2 * x1 * du126) *
chain0 * (-- &4 * x3 * u234 * &2 * y2) in
let PrrC = &2 * Prhoy1 * m4 + rhoy1 * Pm4 + rhoy2 * Pm5 + rhoy3 * Pm6 in
let P2tauNum = (PrrC) + (-- Luf - #0.5 * Ldel) * rr in
let P2tau_uf = P2tauNum/ (safesqrt(del)) in
P2tau_uf`;;