Skip to content

Commit

Permalink
correction testNF
Browse files Browse the repository at this point in the history
  • Loading branch information
rveltz committed Nov 22, 2024
1 parent f48bb60 commit 997cae3
Show file tree
Hide file tree
Showing 2 changed files with 16 additions and 16 deletions.
10 changes: 5 additions & 5 deletions src/continuation/Palc.jl
Original file line number Diff line number Diff line change
Expand Up @@ -388,7 +388,7 @@ function newton_palc(iter::AbstractContinuationIterable,
ϵ = getdelta(prob)
paramlens = getlens(iter)
contparams = getcontparams(iter)
T = eltype(iter)
𝒯 = eltype(iter)
θ = getθ(iter)

z0 = getsolution(state)
Expand All @@ -415,13 +415,13 @@ function newton_palc(iter::AbstractContinuationIterable,
res_f = residual(prob, x, set(par, paramlens, p)); res_n = N(x, p)

dX = _copy(res_f)
dp = zero(T)
up = zero(T)
dp = zero(𝒯)
up = zero(𝒯)

# dFdp = (F(x, p + ϵ) - res_f) / ϵ
dFdp = _copy(residual(prob, x, set(par, paramlens, p + ϵ)))
minus!(dFdp, res_f) # dFdp = dFdp - res_f
rmul!(dFdp, one(T) / ϵ)
rmul!(dFdp, one(𝒯) / ϵ)

res = normAC(res_f, res_n)
residuals = [res]
Expand All @@ -436,7 +436,7 @@ function newton_palc(iter::AbstractContinuationIterable,
while (step < max_iterations) && (res > tol) && line_step && compute
# dFdp = (F(x, p + ϵ) - F(x, p)) / ϵ)
copyto!(dFdp, residual(prob, x, set(par, paramlens, p + ϵ)))
minus!(dFdp, res_f); rmul!(dFdp, one(T) / ϵ)
minus!(dFdp, res_f); rmul!(dFdp, one(𝒯) / ϵ)

# compute jacobian
J = jacobian(prob, x, set(par, paramlens, p))
Expand Down
22 changes: 11 additions & 11 deletions test/testNF.jl
Original file line number Diff line number Diff line change
Expand Up @@ -156,19 +156,19 @@ let
@test norm(bp2d.nf.b2, Inf) < 3e-6
@test norm(bp2d.nf.b1 - prob2d.params.α * 3.23 * I, Inf) < 1e-9
@test norm(bp2d.nf.a, Inf) < 1e-6

# same but when the eigenvalues are not saved in the branch but computed on the fly instead
br_noev = BK.continuation( prob2d, PALC(), ContinuationPar(opts_br; n_inversion = 2, save_eigenvectors = false); normC = norminf)
bp2d = BK.get_normal_form(br_noev, 1; ζs = [[1, 0, 0.], [0, 1, 0.]]);
@test abs(bp2d.nf.b3[1,1,1,1] / 6 - -prob2d.params.α * 0.123) < 1e-15
@test abs(bp2d.nf.b3[1,1,2,2] / 2 - -prob2d.params.α * 0.234) < 1e-15
@test abs(bp2d.nf.b3[1,1,1,2] / 2 - -prob2d.params.α * 0.0) < 1e-15
@test abs(bp2d.nf.b3[2,1,1,2] / 2 - -prob2d.params.α * 0.456) < 1e-15
@test norm(bp2d.nf.b2, Inf) < 3e-15
@test norm(bp2d.nf.b1 - prob2d.params.α * 3.23 * I, Inf) < 1e-9
@test norm(bp2d.nf.a, Inf) < 1e-15
end
end
##############################
# same but when the eigenvalues are not saved in the branch but computed on the fly instead
br_noev = BK.continuation( prob2d, PALC(), ContinuationPar(opts_br; n_inversion = 2, save_eigenvectors = false); normC = norminf)
bp2d = BK.get_normal_form(br_noev, 1; ζs = [[1, 0, 0.], [0, 1, 0.]]);
@test abs(bp2d.nf.b3[1,1,1,1] / 6 - -prob2d.params.α * 0.123) < 1e-15
@test abs(bp2d.nf.b3[1,1,2,2] / 2 - -prob2d.params.α * 0.234) < 1e-15
@test abs(bp2d.nf.b3[1,1,1,2] / 2 - -prob2d.params.α * 0.0) < 1e-15
@test abs(bp2d.nf.b3[2,1,1,2] / 2 - -prob2d.params.α * 0.456) < 1e-15
@test norm(bp2d.nf.b2, Inf) < 3e-15
@test norm(bp2d.nf.b1 - prob2d.params.α * 3.23 * I, Inf) < 1e-9
@test norm(bp2d.nf.a, Inf) < 1e-15
####################################################################################################
# vector field to test nearby secondary bifurcations
FbpSecBif(u, p) = @. -u * (p + u * (2-5u)) * (p -.15 - u * (2+20u))
Expand Down

0 comments on commit 997cae3

Please sign in to comment.