diff options
| -rw-r--r-- | schofield.wl | 52 |
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, β} ] |