Merge branch 'master' of git:research/first_order_singularities/mma

1 files changed, 28 insertions, 36 deletions

diff --git a/schofield.wl b/schofield.wl
index b70b3ac..c60daf5 100644
--- a/schofield.wl
+++ b/schofield.wl
@@ -38,28 +38,18 @@ t[θ_] := ((θ/1)^2 - 1)
 h[n_][θ_] := (1 - (θ/θc)^2) Sum[gC[i] θ^(2*i+1), {i, 0, n}]
 η[g_][θ_] := t[θ] / (g[θ] / I)^(1 / Δ[2])
 
-RFLow[B_, θc_][θ_] := (1/\[Pi])(2 E^(1/(
-   B \[Theta]c)) \[Theta]c ExpIntegralEi[-(1/(B \[Theta]c))] + 
-  E^(1/(B (-\[Theta] + \[Theta]c))) (\[Theta] - \[Theta]c) \
-ExpIntegralEi[1/(B \[Theta] - B \[Theta]c)] - 
-  E^(1/(B \[Theta] + 
-    B \[Theta]c)) (\[Theta] + \[Theta]c) ExpIntegralEi[-(1/(
-     B \[Theta] + B \[Theta]c))])
+fLow[B_, θc_][θ_] := (θc Exp[1/(B θc)] ExpIntegralEi[-1/(B θc)] + (θ - θc) Exp[-1/(B (θ - θc))] ExpIntegralEi[1/(B (θ - θc))]) / π
+
+RFLow[B_, θc_][θ_] := fLow[B, θc][θ] + fLow[B, θc][-θ]
 RFHigh[ξ0_][ξ_] := (ξ^2+ξ0^2)^(5/6) - (ξ0^2)^(5/6)
 
-RF[n_][θ_] := AL RFLow[B, θc][θ] + AH RFHigh[θ0][θ] + Sum[A[i] θ^(2 i), {i, 1, n}]
-RFReg[n_][θ_] := AL (1/\[Pi])(2 E^(1/(
-   B \[Theta]c)) \[Theta]c ExpIntegralEi[-(1/(B \[Theta]c))] - 
-  E^(1/(B \[Theta] + 
-    B \[Theta]c)) (\[Theta] + \[Theta]c) ExpIntegralEi[-(1/(
-     B \[Theta] + B \[Theta]c))]) + AH RFHigh[θ0][θ] + Sum[A[i] θ^(2 i), {i, 1, n}]
-dRFc[n_][m_] := Piecewise[{{AL m! Gamma[m - 1] B^(m - 1) / π, m>1}, {0, True}}] + D[RFReg[n][θ], {θ, m}] /. θ -> θc
+RFReg[n_][θ_] := AL fLow[B, θc][-θ] + AH RFHigh[θ0][θ] + Sum[A[i] θ^(2 i), {i, 1, n}]
+RF[n_][θ_] := RFReg[n][θ] + AL fLow[B, θc][θ]
+dfLow[AL_, B_][m_] := Piecewise[{{AL m! Gamma[m - 1] B^(m - 1) / π, m > 1}, {AL θc Exp[1/(B θc)] ExpIntegralEi[-1/(B θc)] / π, m == 0}, {0, True}}]
+dRFc[n_][m_] := dfLow[AL, B][m] + D[RFReg[n][θ], {θ, m}] /. θ -> θc
 
 RFC[n_][θ_] := RF[n][θ] + AL I Sign[Im[θ]] ((θ-θc)Exp[-1/(B(θ-θc))]-(-θ-θc)Exp[-1/(B(-θ-θc))])
 
