Integration of TaylorN series (#136)

* Implement integrate for TaylorN variables with tests

* Add docs for integrate, fix one case, and a test

docs/src/userguide.md

@@ -317,7 +317,7 @@ exy[8][6] # get the 6th coeff of the 8th order term
 
 Partial differentiation is also implemented for [`TaylorN`](@ref) objects,
 through the function [`derivative`](@ref), specifying the number
-of the variable as the second argument; integration is yet to be implemented.
+of the variable as the second argument.
 
 ```@repl userguide
 p = x^3 + 2x^2 * y - 7x + 2
@@ -333,6 +333,19 @@ an error is thrown.
 derivative( q, 3 )   # error, since we are dealing with 2 variables
 ```
 
+Integration with respect to the r-th variable for
+`HomogeneousPolynomial`s and `TaylorN` objects is obtained
+using [`integrate`](@ref). Note that `integrate` for `TaylorN`
+objects allows to specify a constant of integration, which must
+be independent from the integrated variable.
+
+```@repl userguide
+integrate( derivative( p, 1 ), 1) # integrate with respect to the first variable
+integrate( derivative( p, 1 ), 1, 2) # integration constant is 2
+integrate( derivative( q, 2 ), 2, -x^4) == q
+integrate( derivative( q, 2 ), 2, y)
+```
+
 [`evaluate`](@ref) can also be used for [`TaylorN`](@ref) objects, using
 it on vectors of
 numbers (`Real` or `Complex`); the length of the vector must coincide with the

src/calculus.jl

@@ -15,7 +15,6 @@ The last coefficient is set to zero.
 """
 function derivative(a::Taylor1)
     coeffs = zero(a.coeffs)
-    @inbounds coeffs[1] = a[1]
     @inbounds for i = 1:a.order
         coeffs[i] = i*a[i]
     end
@@ -230,4 +229,69 @@ hessian!(hes::Array{T,2}, f::TaylorN{T}) where {T<:Number} =
     jacobian!(hes, gradient(f))
 
 
-## TODO: Integration...
+##Integration
+"""
+    integrate(a, r)
+
+Integrate the `a::HomogeneousPolynomial` with respect to the `r`-th
+variable. The returned `HomogeneousPolynomial` has no added constant of
+integration. If the order of a corresponds to `get_order()`, a zero
+`HomogeneousPolynomial` of 0-th order is returned.
+
+"""
+function integrate(a::HomogeneousPolynomial, r::Int)
+    @assert 1 ≤ r ≤ get_numvars()
+
+    order_max = get_order()
+    T = promote_type(eltype(a), eltype(a[1]/1))
+    a.order == order_max && return HomogeneousPolynomial(zero(T), 0)
+
+    @inbounds posTb = pos_table[a.order+2]
+    @inbounds num_coeffs = size_table[a.order+2]
+    coeffs = zeros(T, size_table[num_coeffs])
+    @inbounds num_coeffs = size_table[a.order+1]
+
+    @inbounds for i = 1:num_coeffs
+        iind = coeff_table[a.order+1][i]
+        n = iind[r]
+        n == order_max && continue
+        iind[r] += 1
+        kdic = in_base(get_order(), iind)
+        pos = posTb[kdic]
+        coeffs[pos] = a[i] / (n+1)
+        iind[r] -= 1
+    end
+
+    return HomogeneousPolynomial{T}(coeffs, a.order+1)
+end
+
+
+"""
+    integrate(a, r, [x0])
+
+Integrate the `a::TaylorN` series with respect to the `r`-th variable,
+where `x0` the integration constant and must be independent
+of the `r`-th variable; if `x0` is ommitted, it is taken as zero.
+"""
+function integrate(a::TaylorN, r::Int)
+    T = promote_type(eltype(a), eltype(a[0]/1))
+    order_max = min(get_order(), a.order+1)
+    coeffs = zeros(HomogeneousPolynomial{T}, order_max)
+
+    @inbounds for ord = 0:order_max-1
+        coeffs[ord+1] = integrate( a[ord], r)
+    end
+
+    return TaylorN(coeffs)
+end
+function integrate(a::TaylorN, r::Int, x0::TaylorN)
+    # Check constant of integration is independent of re
+    @assert derivative(x0, r) == 0.0 """
+    The integration constant ($x0) must be independent of the
+    $(_params_TaylorN_.variable_names[r]) variable"""
+
+    res = integrate(a, r)
+    return x0+res
+end
+integrate(a::TaylorN, r::Int, x0::NumberNotSeries) =
+    integrate(a,r,TaylorN(HomogeneousPolynomial([convert(eltype(a),x0)], 0)))

test/manyvariables.jl

@@ -147,6 +147,18 @@ end
     vr = rand(2)
     @test hp(vr) == evaluate(hp, vr)
 
+    p = (xT-yT)^6
+    @test integrate(derivative(p, 1), 1, yT^6) == p
+    @test derivative(integrate(p, 2), 2) == p
+    @test derivative(TaylorN(1.0)) == 0.0
+    @test integrate(TaylorN(6.0), 1) == 6xT
+    @test integrate(TaylorN(0.0), 2) == 0.0
+    @test integrate(TaylorN(0.0), 2, xT) == xT
+    @test integrate(xT^17, 2) == 0.0
+    @test integrate(xT^17, 1, yT) == yT
+    @test integrate(xT^17, 1, 2.0) == 2.0
+    @test_throws AssertionError integrate(xT, 1, xT)
+
     @test derivative(2xT*yT^2,1) == 2yT^2
     @test xT*xT^3 == xT^4
     txy = 1.0 + xT*yT - 0.5*xT^2*yT + (1/3)*xT^3*yT + 0.5*xT^2*yT^2