Skip to content

Commit

Permalink
Restrict domain consistently for some functions of Taylor series with…
Browse files Browse the repository at this point in the history 
… interval coeffs (#312)

* Restrict domain consistently of some functions of Taylor series with interval coeffs

* Fix ambiguity

* Add tests and methods for pow and TaylorN

* Bump minor version
  • Loading branch information
@lbenet
lbenet authored Feb 7, 2023
1 parent 16fe479 commit 3a96af4
Show file tree
Hide file tree
Showing 4 changed files with 182 additions and 12 deletions.
2 changes: 1 addition & 1 deletion Project.toml
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
name = "TaylorSeries"
uuid = "6aa5eb33-94cf-58f4-a9d0-e4b2c4fc25ea"
repo = "https://github.com/JuliaDiff/TaylorSeries.jl.git"
version = "0.12.2"
version = "0.13.0"

[deps]
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
Expand Down
148 changes: 142 additions & 6 deletions src/intervals.jl
View file
Original file line number Diff line number Diff line change
@@ -1,15 +1,151 @@
using .IntervalArithmetic

# Method used for Taylor1{Interval{T}}^n
function ^(a::Taylor1{T}, n::Integer) where {T<:Interval}
n == 0 && return one(a)
n == 1 && return copy(a)
n == 2 && return square(a)
n < 0 && return a^float(n)
return power_by_squaring(a, n)
for T in (:Taylor1, :TaylorN)
@eval function ^(a::$T{Interval{S}}, n::Integer) where {S<:Real}
n == 0 && return one(a)
n == 1 && return copy(a)
n == 2 && return square(a)
n < 0 && return a^float(n)
return power_by_squaring(a, n)
end
end

function ^(a::Taylor1{Interval{T}}, r::S) where {T<:Real, S<:Real}
a0 = constant_term(a) ∩ Interval(zero(T), T(Inf))
aux = one(a0)^r

iszero(r) && return Taylor1(aux, a.order)
aa = one(aux) * a
aa[0] = one(aux) * a0
r == 1 && return aa
r == 2 && return square(aa)
r == 1/2 && return sqrt(aa)

l0 = findfirst(a)
lnull = trunc(Int, r*l0 )
if (a.order-lnull < 0) || (lnull > a.order)
return Taylor1( zero(aux), a.order)
end
c_order = l0 == 0 ? a.order : min(a.order, trunc(Int,r*a.order))
c = Taylor1(zero(aux), c_order)
for k = 0:c_order
pow!(c, aa, r, k)
end

return c
end
function ^(a::TaylorN{Interval{T}}, r::S) where {T<:Real, S<:Real}
a0 = constant_term(a) ∩ Interval(zero(T), T(Inf))
a0r = a0^r
aux = one(a0r)

iszero(r) && return TaylorN(aux, a.order)
aa = aux * a
aa[0] = aux * a0
r == 1 && return aa
r == 2 && return square(aa)
r == 1/2 && return sqrt(aa)
isinteger(r) && return aa^round(Int,r)

# @assert !iszero(a0)
iszero(a0) && throw(DomainError(a,
"""The 0-th order TaylorN coefficient must be non-zero
in order to expand `^` around 0."""))

c = TaylorN( a0r, a.order)
for ord in 1:a.order
pow!(c, aa, r, ord)
end

return c
end


for T in (:Taylor1, :TaylorN)
@eval function log(a::$T{Interval{S}}) where {S<:Real}
iszero(constant_term(a)) && throw(DomainError(a,
"""The 0-th order coefficient must be non-zero in order to expand `log` around 0."""))
a0 = constant_term(a) ∩ Interval(S(0.0), S(Inf))
order = a.order
aux = log(a0)
aa = one(aux) * a
aa[0] = one(aux) * a0
c = $T( aux, order )
for k in eachindex(a)
log!(c, aa, k)
end
return c
end

@eval function asin(a::$T{Interval{S}}) where {S<:Real}
a0 = constant_term(a) ∩ Interval(-one(S), one(S))
a0^2 == one(a0) && throw(DomainError(a,
"""Series expansion of asin(x) diverges at x = ±1."""))

order = a.order
aux = asin(a0)
aa = one(aux) * a
aa[0] = one(aux) * a0
c = $T( aux, order )
r = $T( sqrt(1 - a0^2), order )
for k in eachindex(a)
asin!(c, aa, r, k)
end
return c
end

@eval function acos(a::$T{Interval{S}}) where {S<:Real}
a0 = constant_term(a) ∩ Interval(-one(S), one(S))
a0^2 == one(a0) && throw(DomainError(a,
"""Series expansion of asin(x) diverges at x = ±1."""))

order = a.order
aux = acos(a0)
aa = one(aux) * a
aa[0] = one(aux) * a0
c = $T( aux, order )
r = $T( sqrt(1 - a0^2), order )
for k in eachindex(a)
acos!(c, aa, r, k)
end
return c
end

@eval function acosh(a::$T{Interval{S}}) where {S<:Real}
a0 = constant_term(a) ∩ Interval(one(S), S(Inf))
a0^2 == one(a0) && throw(DomainError(a,
"""Series expansion of acosh(x) diverges at x = ±1."""))

order = a.order
aux = acosh(a0)
aa = one(aux) * a
aa[0] = one(aux) * a0
c = $T( aux, order )
r = $T( sqrt(a0^2 - 1), order )
for k in eachindex(a)
acosh!(c, aa, r, k)
end
return c
end

@eval function atanh(a::$T{Interval{S}}) where {S<:Real}
order = a.order
a0 = constant_term(a) ∩ Interval(-one(S), one(S))
aux = atanh(a0)
aa = one(aux) * a
aa[0] = one(aux) * a0
c = $T( aux, order)
r = $T(one(aux) - a0^2, order)
iszero(constant_term(r)) && throw(DomainError(a,
"""Series expansion of atanh(x) diverges at x = ±1."""))

for k in eachindex(a)
atanh!(c, aa, r, k)
end
return c
end
end

function evaluate(a::Taylor1, dx::Interval)
order = a.order
Expand Down
6 changes: 1 addition & 5 deletions src/power.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -85,11 +85,8 @@ end

## Real power ##
function ^(a::Taylor1, r::S) where {S<:Real}
# println()
a0 = constant_term(a)
# @show(a, a0)
aux = one(a0)^r
# @show(aux)

iszero(r) && return Taylor1(aux, a.order)
aa = aux*a
Expand All @@ -103,12 +100,11 @@ function ^(a::Taylor1, r::S) where {S<:Real}
return Taylor1( zero(aux), a.order)
end
c_order = l0 == 0 ? a.order : min(a.order, trunc(Int,r*a.order))
# @show(c_order)
c = Taylor1(zero(aux), c_order)
for k = 0:c_order
pow!(c, aa, r, k)
end
# println()

return c
end

Expand Down
38 changes: 38 additions & 0 deletions test/intervals.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -99,4 +99,42 @@ eeuler = Base.MathConstants.e
displayBigO(true)
@test string(Taylor1(vc, 5)) ==
" ( [1.5, 2] + [0, 0]im ) - ( [1, 2] + [-1, 1]im ) t + ( [-1, 1.5] + [-1, 1.5]im ) t² + ( [0, 0] + [-1, 1.5]im ) t³ + 𝒪(t⁶)"

# Tests related to Iss #311
# `sqrt` and `pow` defined on Interval(0,Inf)
@test_throws DomainError sqrt(ti)
@test sqrt(Interval(0.0, 1.e-15) + ti) == sqrt(Interval(-1.e-15, 1.e-15) + ti)
aa = sqrt(sqrt(Interval(0.0, 1.e-15) + ti))
@test aa == sqrt(sqrt(Interval(-1.e-15, 1.e-15) + ti))
bb = (Interval(0.0, 1.e-15) + ti)^(1/4)
@test bb == (Interval(-1.e-15, 1.e-15) + ti)^(1/4)
@test all(aa.coeffs[2:end] .⊂ bb.coeffs[2:end])
@test_throws DomainError sqrt(x)
@test sqrt(Interval(-1,1)+x) == sqrt(Interval(0,1)+x)
@test (Interval(-1,1)+x)^(1/4) == (Interval(0,1)+x)^(1/4)

# `log` defined on Interval(0,Inf)
@test_throws DomainError log(ti)
@test log(Interval(0.0, 1.e-15) + ti) == log(Interval(-1.e-15, 1.e-15) + ti)
@test_throws DomainError log(y)
@test log(Interval(0.0, 1.e-15) + y) == log(Interval(-1.e-15, 1.e-15) + y)
# `asin` and `acos` defined on Interval(-1,1)
@test_throws DomainError asin(Interval(1.0 .. 2.0) + ti)
@test asin(Interval(-2.0 .. 0.0) + ti) == asin(Interval(-1,0) + ti)
@test_throws DomainError acos(Interval(1.0 .. 2.0) + ti)
@test acos(Interval(-2.0 .. 0.0) + ti) == acos(Interval(-1,0) + ti)
@test_throws DomainError asin(Interval(1.0 .. 2.0) + x)
@test asin(Interval(-2.0 .. 0.0) + x) == asin(Interval(-1,0) + x)
@test_throws DomainError acos(Interval(1.0 .. 2.0) + x)
@test acos(Interval(-2.0 .. 0.0) + x) == acos(Interval(-1,0) + x)
# acosh defined on Interval(1,Inf)
@test_throws DomainError acosh(Interval(0.0 .. 1.0) + ti)
@test acosh(Interval(0.0 .. 2.0) + ti) == acosh(Interval(1.0 .. 2.0) + ti)
@test_throws DomainError acosh(Interval(0.0 .. 1.0) + x)
@test acosh(Interval(0.0 .. 2.0) + x) == acosh(Interval(1.0 .. 2.0) + x)
# atanh defined on Interval(-1,1)
@test_throws DomainError atanh(Interval(1.0 .. 1.0) + ti)
@test atanh(Interval(-2.0 .. 0.0) + ti) == atanh(Interval(-1.0 .. 0.0) + ti)
@test_throws DomainError atanh(Interval(1.0 .. 1.0) + y)
@test atanh(Interval(-2.0 .. 0.0) + y) == atanh(Interval(-1.0 .. 0.0) + y)
end

2 comments on commit 3a96af4

@lbenet
Copy link
Member Author

@lbenet lbenet commented on 3a96af4 Feb 7, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@JuliaRegistrator register

@JuliaRegistrator
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Registration pull request created: JuliaRegistries/General/77232

After the above pull request is merged, it is recommended that a tag is created on this repository for the registered package version.

This will be done automatically if the Julia TagBot GitHub Action is installed, or can be done manually through the github interface, or via:

git tag -a v0.13.0 -m "<description of version>" 3a96af40d73f381d6079aeea634336f328272744
git push origin v0.13.0

Please sign in to comment.