-ddξ[h_][f_] := D[f, θ] / D[h[θ] / RealAbs[t[θ]]^Δ[2], θ]
-ddη[h_][f_] := D[f, θ] / D[t[θ] / h[θ]^(1 / Δ[2]), θ]
-
 invDerivativeList[n_][f_][x_] := Module[
     {xp, dfs, fp, Pns},
   dfs = Rest[NestList[D[#, xp] &, f[xp], n]] /. xp -> x;
@@ -69,10 +59,10 @@ invDerivativeList[n_][f_][x_] := Module[
   MapIndexed[{Pn, i} \[Function] Pn/dfs[[1]]^(2 i[[1]] - 1), Pns]
   ]
 
-dFdξLowList[n_, h_][m_] := Module[
+dGdξLowList[n_, h_][m_] := Module[
   { ds, dF, df },
   ds = invDerivativeList[m+1][Function[θ, h[θ] / t[θ]^Δ[2]]][θc];
-  dF = NestList[Function[f, D[f, θ]], RFReg[n][θ], m] + Table[Piecewise[{{AL k! Gamma[k - 1] B^(k - 1)/\[Pi], k > 1}, {0, True}}], {k, 0, m}] /. θ -> θc;
+  dF = NestList[Function[f, D[f, θ]], RFReg[n][θ], m] + Table[dfLow[AL, B][k], {k, 0, m}] /. θ -> θc;
   df = NestList[D[#, \[Theta]] &, 
       fp[\[Theta]]/t[\[Theta]]^2 - 1/(8 \[Pi]) Log[t[\[Theta]]^2], 
       m] /. 
@@ -81,22 +71,29 @@ dFdξLowList[n_, h_][m_] := Module[
   Table[Sum[df[[k+1]] BellY[j, k, ds[[;; j - k + 1]]], {k, 0, j}]/(j!), {j, 0, m}]
 ]
 
-dFdξHighList[n_, h_][m_] := Module[
+dGdξList[n_, h_][m_, θp_] := Module[
   { ds, dF, df },
-  ds = invDerivativeList[m+1][Function[θ, h[θ] / (-t[θ])^Δ[2]]][0];
-  dF = NestList[Function[f, D[f, θ]], RF[n][θ], m] /. θ -> 0;
+  ds = invDerivativeList[m+1][Function[θ, h[θ] / RealAbs[t[θ]]^Δ[2]]][θp];
+  dF = NestList[Function[f, D[f, θ]], RF[n][θ], m] /. θ -> θp;
   df = NestList[D[#, \[Theta]] &, 
       fp[\[Theta]]/t[\[Theta]]^2 - 1/(8 \[Pi]) Log[t[\[Theta]]^2], 
       m] /. 
    Map[Derivative[#][fp][\[Theta]] -> dF[[# + 1]] &, 
-    Range[0, m]] /. θ -> 0;
+    Range[0, m]] /. θ -> θp;
   Table[Sum[df[[k+1]] BellY[j, k, ds[[;; j - k + 1]]], {k, 0, j}]/(j!), {j, 0, m}]
 ]
 
-dFdξLow[n_, h_][m_] := Module[{ff, hh}, Nest[ddξ[hh], ff[θ] / t[θ]^2 - Log[t[θ]^2] / (8 π), m] /. θ -> θc /. Map[Derivative[#][ff][θc] -> dRFc[n][#] &, Range[0, m]] /. hh -> h]
-dFdξHigh[n_, h_][m_] := Module[{ff, hh}, Nest[ddξ[hh], ff[θ] / t[θ]^2 - Log[t[θ]^2] / (8 π), m] /. θ -> 0 /. Map[Derivative[#][ff][0] -> eqHighRHS[RF[n]][#] &, Range[0, m]] /. hh -> h]
-dFdη[n_, h_][m_][tt_] := Module[{ff, hh}, Nest[ddη[hh], h[θ]^(-2 / Δ[]) (ff[θ] - t[θ]^2 Log[hh[θ]^2] / (8 π Δ[])), m] /. θ -> tt /. Map[Derivative[#][ff][tt] -> Derivative[#][RF[n]][tt] &, Range[0, m]] /. hh -> h]
-dFdηList[n_, h_][m_][tt_] := Module[{ff, hh}, NestList[ddη[hh], h[θ]^(-2 / Δ[2]) (ff[θ] - t[θ]^2 Log[hh[θ]^2] / (8 π Δ[2])), m] /. θ -> tt /. Map[Derivative[#][ff][tt] -> Derivative[#][RF[n]][tt] &, Range[0, m]] /. hh -> h]
+dΦdηList[n_, h_][m_, θp_] := Module[
+  { ds, dF, df },
+  ds = invDerivativeList[m+1][Function[θ, t[θ] / h[θ]^(1 / Δ[2])]][θp];
+  dF = NestList[Function[f, D[f, θ]], RF[n][θ], m] /. θ -> θp;
+  df = NestList[D[#, \[Theta]] &, 
+      hh[θ]^(-2 / Δ[]) (fp[θ] - t[θ]^2 Log[hh[θ]^2] / (8 π Δ[])), 
+      m] /. 
+   Map[Derivative[#][fp][\[Theta]] -> dF[[# + 1]] &, 
+    Range[0, m]] /. hh -> h /. θ -> θp;
+  Table[Sum[df[[k+1]] BellY[j, k, ds[[;; j - k + 1]]], {k, 0, j}]/(j!), {j, 0, m}]
+]
 
 ruleB[g_] := B -> (2 * OverBar[s] / π) * (- g'[θc] / t[θc]^Δ[2])
 ruleθ0[g_] := Simplify[g[I θ0]/(-t[I θ0])^Δ[2]/I] - Around[0.18930, 0.00005]
@@ -122,8 +119,6 @@ eqMid[F_, h_][m_] := D[
     /. η -> t[θ] / h[θ]^(1 / Δ[2]),
   {θ, m} ] / m! /. θ -> 1
 
-δ0 = 10^(-14);
-
 Φs = {
   -1.19773338379799339,
   -0.31881012489061,
@@ -180,10 +175,7 @@ dRule[sym_][f_, i_] := Derivative[i[[1]] - 1][sym][0] -> f (i[[1]] - 1)!
 GlRules = MapIndexed[dRule[Gl], Gls];
 GhRules = MapIndexed[dRule[Gh], Ghs];
 
-ClearAll[rules]
-rules[g_] := Join[ΦRules, GlRules, GhRules, {ruleAL[g], ruleB[g], gC[0]->1}]
-
-eq[n_, g_][m_, p_, q_] := Flatten[Join[{ruleθ0[g], ruleAH[g], g'[0] θc - 1}, eqLow[n, g][#] & /@ Range[0, m],eqMid[RF[n], g][#] & /@ Range[0, p], eqHigh[n, g] /@ Range[2, q, 2]]] //. rules[g] /. Around[x_, _] :> x
+rules[g_] := Join[ΦRules, GlRules, GhRules, {ruleAL[g], ruleB[g], gC[0] -> 1}]
 
 eqAround[n_, g_][m_, p_, q_] := Flatten[Join[{ruleAH[g], ruleθ0[g]}, eqLow[n, g][#] & /@ Range[0, m],eqMid[RF[n], g][#] & /@ Range[0, p], eqHigh[n, g] /@ Range[2, q, 2]]] //. rules[g]
 
@@ -191,9 +183,9 @@ formResiduals[data_, functions_, δ_:10^(-15)] := If[Head[#1]===Around,
   (#2-#1["Value"]) / Max[#1["Uncertainty"], δ],
   (#2-#1) / δ] & @@@ Thread[{data, functions}]
 
-resLow[n_, g_, δ_][m_] := formResiduals[Gls[[;;m+1]], dFdξLowList[n, g][m], δ]
+resLow[n_, g_, δ_][m_] := formResiduals[Gls[[;;m+1]], dGdξLowList[n, g][m], δ]
 
-resHigh[n_, g_, δ_][m_] := Rest[formResiduals[Ghs[[;;m+1]], dFdξHighList[n, g][m], δ][[;;;;2]]]
+resHigh[n_, g_, δ_][m_] := Rest[formResiduals[Ghs[[;;m+1]], dGdξList[n, g][m, 0], δ][[;;;;2]]]
 
 res[F_, g_, δ_][m_] := Join[resLow[F, g, δ][m], resHigh[F, g, δ][m], {ruleθ0[g] / 0.00005, ruleAH[g] / 0.02} /. Around[x_, _] :> x]
 chiSquared[F_, g_, δ_][m_] := Total[res[F, g, δ][m]^2]
@@ -245,7 +237,7 @@ levenbergMarquardtHelper[Δ_, rf_, Jf_, β0_, λ0_ : 1, ν_ : 2, ε_ : 10^-15, a
     λ = λ0,
     newβ = β0,
     x, oldJ, oldr, newr, M, g, δ, oldC, newC, vlast,
-    h = 0.1, fvv, α = 0.75, a, v, U, nSteps = 0, reject, μ =3
+    h = 0.1, fvv, α = 0.1, a, v, U, nSteps = 0, reject, μ =3
    },
   PrintTemporary["Beginning the algorithm."];
   oldr = rf[Append[β[[All, 2]], Δ]];