RFC: Return n-th derivative of a Taylor1 as a Taylor1 (#137)

* Add derivative method
derivative(n,a) returns a Taylor1 variable; renamed existing derivative(n,a)
method as derivativeval

* Export derivativeval; add tests

* Fix & add derivative tests

* derivative(n,a) returns value; derivative(a,n) returns polynomial

* Improve performance with in-place operations with suggestions by @lbenet

* Fix docstrings

* Fix tests for julia v0.7-DEV

* Remove comment; add derivative test; add jacobian test for mixtures

src/calculus.jl

@@ -14,11 +14,62 @@ Return the `Taylor1` polynomial of the differential of `a::Taylor1`.
 The last coefficient is set to zero.
 """
 function derivative(a::Taylor1)
-    coeffs = zero(a.coeffs)
-    @inbounds for i = 1:a.order
-        coeffs[i] = i*a[i]
+    res = zero(a)
+    @inbounds for ord = 0:a.order-1
+        derivative!(res, a, ord)
+    end
+    return res
+end
+
+"""
+    derivative!(res, a) --> nothing
+
+In-place version of `derivative`. Compute the `Taylor1` polynomial of the
+differential of `a::Taylor1` and save it into `res`. The last coefficient is
+set to zero.
+"""
+function derivative!(res::Taylor1, a::Taylor1)
+    @inbounds for ord = 0:a.order-1
+        derivative!(res, a, ord)
+    end
+    res[a.order] = zero(a[0])
+    nothing
+end
+
+doc"""
+    derivative!(p, a, k) --> nothing
+
+Update in-place the `k-th` expansion coefficient `p[k]` of `p = derivative(a)`
+for both `p` and `a` `Taylor1`.
+
+The coefficients are given by
+
+```math
+\begin{equation*}
+p_k = (k+1)a_{k+1}.
+\end{equation*}
+```
+
+"""
+derivative!(p::Taylor1, a::Taylor1, k::Int) = (p[k] = (k+1)*a[k+1])
+
+"""
+    derivative(a, n)
+
+Compute recursively the `Taylor1` polynomial of the n-th derivative of
+`a::Taylor1`.
+"""
+function derivative(a::Taylor1{T}, n::Int) where {T <: Number}
+    @assert a.order ≥ n ≥ 0
+    if n==0
+        return a
+    else
+        res = deepcopy(a)
+        for i in 1:n
+            derivative!(res, res)
+        end
+        return res
     end
-    return Taylor1(coeffs, a.order)
 end
 
 """

test/manyvariables.jl

@@ -250,7 +250,11 @@ end
     @test q[end-1:end] == pol[end-1:end]
 
 
-    @test_throws DomainError yT^(-2)
+    if VERSION < v"0.7.0-DEV"
+        @test_throws DomainError yT^(-2)
+    else
+        @test_throws AssertionError yT^(-2)
+    end
     @test_throws DomainError yT^(-2.0)
     @test (1+xT)^(3//2) == ((1+xT)^0.5)^3
     @test real(xH) == xH

test/mixtures.jl

@@ -259,4 +259,8 @@ end
     @test p1N ≈ p1N
     @test p1N ≈ q1N
 
+    Pv = [rndTN(get_order(), 3), rndTN(get_order(), 3)]
+    Qv = convert.(Taylor1{TaylorN{Float64}}, Pv)
+
+    @test jacobian(Pv) == jacobian(Qv)
 end

test/onevariable.jl

@@ -301,9 +301,25 @@ end
     @test evaluate(v) == vv
     @test evaluate(v, complex(0.0,0.2)) ==
         [complex(0.0,sinh(0.2)),complex(cos(0.2),sin(-0.2))]
+
+    @test derivative(exp(ta(1.0)), 0) == exp(ta(1.0))
+    expected_result_approx = Taylor1(convert(Vector{Float64},exp(ta(1.0))[0:10]))
+    @test derivative(exp(ta(1.0)), 5) ≈ expected_result_approx atol=eps() rtol=0.0
+    expected_result_approx = Taylor1(convert(Vector{Float64},exp(ta(1.0pi))[0:12]),15)
+    @test derivative(exp(ta(1.0pi)), 3) ≈ expected_result_approx atol=eps(16.0) rtol=0.0
+    expected_result_approx = Taylor1(convert(Vector{Float64},exp(ta(1.0pi))[0:5]),15)
+    @test derivative(exp(ta(1.0pi)), 10) ≈ expected_result_approx atol=eps(64.0) rtol=0.0
+
+
+
+    @test derivative(exp(ta(1.0)), 5)() == exp(1.0)
+    @test derivative(exp(ta(1.0pi)), 3)() == exp(1.0pi)
+    @test isapprox(derivative(exp(ta(1.0pi)), 10)() , exp(1.0pi) )
+
     @test derivative(5, exp(ta(1.0))) == exp(1.0)
     @test derivative(3, exp(ta(1.0pi))) == exp(1.0pi)
     @test isapprox(derivative(10, exp(ta(1.0pi))) , exp(1.0pi) )
+
     @test integrate(derivative(exp(t)),1) == exp(t)
     @test integrate(cos(t)) == sin(t)
 
@@ -318,7 +334,11 @@ end
     @test_throws ArgumentError 1/t
     @test_throws ArgumentError zt/zt
     @test_throws ArgumentError t^1.5
-    @test_throws DomainError t^(-2)
+    if VERSION < v"0.7.0-DEV"
+        @test_throws DomainError t^(-2)
+    else
+        @test_throws ArgumentError t^(-2)
+    end
     @test_throws ArgumentError sqrt(t)
     @test_throws ArgumentError log(t)
     @test_throws ArgumentError cos(t)/sin(t)