diff options
Diffstat (limited to 'schofield.wl')
| -rw-r--r-- | schofield.wl | 64 |
1 files changed, 39 insertions, 25 deletions
diff --git a/schofield.wl b/schofield.wl index 1c9482c..b26b662 100644 --- a/schofield.wl +++ b/schofield.wl @@ -30,27 +30,36 @@ $Assumptions = {θc > 0, θc > 1, gC[_] ∈ Reals, B > 0, γ > 0, ξ0 > 0} } ] +OverBar[s] := 1.357838341706595496 + Δ[D_:2] := β[D] δ[D] t[θ_] := ((θ)^2 - 1) + (* hBasis = Orthogonalize[x^# & /@ Range[0, 20], Function[{f, g}, Integrate[f g (1 - x^2)^2, {x, -1, 1}]]] + h[n_][θ_] := (1 - (θ/θc)^2) Sum[gC[i] hBasis[[2 * i + 2]], {i, 0, n}] /. x -> θ / θc -*) h[n_][θ_] := (1 - (θ/θc)^2) Sum[gC[i]θ^(2 * i + 1), {i, 0, n}] +*) +h[n_][θ_] := (1 - (θ/θc)^2) Sum[gC[i]LegendreP[(2 * i + 1), θ/θc], {i, 0, n}] -RFLow[B_, θc_][θ_] := (2 E^(1/( +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))]) -RFHigh[ξ0_][ξ_] := (1+(ξ/ξ0)^2)^(5/6) - 1 +RFHigh[ξ0_][ξ_] := (ξ^2+ξ0^2)^(5/6) -RF[n_][θ_] := AL RFLow[B, θc][θ] + AH RFHigh[θ0][θ] + Sum[A[i] θ^(2(i+1)), {i, 0, n}] +RF[n_][θ_] := AL RFLow[B, θc][θ] + AH RFHigh[θ0][θ] + Sum[A[i] LegendreP[2i, θ], {i, 0, n}] -ruleB[g_] := B - π / (2 1.3578383417066) (- g'[θc] / t[θc]^Δ[2]) +RFC[n_][θ_] := RF[n][θ] + AL I Sign[Im[θ]] ((θ-θc)Exp[-1/(B(θ-θc))]-(-θ-θc)Exp[-1/(B(-θ-θc))]) + +ruleB[g_] := B - (2 * OverBar[s] / π) * (- g'[θc] / t[θc]^Δ[2]) ruleθ0[g_] := Simplify[g[I θ0]/(-t[I θ0])^Δ[2]/I] - 0.18930 +ruleAL[g_] := AL + t[θc]^2 OverBar[s] / (2 π) * (- g'[θc] / t[θc]^Δ[2]) +ruleg0[g_] := 1 - g'[0] eqLowRHS[F_][m_] := SeriesCoefficient[F[θ], {θ, θc, m}, Assumptions -> Join[{θ < θc, θ > 1}, $Assumptions]] @@ -71,6 +80,8 @@ eqMid[F_, h_][m_] := D[ /. η -> t[θ] / h[θ]^(1 / Δ[2]), {θ, m} ] /. θ -> 1 +δ0 = 10^-16; + Φs = { -1.197733383797993, -0.318810124891, @@ -83,15 +94,15 @@ eqMid[F_, h_][m_] := D[ } Gls = { - 0, - -1.3578383417066, - -0.048953289720, - 0.038863932, - -0.068362119, - 0.18388370, - -0.6591714, - 2.937665, - -15.61, + Around[0, δ0], + Around[-OverBar[s], δ0], + Around[−0.048953289720, 2 10^(-12)], + Around[ 0.0388639290, 1 10^(-10)], + Around[-0.068362121, 1 10^(-9)], + Around[ 0.18388371, 1 10^(-8)], + Around[-0.659170, 1 10^(-6)], + Around[ 2.937665, 3 10^(-6)], + Around[-15.61, 10^(-2)], 96.76, -679, 5340 @@ -100,15 +111,15 @@ Gls = { Ghs = { 0, 0, - -1.8452280782328, + Around[ -1.84522807823, 10^(-11)], 0, - 8.333711750, + Around[ 8.3337117508, 10^(-10)], 0, - -95.16896, + Around[-95.16897, 10^(-5)], 0, - 1457.62, + Around[1457.62, 3 10^(-2)], 0, - -25891, + Around[-25891, 2], 0, 502000 } @@ -119,20 +130,23 @@ dRule[sym_][f_, i_] := Derivative[i[[1]] - 1][sym][0] -> f (i[[1]] - 1)! GlRules = MapIndexed[dRule[Gl], Gls]; GhRules = MapIndexed[dRule[Gh], Ghs]; -rules[g_] := Join[ΦRules, GlRules, GhRules] +rules := Join[ΦRules, GlRules, GhRules] (*ξ0 := 0.18930*) -gC[0] := 1 +(*gC[0] := 1*) tC[0] := 1 -eq[F_, g_][m_] := Flatten[Join[{ruleB[g], ruleθ0[g]},{eqLow[F, g][#](*, eqMid[F, g][#]*)} & /@ Range[0, m], eqHigh[F, g] /@ Range[2, m, 2]]] //. rules[g] +eq[F_, g_][m_] := Flatten[Join[{ruleg0[g], ruleB[g], ruleθ0[g]},{eqLow[F, g][#](*, eqMid[F, g][#]*)} & /@ Range[0, m], eqHigh[F, g] /@ Range[0, m, 2]]] //. rules /. Around[x_, _] :> x + +chiSquaredLow[F_, g_][m_] := Total[(((#[[1]] /. rules)["Value"] - #[[2]])^2 / (#[[1]] /. rules)["Uncertainty"]^2)& /@ Solve[0 == (eqLow[ff, g] /@ Range[0, m]), Derivative[#][Gl][0]& /@ Range[0, m]][[1]]] /. ff -> F +chiSquaredHigh[F_, g_][m_] := Total[(((#[[1]] /. rules)["Value"] - #[[2]])^2 / (#[[1]] /. rules)["Uncertainty"]^2)& /@ Solve[0 == (eqHigh[ff, g] /@ Range[2, m, 2]) /. ff'[0] -> 0, Derivative[#][Gh][0]& /@ Range[2, m, 2]][[1]]] /. ff -> F -newSol[eqs_, oldSol_, newVars_, δ_:0, opts___] := FindRoot[ +newSol[eqs_, oldSol_, newVars_, δ_:0, γ_:0, opts___] := FindRoot[ eqs, Join[ - oldSol /. Rule -> List, + {#1, #2 + γ * RandomVariate[NormalDistribution[]]} & @@@ (oldSol /. Rule -> List), Thread[{newVars, δ * RandomVariate[NormalDistribution[], Length[newVars]]}] ], - MaxIterations -> 500, + MaxIterations -> 1000, opts ] |