OP tests pass

Author: dlfivefifty

src/Domains/PiecewiseSegment.jl

@@ -1,3 +1,5 @@
+export PiecewiseSegment
+
 struct PiecewiseSegment{T} <: Domain{T}
     points::Vector{T}
     PiecewiseSegment{T}(d::Vector{T}) where {T} = new{T}(d)

src/Operators/banded/ToeplitzOperator.jl

@@ -161,17 +161,6 @@ function Base.inv(T::ToeplitzOperator)
     ToeplitzOperator(ai[2:end],ai[1:1])
 end
 
-function Fun(T::ToeplitzOperator)
-   if length(T.nonnegative)==1
-      Fun(Taylor(),[T.nonnegative;T.negative])
-    elseif length(T.negative)==0
-        Fun(Hardy{false}(),T.nonnegative)
-    else
-        Fun(Laurent(Circle()),interlace(T.nonnegative,T.negative))
-    end
-end
-
-
 
 ####
 # Toeplitz/Hankel

src/Spaces/HeavisideSpace.jl

@@ -18,12 +18,10 @@ convert(::Type{HeavisideSpace},d::AbstractVector) =
     HeavisideSpace(PiecewiseSegment(sort(d)))
 
 spacescompatible(a::SplineSpace{λ},b::SplineSpace{λ}) where {λ} = domainscompatible(a,b)
-canonicalspace(sp::HeavisideSpace) = PiecewiseSpace(map(Chebyshev,components(domain(sp))))
 
 
-function evaluate(f::Fun{HeavisideSpace{T,R}},x::Real) where {T<:Real,R}
-    p = domain(f).points
-    c = f.coefficients
+function evaluate(c::AbstractVector{T}, s::HeavisideSpace{<:Real}, x::Real) where T
+    p = domain(s).points
     for k=1:length(c)
         if p[k] ≤ x ≤ p[k+1]
             return c[k]
@@ -33,7 +31,7 @@ function evaluate(f::Fun{HeavisideSpace{T,R}},x::Real) where {T<:Real,R}
 end
 
 
-function evaluate(f::Fun{SplineSpace{1,T,R}},x::Real) where {T<:Real,R}
+function evaluate(c::AbstractVector{T}, s::SplineSpace{1,<:Real}, x::Real) where T
     p = domain(f).points
     c = f.coefficients
     for k=1:length(p)-1

src/Spaces/ProductSpaceOperators.jl

@@ -1,4 +1,4 @@
-
+export continuity
 
 ## Space promotion for InterlaceOperator
 # It's here because we need DirectSumSpace

src/specialfunctions.jl

@@ -551,4 +551,42 @@ function roots(f::Fun{P}) where P<:PiecewiseSpace
     rts
 end
 
-roots(f::Fun{S,T}) where {S<:PointSpace,T} = space(f).points[values(f) .== 0]
+roots(f::Fun{S,T}) where {S<:PointSpace,T} = space(f).points[values(f) .== 0]
+
+
+
+#
+# These formulæ, appearing in Eq. (2.5) of:
+#
+# A.-K. Kassam and L. N. Trefethen, Fourth-order time-stepping for stiff PDEs, SIAM J. Sci. Comput., 26:1214--1233, 2005,
+#
+# are derived to implement ETDRK4 in double precision without numerical instability from cancellation.
+#
+
+expα_asy(x) = (exp(x)*(4-3x+x^2)-4-x)/x^3
+expβ_asy(x) = (exp(x)*(x-2)+x+2)/x^3
+expγ_asy(x) = (exp(x)*(4-x)-4-3x-x^2)/x^3
+
+# TODO: General types
+
+expα_taylor(x::Union{Float64,ComplexF64}) = @evalpoly(x,1/6,1/6,3/40,1/45,5/1008,1/1120,7/51840,1/56700,1/492800,1/4790016,11/566092800,1/605404800,13/100590336000,1/106748928000,1/1580833013760,1/25009272288000,17/7155594141696000,1/7508956815360000)
+expβ_taylor(x::Union{Float64,ComplexF64}) = @evalpoly(x,1/6,1/12,1/40,1/180,1/1008,1/6720,1/51840,1/453600,1/4435200,1/47900160,1/566092800,1/7264857600,1/100590336000,1/1494484992000,1/23712495206400,1/400148356608000,1/7155594141696000,1/135161222676480000)
+expγ_taylor(x::Union{Float64,ComplexF64}) = @evalpoly(x,1/6,0/1,-1/120,-1/360,-1/1680,-1/10080,-1/72576,-1/604800,-1/5702400,-1/59875200,-1/691891200,-1/8717829120,-1/118879488000,-1/1743565824000,-1/27360571392000,-1/457312407552000,-1/8109673360588800)
+
+expα(x::Float64) = abs(x) < 17/16 ? expα_taylor(x) : expα_asy(x)
+expβ(x::Float64) = abs(x) < 19/16 ? expβ_taylor(x) : expβ_asy(x)
+expγ(x::Float64) = abs(x) < 15/16 ? expγ_taylor(x) : expγ_asy(x)
+
+expα(x::ComplexF64) = abs2(x) < (17/16)^2 ? expα_taylor(x) : expα_asy(x)
+expβ(x::ComplexF64) = abs2(x) < (19/16)^2 ? expβ_taylor(x) : expβ_asy(x)
+expγ(x::ComplexF64) = abs2(x) < (15/16)^2 ? expγ_taylor(x) : expγ_asy(x)
+
+expα(x) = expα_asy(x)
+expβ(x) = expβ_asy(x)
+expγ(x) = expγ_asy(x)
+
+
+for f in (:(exp),:(expm1),:expα,:expβ,:expγ)
+    @eval $f(op::Operator) = OperatorFunction(op,$f)
+end
+

test/ETDRK4Test.jl

@@ -0,0 +1,16 @@
+using ApproxFunBase, Test
+import ApproxFunBase: expα, expβ, expγ
+
+#
+# This tests the numerical stability of the formulæ appearing in Eq. (2.5) of:
+#
+# A.-K. Kassam and L. N. Trefethen, Fourth-order time-stepping for stiff PDEs, SIAM J. Sci. Comput., 26:1214--1233, 2005.
+#
+
+@testset "ETDRK4" begin
+    x = 10 .^ range(-15, stop=2, length=18)
+
+    @test norm(expα.(x)./expα.(big.(x)).-1,Inf) < eps()
+    @test norm(expβ.(x)./expβ.(big.(x)).-1,Inf) < eps()
+    @test norm(expγ.(x)./expγ.(big.(x)).-1,Inf) < eps()
+end

test/runtests.jl

@@ -47,3 +47,22 @@ end
 
 @time include("MatrixTest.jl")
 @time include("SpacesTest.jl")
+
+
+@testset "blockbandwidths for FiniteOperator of pointscompatibleace bug" begin
+    S = ApproxFunBase.PointSpace([1.0,2.0])
+    @test ApproxFunBase.blockbandwidths(FiniteOperator([1 2; 3 4],S,S)) == (0,0)
+end
+
+@testset "DiracDelta sampling" begin
+    δ = 0.3DiracDelta(0.1) + 3DiracDelta(2.3)
+    Random.seed!(0)
+    for _=1:10
+        @test sample(δ) ∈ [0.1, 2.3]
+    end
+    Random.seed!(0)
+    r = sample(δ, 10_000)
+    @test count(i -> i == 0.1, r)/length(r) ≈ 0.3/(3.3) atol=0.01
+end
+
+@time include("ETDRK4Test.jl")