WIP: Integration for TaylorModelN (#82)

* add integration for TaylorModelN

* add bound_integration for TMN and a new method for integrate

* minor change to integrate for TMN

* add tests for integrate

* add tests for TMN integration

* add more tests for TMN integration

* add more tests for TMN integration

* improve tests readability

* modify check_containment to re-use get_random_point function

* remove x0 parameter from integrate function

* add docstring to integrate for TMN

src/integration.jl

@@ -43,6 +43,25 @@ function integrate(a::TaylorModel1{TaylorModelN{N,T,S},S},
     return TaylorModel1( integ_pol, ΔN, a.x0, a.dom )
 end
 
+"""
+    integrate(fT, which)
+
+Integrates a `TaylorModelN` with respect to `which` variable.
+The returned `TaylorModelN` corresponds to the Taylor Model
+of the definite integral ∫f(x) - ∫f(expansion_point).
+"""
+function integrate(fT::TaylorModelN, which=1)
+    p̂ = integrate(fT.pol, which)
+    order = fT.pol.order
+    r = TaylorN(p̂.coeffs[1:order+1])
+    s = TaylorN(p̂.coeffs[order+2:end])
+    Δ = bound_integration(fT, s, which)
+    return TaylorModelN(r, Δ, fT.x0, fT.dom)
+end
+function integrate(fT::TaylorModelN, s::Symbol)
+    which = TaylorSeries.lookupvar(s)
+    return integrate(fT, which)
+end
 
 """
     bound_integration(xTM1::TaylorModel1{Interval{S},S}, δt::Interval{S})
@@ -66,8 +85,10 @@ function bound_integration(a::Vector{TaylorModel1{T,S}}, δ) where {T,S}
     Δ = δ .* (remainder.(a) .+ getcoeff.(polynomial.(a), order) .* aux)
     return IntervalBox(Δ)
 end
-
-
+function bound_integration(fT::TaylorModelN, s::TaylorN, which)
+    Δ = s(fT.dom - fT.x0) + fT.rem * (fT.dom[which] - fT.x0[which])
+    return Δ
+end
 
 function picard_lindelof(f!, dxTM1TMN::Vector{TaylorModel1{T,S}},
         xTM1TMN::Vector{TaylorModel1{T,S}}, t, params) where {T,S}

test/TMN.jl

@@ -8,10 +8,14 @@ const _num_tests = 1000
 
 setformat(:full)
 
+function get_random_point(ib0::IntervalBox{N, T}) where {N, T}
+    xmid = mid(ib0)
+    return diam.(ib0) .* (rand(N) .- 0.5) .+ xmid
+end
 
 function check_containment(ftest, xx::TaylorModelN{N,T,S}, tma::TaylorModelN{N,T,S}) where {N,T,S}
     x0 = expansion_point(tma)
-    xfp = diam.(tma.dom) .* (rand(N) .- 0.5) .+ mid(x0)
+    xfp = get_random_point(tma.dom)
     xbf = [big(xfp[i]) for i=1:N]
     ib = IntervalBox([@interval(xfp[i]) for i=1:N]...)
     range = evaluate(tma, ib-x0)
@@ -243,6 +247,105 @@ set_variables(Interval{Float64}, [:x, :y], order=_order_max)
         @test sup(norm(tmb.pol - (1+ym).pol, Inf)) < 1.0e-15
     end
 
+    @testset "Tests for integrate" begin
+        ib0 = (0. .. 1.) × (0. .. 1.)
+        b0 = (0.5 .. 0.5) × (0.5 .. 0.5)
+        xm = TaylorModelN(1, _order, b0, ib0)
+        ym = TaylorModelN(2, _order, b0, ib0)
+
+        f(x, y) = cos(x)
+        ∫fdx(x, y) = sin(x)
+        ∫fdy(x, y) = cos(x) * y
+        fT = f(xm, ym)
+        ∫fTdx = integrate(fT, :x)
+        ∫fTdy = integrate(fT, :y)
+
+        for ind in 1:_num_tests
+            xtest = get_random_point(ib0)
+            cx = [mid(ib0[1]), xtest[2]]
+            cy = [xtest[1], mid(ib0[2])]
+            aux = IntervalBox(xtest) - b0
+            @test (∫fdx(xtest...) - ∫fdx(cx...)) ∈ ∫fTdx(aux)
+            @test (∫fdy(xtest...) - ∫fdy(cy...)) ∈ ∫fTdy(aux)
+        end
+
+        f(x, y) = sin(x) * cos(y)
+        ∫fdx(x, y) = -cos(x) * cos(y)
+        ∫fdy(x, y) = sin(x) * sin(y)
+        fT = f(xm, ym)
+        ∫fTdx = integrate(fT, :x)
+        ∫fTdy = integrate(fT, :y)
+
+        for ind in 1:_num_tests
+            xtest = get_random_point(ib0)
+            cx = [mid(ib0[1]), xtest[2]]
+            cy = [xtest[1], mid(ib0[2])]
+            aux = IntervalBox(xtest) - b0
+            @test (∫fdx(xtest...) - ∫fdx(cx...)) ∈ ∫fTdx(aux)
+            @test (∫fdy(xtest...) - ∫fdy(cy...)) ∈ ∫fTdy(aux)
+        end
+
+        f(x, y) = exp(x)
+        ∫fdx(x, y) = exp(x)
+        ∫fdy(x, y) = exp(x) * y
+        fT = f(xm, ym)
+        ∫fTdx = integrate(fT, :x)
+        ∫fTdy = integrate(fT, :y)
+
+        for ind in 1:_num_tests
+            xtest = get_random_point(ib0)
+            cx = [mid(ib0[1]), xtest[2]]
+            cy = [xtest[1], mid(ib0[2])]
+            aux = IntervalBox(xtest) - b0
+            @test (∫fdx(xtest...) - ∫fdx(cx...)) ∈ ∫fTdx(aux)
+            @test (∫fdy(xtest...) - ∫fdy(cy...)) ∈ ∫fTdy(aux)
+        end
+
+        f(x, y) = log(x) * x^2 + cos(x * y) + sin(x * y)
+        ∫fdx(x, y) = (x^3 * y * (3log(x) - 1) + 9sin(x * y) - 9cos(x * y)) / 9y
+        ∫fdy(x, y) = (x^3 * y * log(x) + sin(x * y) - cos(x * y)) / x
+        fT = f(xm, ym)
+        ∫fTdx = integrate(fT, :x)
+        ∫fTdy = integrate(fT, :y)
+
+        for ind in 1:_num_tests
+            xtest = get_random_point(ib0)
+            cx = [mid(ib0[1]), xtest[2]]
+            cy = [xtest[1], mid(ib0[2])]
+            aux = IntervalBox(xtest) - b0
+            @test (∫fdx(xtest...) - ∫fdx(cx...)) ∈ ∫fTdx(aux)
+            @test (∫fdy(xtest...) - ∫fdy(cy...)) ∈ ∫fTdy(aux)
+        end
+
+        f(x, y) = exp(-0.5 * (x^2 + y^2)) * x
+        ∫fdx(x, y) = -exp(-0.5 * (x^2 + y^2))
+        fT = f(xm, ym)
+        ∫fTdx = integrate(fT, :x)
+
+        for ind in 1:_num_tests
+            xtest = get_random_point(ib0)
+            cx = [mid(ib0[1]), xtest[2]]
+            aux = IntervalBox(xtest) - b0
+            @test (∫fdx(xtest...) - ∫fdx(cx...)) ∈ ∫fTdx(aux)
+        end
+
+        f(x, y) = x * cos(y) * exp(x + y)
+        ∫fdx(x, y) = (x - 1) * exp(x + y) * cos(y)
+        ∫fdy(x, y) = 0.5 * x * exp(x + y) * (sin(y) + cos(y))
+        fT = f(xm, ym)
+        ∫fTdx = integrate(fT, :x)
+        ∫fTdy = integrate(fT, :y)
+
+        for ind in 1:_num_tests
+            xtest = get_random_point(ib0)
+            cx = [mid(ib0[1]), xtest[2]]
+            cy = [xtest[1], mid(ib0[2])]
+            aux = IntervalBox(xtest) - b0
+            @test (∫fdx(xtest...) - ∫fdx(cx...)) ∈ ∫fTdx(aux)
+            @test (∫fdy(xtest...) - ∫fdy(cy...)) ∈ ∫fTdy(aux)
+        end
+    end
+
     @testset "Display" begin
         xm = TaylorModelN(1, _order, b1, ib1)
         ym = TaylorModelN(2, _order, b1, ib1)