summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--schofield.wl52
1 files changed, 32 insertions, 20 deletions
diff --git a/schofield.wl b/schofield.wl
index 908ca94..77dce96 100644
--- a/schofield.wl
+++ b/schofield.wl
@@ -63,10 +63,10 @@ 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} |->
+ Pns = FoldList[{Pm, m} \[Function]
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]
+ MapIndexed[{Pn, i} \[Function] Pn/dfs[[1]]^(2 i[[1]] - 1), Pns]
]
dFdξLowList[n_, h_][m_] := Module[
@@ -99,7 +99,7 @@ dFdη[n_, h_][m_][tt_] := Module[{ff, hh}, Nest[ddη[hh], h[θ]^(-2 / Δ[]) (ff[
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]
+ruleθ0[g_] := Simplify[g[I θ0]/(-t[I θ0])^Δ[2]/I] - Around[0.18930, 0.00005]
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]
@@ -188,8 +188,8 @@ eq[n_, g_][m_, p_, q_] := Flatten[Join[{ruleθ0[g], ruleAH[g], g'[0] θc - 1}, e
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_, δ_:10^(-15)] := If[Head[#1]===Around,
- (#1["Value"] - #2) / Max[#1["Uncertainty"], δ],
- (#1 - #2) / δ] & @@@ Thread[{data, functions}]
+ (#2-#1["Value"]) / Max[#1["Uncertainty"], δ],
+ (#2-#1) / δ] & @@@ Thread[{data, functions}]
resLow[n_, g_, δ_][m_] := formResiduals[Gls[[;;m+1]], dFdξLowList[n, g][m], δ]
@@ -208,7 +208,16 @@ newSol[eqs_, oldSol_, newVars_, δ_:0, γ_:0, opts___] := FindRoot[
opts
]
-levenbergMarquardtStep[M_, λ_, g_] := Quiet[LinearSolve[M + λ DiagonalMatrix[Diagonal[M]], g], LinearSolve::luc]
+levenbergMarquardtStep[M_, λ_, g_] := Quiet[
+ - LinearSolve[M + λ DiagonalMatrix[Diagonal[M]], g]
+ , LinearSolve::luc
+]
+
+levenbergMarquardtAcceleration[M_, λ_, v_, J_, rold_, rf_, h_, β_] := Quiet[
+ - Re[LinearSolve[M + λ DiagonalMatrix[Diagonal[M]], Transpose[J].(2/h((rf[β[[All,2]] + h v] - rold) / h - J . v))] / 2]
+ , LinearSolve::luc
+]
+
levenburgMarquardt[r_, β0_, λ0_ : 1, ν_ : 2, ε_ : 10^-15] :=
Module[
@@ -217,37 +226,37 @@ levenburgMarquardt[r_, β0_, λ0_ : 1, ν_ : 2, ε_ : 10^-15] :=
β = β0,
λ = λ0,
newβ = β0,
- x, rf, Jf, oldJ, oldr, newr, M, g, δ, oldC, newC,
- h = 0.1, fvv, α, a, v
+ x, rf, Jf, oldJ, oldr, newr, M, g, δ, oldC, newC, δlast,
+ h = 0.1, fvv, α = 0.75, a, v, U, nSteps = 0
},
PrintTemporary["Compiling the Jacobian..."];
rf = Quiet[Compile[{{x, _Real, 1}},
Evaluate[r /. Thread[Rule[First /@ β, Part[x, #] & /@ Range[n]]]]
], Part::partd];
Jf = Quiet[Compile[{{x, _Real, 1}},
- Evaluate[-D[r, {First /@ β}] /. Thread[Rule[First /@ β, Part[x, #] & /@ Range[n]]]]
+ Evaluate[D[r, {First /@ β}] /. Thread[Rule[First /@ β, Part[x, #] & /@ Range[n]]]]
], Part::partd];
PrintTemporary["Beginning the algorithm."];
oldr = rf[β[[All, 2]]];
oldJ = Jf[β[[All, 2]]];
+ U = First[SingularValueDecomposition[oldJ]];
oldC = Re[Total[oldr^2]];
newC = oldC;
g = Re[Transpose[oldJ] . oldr];
M = Re[Transpose[oldJ] . oldJ];
v = levenbergMarquardtStep[M, λ / ν, g];
- a = LinearSolve[oldJ, -2/h((rf[β[[All,2]] + h v] - oldr) / h - oldJ . v)];
- If[2 Norm[a] / Norm[v] > α, a /= (2 / Norm[v])];
- δ = v + a / 2;
- PrintTemporary["Current cost value: ", Dynamic[oldC]]
+ a = levenbergMarquardtAcceleration[M, λ / ν, v, oldJ, oldr, rf, h, β];
+ δ = v + a;
+ δlast = δ;
+ PrintTemporary["Current cost value: ", Dynamic[oldC]];
While[Norm[δ] > ε,
newβ[[All, 2]] += δ;
newr = rf[newβ[[All, 2]]];
newC = Re[Total[newr^2]];
- While[newC > oldC,
+ While[((1 - δ . δlast / Norm[δ] / Norm[δlast])newC > oldC) || (2 Norm[a] / Norm[v] > α),
v = levenbergMarquardtStep[M, λ, g];
- a = LinearSolve[oldJ, -2/h((rf[β[[All,2]] + h v] - oldr) / h - oldJ . v)];
- If[2 Norm[a] / Norm[v] > α, a /= (2 / Norm[v])];
- δ = v + a / 2;
+ a = levenbergMarquardtAcceleration[M, λ, v, oldJ, oldr, rf, h, β];
+ δ = v + a;
newβ = β;
newβ[[All, 2]] += δ;
newr = rf[newβ[[All, 2]]];
@@ -259,13 +268,16 @@ levenburgMarquardt[r_, β0_, λ0_ : 1, ν_ : 2, ε_ : 10^-15] :=
oldr = newr;
β = newβ;
oldJ = Jf[β[[All, 2]]];
+ U = First[SingularValueDecomposition[oldJ]];
g = Re[Transpose[oldJ] . oldr];
M = Re[Transpose[oldJ] . oldJ];
v = levenbergMarquardtStep[M, λ / ν, g];
- a = LinearSolve[oldJ, -2/h((rf[β[[All,2]] + h v] - oldr) / h - oldJ . v)];
- If[2 Norm[a] / Norm[v] > α, a /= (2 / Norm[v])];
- δ = v + a / 2;
+ a = levenbergMarquardtAcceleration[M, λ / ν, v, oldJ, oldr, rf, h, β];
+ δlast = δ;
+ δ = v + a;
+ nSteps += 1;
];
+ Print[nSteps];
{newC, β}
]