1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
|
BeginPackage["Schofield`"]
$Assumptions = {θc > 0, θc > 1, gC[_] ∈ Reals, B > 0}
β[D_:2] := Piecewise[
{
{1/8, D == 2},
{0.326419, D == 3},
{1/2, D == 4},
{β, True}
}
]
δ[D_:2] := Piecewise[
{
{15, D == 2},
{4.78984, D == 3},
{3, D == 4},
{δ, True}
}
]
α[D_:2] := Piecewise[
{
{0, D == 2},
{0.11008, D == 3},
{0, D == 4},
{α, True}
}
]
Δ[D_:2] := β[D] δ[D]
OverBar[s] := 2^(1/12) Exp[-1/8] Glaisher^(3/2)
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))])
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
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;
Pns = FoldList[{Pm, m} |->
fp'[xp] D[Pm, xp] - (2 m - 1) fp''[xp] Pm, 1, Range[n - 1]] /.
Derivative[m_][fp][xp] :> dfs[[m]];
MapIndexed[{Pn, i} |-> Pn/dfs[[1]]^(2 i[[1]] - 1), Pns]
]
dFdξ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[D[#, \[Theta]] &,
fp[\[Theta]]/t[\[Theta]]^2 - 1/(8 \[Pi]) Log[t[\[Theta]]^2],
m] /.
Map[Derivative[#][fp][\[Theta]] -> dF[[# + 1]] &,
Range[0, m]] /. θ -> θc;
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[
{ ds, dF, df },
ds = invDerivativeList[m+1][Function[θ, h[θ] / (-t[θ])^Δ[2]]][0];
dF = NestList[Function[f, D[f, θ]], RF[n][θ], m] /. θ -> 0;
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;
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]
ruleB[g_] := B -> (2 * OverBar[s] / π) * (- g'[θc] / t[θc]^Δ[2])
ruleθ0[g_] := Around[0.18930, 0.00005] - Simplify[g[I θ0]/(-t[I θ0])^Δ[2]/I]
ruleAL[g_] := AL -> Exp[Δ[2] t[θc]^(Δ[2] - 1) t'[θc] / (2 OverBar[s] / π g'[θc]) - t[θc]^Δ[2] g''[θc] / (4 OverBar[s] / π g'[θc]^2)] t[θc]^(1/8) OverBar[s] / (2 π) * g'[θc]
ruleAH[g_] := AH / ((g[I θ0]/ I)^(2 / Δ[2]) * (-η[g]'[I θ0] / (2 θ0 I))^(5/6)) + Around[1.37, 0.02]
eqLowRHSReg[n_][m_] := dRFc[n][m]
eqLowLHS[h_][m_] :=D[
t[θ]^2 (Gl[h[θ] t[θ]^-Δ[2]] + Log[t[θ]^2]/(8 π)),
{θ, m} ] /. θ -> θc
eqLow[n_, h_][m_] := (eqLowRHSReg[n][m] - eqLowLHS[h][m]) / m!
eqHighRHS[F_][m_] := D[F[θ], {θ, m} ] /. θ -> 0
eqHighLHS[h_][m_] := D[(-t[θ])^2 (Gh[h[θ] (-t[θ])^-Δ[2]] + Log[(-t[θ])^2]/(8 π)), {θ, m} ] /. θ -> 0
eqHigh[n_, h_][m_] := (eqHighRHS[RF[n]][m] - eqHighLHS[h][m]) / m!
eqMid[F_, h_][m_] := D[
F[θ] - t[θ]^2 Log[h[θ]^2]/(8 Δ[2]π) - h[θ]^((2-α[2])/Δ[2]) Φ[η]
/. η -> t[θ] / h[θ]^(1 / Δ[2]),
{θ, m} ] / m! /. θ -> 1
δ0 = 10^(-14);
Φs = {
-1.19773338379799339,
-0.31881012489061,
0.110886196683,
0.01642689465,
-2.639978 10^-4,
-5.140526 10^-4,
2.08856 10^-4,
-4.4819 10^-5,
3.16 10^-7,
4.31 10^-6,
-1.99 10^-6
}
Gls = {
0,
-OverBar[s],
−1.000960328725262189480934955172097320572505951770117 Sqrt[2]/((2 )^(-7/8) (2^(3/16)/OverBar[s])^2)/2/(12 \[Pi]),
Around[ 0.038863932, 3.0 10^(-9)],
Around[−0.068362119, 2.0 10^(-9)],
Around[ 0.18388371, 1.0 10^(-8)],
Around[-0.659170, 1.0 10^(-6)],
Around[ 2.937665, 3.0 10^(-6)],
Around[-15.61, 1.0 10^(-2)],
Around[ 96.76, 1.0 10^(-2)],
-6.79 10^2,
5.34 10^3,
-4.66 10^4,
4.46 10^5,
-4.66 10^6
}
Ghs = {
0,
0,
-1.000815260440212647119476363047210236937534925597789 Sqrt[2]/((2 )^(-7/8) (2^(3/16)/OverBar[s])^2)/2,
0,
Around[ 8.333711750, 5.0 10^(-9)],
0,
Around[-95.16897, 3.0 10^(-5)],
0,
Around[1457.62, 3.0 10^(-2)],
0,
Around[-2.5891 10^4, 2.0],
0,
5.02 10^5,
0,
-1.04 10^7
}
dRule[sym_][f_, i_] := Derivative[i[[1]] - 1][sym][0] -> f (i[[1]] - 1)!
ΦRules = MapIndexed[dRule[Φ], Φs];
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
eqAround[n_, g_][m_, p_, q_] := Flatten[Join[{ruleθ0[g], ruleAH[g]}, eqLow[n, g][#] & /@ Range[0, m],eqMid[RF[n], g][#] & /@ Range[0, p], eqHigh[n, g] /@ Range[2, q, 2]]] //. rules[g]
formResiduals[data_, functions_, δ_:0] := If[Head[#1]===Around,
(#1["Value"] - #2) / Max[#1["Uncertainty"], δ],
(#1 - #2) / δ] & @@@ Thread[{data, functions}]
resLow[n_, g_, δ_][m_] := formResiduals[Gls[[;;m+1]], dFdξLowList[n, g][m], δ]
resHigh[n_, g_, δ_][m_] := Rest[formResiduals[Ghs[[;;m+1]], dFdξHighList[n, g][m], δ][[;;;;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]
newSol[eqs_, oldSol_, newVars_, δ_:0, γ_:0, opts___] := FindRoot[
eqs,
Join[
{#1, #2 + γ * RandomVariate[NormalDistribution[]]} & @@@ (oldSol /. Rule -> List),
Thread[{newVars, δ * RandomVariate[NormalDistribution[], Length[newVars]]}]
],
MaxIterations -> 500,
opts
]
levenburgMarquardt[r_, \[Beta]0_, \[Lambda]0_ : 1, \[Nu]_ : 2] :=
Module[
{
n, \[Beta], new\[Beta], rc, J, I, oldJ, newJ, oldr, newr, M,
g, \[Delta], newcost, \[Lambda], oldcost, temp
},
n = Length[\[Beta]0]; \[Lambda] = \[Lambda]0; \[Beta] = \[Beta]0;
PrintTemporary["Compiling the Jacobian..."];
rc = Compile[{{x, _Real, 1}},
Evaluate[
r /. Thread[Rule[First /@ \[Beta], Part[x, #] & /@ Range[n]]]]];
J = Compile[{{x, _Real, 1}},
Evaluate[
D[r, {First /@ \[Beta]}] /.
Thread[Rule[First /@ \[Beta], Part[x, #] & /@ Range[n]]]]];
PrintTemporary["Beginning the algorithm."];
oldJ = J[\[Beta][[All, 2]]];
oldr = rc[\[Beta][[All, 2]]];
oldcost = Re[Total[oldr^2]];
new\[Beta] = \[Beta];
g = Re[Transpose[oldJ] . oldr];
M = Re[Transpose[oldJ] . oldJ];
\[Delta] = LinearSolve[M + \[Lambda] DiagonalMatrix[Diagonal[M]], g];
PrintTemporary[Dynamic[oldcost]]
While[Norm[\[Delta]] > 10^-15,
\[Delta] =
LinearSolve[M + \[Lambda]/\[Nu] DiagonalMatrix[Diagonal[M]], g];
new\[Beta][[All, 2]] -= \[Delta];
newr = rc[new\[Beta][[All, 2]]];
newcost = Re[Total[newr^2]];
While[newcost > oldcost,
\[Delta] =
LinearSolve[M + \[Lambda] DiagonalMatrix[Diagonal[M]], g];
new\[Beta] = \[Beta];
new\[Beta][[All, 2]] -= \[Delta];
newr = rc[new\[Beta][[All, 2]]];
newcost = Re[Total[newr^2]];
\[Lambda] *= \[Nu];
];
\[Lambda] /= \[Nu];
oldcost = newcost;
oldr = newr;
\[Beta] = new\[Beta];
oldJ = J[\[Beta][[All, 2]]];
g = Re[Transpose[oldJ] . oldr];
M = Re[Transpose[oldJ] . oldJ];
];
{newcost, \[Beta]}
]
EndPackage[]
|
