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
|
BeginPackage["Schofield`"]
$Assumptions = {θc > 0, θc > 1, gC[_] ∈ Reals, B > 0, γ > 0, ξ0 > 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]
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}]
RFLow[B_, θc_][θ_] := (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
RF[n_][θ_] := AL RFLow[B, θc][θ] + AH RFHigh[θ0][θ] + Sum[A[i] θ^(2(i+1)), {i, 0, n}]
ruleB[g_] := B - π / (2 1.3578383417066) (- g'[θc] / t[θc]^Δ[2])
ruleθ0[g_] := Simplify[g[I θ0]/(-t[I θ0])^Δ[2]/I] - 0.18930
eqLowRHS[F_][m_] := SeriesCoefficient[F[θ], {θ, θc, m}, Assumptions -> Join[{θ < θc, θ > 1}, $Assumptions]]
eqLowLHS[h_][m_] :=D[
t[θ]^2 (Gl[h[θ] t[θ]^-Δ[2]] + Log[t[θ]^2]/(8 π)),
{θ, m} ] / m! /. θ -> θc
eqLow[F_, h_][m_] := eqLowRHS[F][m] - eqLowLHS[h][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[F_, h_][m_] := eqHighRHS[F][m] - eqHighLHS[h][m]
eqMid[F_, h_][m_] := D[
F[θ] - t[θ]^2 Log[h[θ]^2]/(8 Δ[2]π) - h[θ]^((2-α[2])/Δ[2]) Φ[η]
/. η -> t[θ] / h[θ]^(1 / Δ[2]),
{θ, m} ] /. θ -> 1
Φs = {
-1.197733383797993,
-0.318810124891,
0.110886196683,
0.01642689465,
-2.639978 10^-4,
-5.140526 10^-4,
2.08856 10^-4,
-4.4819 10^-5
}
Gls = {
0,
-1.3578383417066,
-0.048953289720,
0.038863932,
-0.068362119,
0.18388370,
-0.6591714,
2.937665,
-15.61,
96.76,
-679,
5340
}
Ghs = {
0,
0,
-1.8452280782328,
0,
8.333711750,
0,
-95.16896,
0,
1457.62,
0,
-25891,
0,
502000
}
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];
rules[g_] := Join[ΦRules, GlRules, GhRules]
(*ξ0 := 0.18930*)
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]
newSol[eqs_, oldSol_, newVars_, δ_:0, opts___] := FindRoot[
eqs,
Join[
oldSol /. Rule -> List,
Thread[{newVars, δ * RandomVariate[NormalDistribution[], Length[newVars]]}]
],
MaxIterations -> 500,
opts
]
EndPackage[]
|
