Skip to content

Commit

Permalink
Merge pull request #56 from huiyuxie/main
Browse files Browse the repository at this point in the history
Fix Issue #34: Add new method Alefeld, Potra, Shi (1995)
  • Loading branch information
ChrisRackauckas authored Mar 30, 2023
2 parents 622be85 + c50cd34 commit 58ffb58
Show file tree
Hide file tree
Showing 4 changed files with 200 additions and 3 deletions.
2 changes: 1 addition & 1 deletion lib/SimpleNonlinearSolve/Project.toml
Original file line number Diff line number Diff line change
Expand Up @@ -43,4 +43,4 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[targets]
test = ["BenchmarkTools", "SafeTestsets", "Pkg", "Test", "StaticArrays", "NNlib"]
test = ["BenchmarkTools", "SafeTestsets", "Pkg", "Test", "StaticArrays", "NNlib"]
5 changes: 3 additions & 2 deletions lib/SimpleNonlinearSolve/src/SimpleNonlinearSolve.jl
Original file line number Diff line number Diff line change
Expand Up @@ -38,6 +38,7 @@ include("brent.jl")
include("dfsane.jl")
include("ad.jl")
include("halley.jl")
include("alefeld.jl")

import SnoopPrecompile

Expand All @@ -59,13 +60,13 @@ SnoopPrecompile.@precompile_all_calls begin for T in (Float32, Float64)
=#

prob_brack = IntervalNonlinearProblem{false}((u, p) -> u * u - p, T.((0.0, 2.0)), T(2))
for alg in (Bisection, Falsi, Ridder, Brent)
for alg in (Bisection, Falsi, Ridder, Brent, Alefeld)
solve(prob_brack, alg(), abstol = T(1e-2))
end
end end

# DiffEq styled algorithms
export Bisection, Brent, Broyden, LBroyden, SimpleDFSane, Falsi, Halley, Klement,
Ridder, SimpleNewtonRaphson, SimpleTrustRegion
Ridder, SimpleNewtonRaphson, SimpleTrustRegion, Alefeld

end # module
184 changes: 184 additions & 0 deletions lib/SimpleNonlinearSolve/src/alefeld.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,184 @@
"""
`Alefeld()`
An implementation of algorithm 4.2 from [Alefeld](https://dl.acm.org/doi/10.1145/210089.210111).
The paper brought up two new algorithms. Here choose to implement algorithm 4.2 rather than
algorithm 4.1 because, in certain sense, the second algorithm(4.2) is an optimal procedure.
"""
struct Alefeld <: AbstractBracketingAlgorithm end

function SciMLBase.solve(prob::IntervalNonlinearProblem,
alg::Alefeld, args...; abstol = nothing,
reltol = nothing,
maxiters = 1000, kwargs...)

f = Base.Fix2(prob.f, prob.p)
a, b = prob.tspan
c = a - (b - a) / (f(b) - f(a)) * f(a)

fc = f(c)
if iszero(fc)
return SciMLBase.build_solution(prob, alg, c, fc;
retcode = ReturnCode.Success,
left = a,
right = b)
end
a, b, d = _bracket(f, a, b, c)
e = zero(a) # Set e as 0 before iteration to avoid a non-value f(e)

# Begin of algorithm iteration
for i in 2:maxiters
# The first bracketing block
f₁, f₂, f₃, f₄ = f(a), f(b), f(d), f(e)
if i == 2 || (f₁ == f₂ || f₁ == f₃ || f₁ == f₄ || f₂ == f₃ || f₂ == f₄ || f₃ == f₄)
c = _newton_quadratic(f, a, b, d, 2)
else
c = _ipzero(f, a, b, d, e)
if (c - a) * (c - b) 0
c = _newton_quadratic(f, a, b, d, 2)
end
end
ē, fc = d, f(c)
(a == c || b == c) &&
return SciMLBase.build_solution(prob, alg, c, fc;
retcode = ReturnCode.FloatingPointLimit,
left = a,
right = b)
iszero(fc) &&
return SciMLBase.build_solution(prob, alg, c, fc;
retcode = ReturnCode.Success,
left = a,
right = b)
ā, b̄, d̄ = _bracket(f, a, b, c)

# The second bracketing block
f₁, f₂, f₃, f₄ = f(ā), f(b̄), f(d̄), f(ē)
if f₁ == f₂ || f₁ == f₃ || f₁ == f₄ || f₂ == f₃ || f₂ == f₄ || f₃ == f₄
c = _newton_quadratic(f, ā, b̄, d̄, 3)
else
c = _ipzero(f, ā, b̄, d̄, ē)
if (c - ā) * (c - b̄) 0
c = _newton_quadratic(f, ā, b̄, d̄, 3)
end
end
fc = f(c)
(ā == c ||== c) &&
return SciMLBase.build_solution(prob, alg, c, fc;
retcode = ReturnCode.FloatingPointLimit,
left = ā,
right = b̄)
iszero(fc) &&
return SciMLBase.build_solution(prob, alg, c, fc;
retcode = ReturnCode.Success,
left = ā,
right = b̄)
ā, b̄, d̄ = _bracket(f, ā, b̄, c)

