Skip to content
This repository was archived by the owner on May 15, 2025. It is now read-only.

format #59

Merged
merged 1 commit into from
Mar 31, 2023
Merged

format #59

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
85 changes: 43 additions & 42 deletions src/alefeld.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,18 +9,17 @@ algorithm 4.1 because, in certain sense, the second algorithm(4.2) is an optimal
struct Alefeld <: AbstractBracketingAlgorithm end

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

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,
retcode = ReturnCode.Success,
left = a,
right = b)
end
Expand All @@ -33,47 +32,47 @@ function SciMLBase.solve(prob::IntervalNonlinearProblem,
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
else
c = _ipzero(f, a, b, d, e)
if (c - a) * (c - b) ≥ 0
c = _newton_quadratic(f, a, b, d, 2)
end
end
end
ē, fc = d, f(c)
(a == c || b == c) &&
(a == c || b == c) &&
return SciMLBase.build_solution(prob, alg, c, fc;
retcode = ReturnCode.FloatingPointLimit,
left = a,
right = b)
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)
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
else
c = _ipzero(f, ā, b̄, d̄, ē)
if (c - ā) * (c - b̄) ≥ 0
c = _newton_quadratic(f, ā, b̄, d̄, 3)
end
end
fc = f(c)
(ā == c || b̄ == c) &&
(ā == c || b̄ == c) &&
return SciMLBase.build_solution(prob, alg, c, fc;
retcode = ReturnCode.FloatingPointLimit,
left = ā,
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)
retcode = ReturnCode.Success,
left = ā,
right = b̄)
ā, b̄, d̄ = _bracket(f, ā, b̄, c)

# The third bracketing block
if abs(f(ā)) < abs(f(b̄))
Expand All @@ -86,17 +85,17 @@ function SciMLBase.solve(prob::IntervalNonlinearProblem,
c = 0.5 * (ā + b̄)
end
fc = f(c)
(ā == c || b̄ == c) &&
(ā == c || b̄ == c) &&
return SciMLBase.build_solution(prob, alg, c, fc;
retcode = ReturnCode.FloatingPointLimit,
left = ā,
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)
retcode = ReturnCode.Success,
left = ā,
right = b̄)
ā, b̄, d = _bracket(f, ā, b̄, c)

# The last bracketing block
if b̄ - ā < 0.5 * (b - a)
Expand All @@ -105,14 +104,14 @@ function SciMLBase.solve(prob::IntervalNonlinearProblem,
e = d
c = 0.5 * (ā + b̄)
fc = f(c)
(ā == c || b̄ == c) &&
(ā == c || b̄ == c) &&
return SciMLBase.build_solution(prob, alg, c, fc;
retcode = ReturnCode.FloatingPointLimit,
left = ā,
left = ā,
right = b̄)
iszero(fc) &&
return SciMLBase.build_solution(prob, alg, c, fc;
retcode = ReturnCode.Success,
retcode = ReturnCode.Success,
left = ā,
right = b̄)
a, b, d = _bracket(f, ā, b̄, c)
Expand All @@ -133,43 +132,45 @@ function SciMLBase.solve(prob::IntervalNonlinearProblem,
end

# Define subrotine function bracket, check fc before bracket to return solution
function _bracket(f::F, a, b, c) where F
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
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
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)
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ᵢ₋₁ = a
else
rᵢ₋₁ = b
end
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ᵢ₋₁ -
(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
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))
Expand All @@ -181,4 +182,4 @@ function _ipzero(f::F, a, b, c, d) where F
Q₃₃ = (D₃₂ - Q₂₂) * f(a) / (f(d) - f(a))

return a + Q₃₁ + Q₃₂ + Q₃₃
end
end