diff --git a/euler-newton.jl b/euler-newton.jl index e3cefd3..7ec2c04 100644 --- a/euler-newton.jl +++ b/euler-newton.jl @@ -11,19 +11,21 @@ module EulerNewton vars = variables(H(t)) # Jacobian of H evaluated at (x,t) JH = [jh(vars=>x) for jh in differentiate(H(t), vars)] - Δx = JH \ -[gg(vars=>x) for gg in H(1)-H(0)] # ∂H/∂t is the same as γG-F=H(1)-H(0) for our choice of homotopy - xp = x .+ Δx * step_size + # ∂H/∂t is the same as γG-F=H(1)-H(0) for our choice of homotopy + Δx = JH \ -[gg(vars=>x) for gg in H(1)-H(0)] + xh = x + Δx * step_size # Corrector step + JHh=differentiate(H(t+step_size), vars) for _ in 1:10 - JH = [jh(vars=>xp) for jh in differentiate(H(t+step_size), vars)] - Δx = JH \ -[h(vars=>xp) for h in H(t+step_size)] - xp = xp .+ Δx - if LinearAlgebra.norm(Δx) < 1e-6 + JH = [jh(vars=>xh) for jh in JHh] + Δx = JH \ -[h(vars=>xh) for h in H(t+step_size)] + xh = xh + Δx + if LinearAlgebra.norm([h(vars=>xh) for h in H(t+step_size)]) < 1e-8 break end end - return xp + return xh end end