# The third bracketing block
if abs(f(ā)) < abs(f(b̄))
u =
else
u =
end
c = u - 2 * (b̄ - ā) / (f(b̄) - f(ā)) * f(u)
if (abs(c - u)) > 0.5 * (b̄ - ā)
c = 0.5 * (ā + b̄)
end
fc = f(c)
(ā == c ||== c) &&
return SciMLBase.build_solution(prob, alg, c, fc;
retcode = ReturnCode.FloatingPointLimit,
left = ā,
right = b̄)
iszero(fc) &&
return SciMLBase.build_solution(prob, alg, c, fc;
retcode = ReturnCode.Success,
left = ā,
right = b̄)
ā, b̄, d = _bracket(f, ā, b̄, c)

# The last bracketing block
if-< 0.5 * (b - a)
a, b, e = ā, b̄, d̄
else
e = d
c = 0.5 * (ā + b̄)
fc = f(c)
(ā == c ||== c) &&
return SciMLBase.build_solution(prob, alg, c, fc;
retcode = ReturnCode.FloatingPointLimit,
left = ā,
right = b̄)
iszero(fc) &&
return SciMLBase.build_solution(prob, alg, c, fc;
retcode = ReturnCode.Success,
left = ā,
right = b̄)
a, b, d = _bracket(f, ā, b̄, c)
end
end

# Reassign the value a, b, and c
if b == c
b = d
elseif a == c
a = d
end
fc = f(c)

# Reuturn solution when run out of max interation
return SciMLBase.build_solution(prob, alg, c, fc; retcode = ReturnCode.MaxIters,
left = a, right = b)
end

# Define subrotine function bracket, check fc before bracket to return solution
function _bracket(f::F, a, b, c) where F
if iszero(f(c))
ā, b̄, d = a, b, c
else
if f(a) * f(c) < 0
ā, b̄, d = a, c, b
elseif f(b) * f(c) < 0
ā, b̄, d = c, b, a
end
end

return ā, b̄, d
end

# Define subrotine function newton quadratic, return the approximation of zero
function _newton_quadratic(f::F, a, b, d, k) where F
A = ((f(d) - f(b)) / (d - b) - (f(b) - f(a)) / (b - a)) / (d - a)
B = (f(b) - f(a)) / (b - a)

if iszero(A)
return a - (1 / B) * f(a)
elseif A * f(a) > 0
rᵢ₋₁ = a
else
rᵢ₋₁ = b
end

for i in 1:k
rᵢ = rᵢ₋₁ - (f(a) + B * (rᵢ₋₁ - a) + A * (rᵢ₋₁ - a) * (rᵢ₋₁ - b)) / (B + A * (2 * rᵢ₋₁ - a - b))
rᵢ₋₁ = rᵢ
end

return rᵢ₋₁
end

# Define subrotine function ipzero, also return the approximation of zero
function _ipzero(f::F, a, b, c, d) where F
Q₁₁ = (c - d) * f(c) / (f(d) - f(c))
Q₂₁ = (b - c) * f(b) / (f(c) - f(b))
Q₃₁ = (a - b) * f(a) / (f(b) - f(a))
D₂₁ = (b - c) * f(c) / (f(c) - f(b))
D₃₁ = (a - b) * f(b) / (f(b) - f(a))
Q₂₂ = (D₂₁ - Q₁₁) * f(b) / (f(d) - f(b))
Q₃₂ = (D₃₁ - Q₂₁) * f(a) / (f(c) - f(a))
D₃₂ = (D₃₁ - Q₂₁) * f(c) / (f(c) - f(a))
Q₃₃ = (D₃₂ - Q₂₂) * f(a) / (f(d) - f(a))

return a + Q₃₁ + Q₃₂ + Q₃₃
end
12 changes: 12 additions & 0 deletions lib/SimpleNonlinearSolve/test/basictests.jl
Original file line number Diff line number Diff line change
Expand Up @@ -210,6 +210,18 @@ for p in 1.1:0.1:100.0
@test ForwardDiff.derivative(g, p) 1 / (2 * sqrt(p))
end

# Alefeld
g = function (p)
probN = IntervalNonlinearProblem{false}(f, typeof(p).(tspan), p)
sol = solve(probN, Alefeld())
return sol.u
end

for p in 1.1:0.1:100.0
@test g(p) sqrt(p)
@test ForwardDiff.derivative(g, p) 1 / (2 * sqrt(p))
end

f, tspan = (u, p) -> p[1] * u * u - p[2], (1.0, 100.0)
t = (p) -> [sqrt(p[2] / p[1])]
p = [0.9, 50.0]
Expand Down

0 comments on commit 58ffb58

Please sign in to comment.