recursive evaluation of derivatives (#5)

Project.toml

@@ -1,6 +1,6 @@
 name = "LegendrePolynomials"
 uuid = "3db4a2ba-fc88-11e8-3e01-49c72059a882"
-version = "0.3.1"
+version = "0.3.2"
 
 [deps]
 OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"

docs/Project.toml

@@ -1,10 +1,8 @@
 [deps]
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
-DualNumbers = "fa6b7ba4-c1ee-5f82-b5fc-ecf0adba8f74"
 HyperDualNumbers = "50ceba7f-c3ee-5a84-a6e8-3ad40456ec97"
 LegendrePolynomials = "3db4a2ba-fc88-11e8-3e01-49c72059a882"
 
 [compat]
 Documenter = "0.26"
-DualNumbers = "0.6"
 HyperDualNumbers = "4"

docs/src/derivatives.md

@@ -1,11 +1,31 @@
 # Derivatives of Legendre Polynomials
 
+## Analytical recursive approach
+
+The Bonnet's recursion formula 
+
+```math
+P_\ell(x) = \left((2\ell-1) x P_{\ell-1}(x) - (\ell-1)P_{\ell - 2}(x)\right)/\ell
+```
+
+may be differentiated an arbitrary number of times analytically to obtain recursion relations for higher derivatives: 
+
+```math
+\frac{d^n P_\ell(x)}{dx^n} = \frac{(2\ell-1)}{\ell} \left(x \frac{d^n P_{\ell-1}(x)}{dx^n} + 
+n \frac{d^{(n-1)} P_{\ell-1}(x)}{dx^{(n-1)}} \right) - \frac{(\ell-1)}{\ell} \frac{d^n P_{\ell-2}(x)}{dx^n}
+```
+
+This provides a simultaneous recursion relation in ``\ell`` as well as ``n``, solving which we may obtain derivatives up to any order. This is the approach used in this package to compute the derivatives of Legendre polynomials.
+
+## Automatic diferentiation
+
 The Julia automatic differentiation framework may be used to compute the derivatives of Legendre polynomials alongside their values. Since the defintions of the polynomials are completely general, they may be called with dual or hyperdual numbers as arguments to evaluate derivarives in one go. 
 We demonstrate one example of this using the package [`HyperDualNumbers.jl`](https://github.com/JuliaDiff/HyperDualNumbers.jl) v4:
 
 ```@meta
 DocTestSetup = quote
 	using LegendrePolynomials
+  using HyperDualNumbers
 end
 ```
 
@@ -84,51 +104,4 @@ Pl_dPl_d2Pl (generic function with 1 method)
 
 julia> Pl_dPl_d2Pl(0.5, lmax = 3)
 ([1.0, 0.5, -0.125, -0.4375], [0.0, 1.0, 1.5, 0.375], [0.0, 0.0, 3.0, 7.5])
-```
-
-# Analytical approach for higher derivatives
-
-Legendre polynomials satisfy the differential equation 
-
-```math
-\frac{d}{dx}\left[(1-x^2)\frac{d P_n}{dx} \right] + n(n+1) P_n(x) = 0
-```
-
-We may rearrange the terms to obtain 
-
-```math
-\frac{d^2 P_n}{dx^2} = \frac{1}{(1-x^2)}\left( 2x \frac{d P_n(x)}{dx} - n(n+1)P_n{x} \right)
-```
-
-We may therefore compute the second derivative from the function and its first derivative. Higher derivatives may further be computed in terms of the lower ones.
-
-We demonstrate the second-derivative computation using the package [`DualNumbers.jl`](https://github.com/JuliaDiff/DualNumbers.jl) v0.5: 
-
-```jldoctest dual
-julia> using DualNumbers
-
-julia> x = 0.5;
-
-julia> xd = Dual(x, one(x));
-
-julia> d2Pl(x, P, dP, n) = (2x * dP - n*(n+1) * P)/(1 - x^2);
-
-julia> function d2Pl(x, n)
-           xd = Dual(x, one(x))
-           y = Pl(xd, n)
-           P, dP = realpart(y), dualpart(y)
-           d2Pl(x, P, dP, n)
-       end;
-
-julia> d2Pl(x, 20)
-32.838787646905985
-```
-
-We may check that this matches the result obtained using `HyperDualNumbers`:
-
-```jldoctest hyperdual
-julia> ε₁ε₂part(Pl(xh, 20))
-32.838787646905985
-```
-
-Unfortunately at this point, higher derivatives need to be evaluated analytically and expresed in terms of lower derivatives.
+```

docs/src/index.md

@@ -7,7 +7,7 @@ end
 
 # Introduction
 
-Compute [Legendre polynomials](https://en.wikipedia.org/wiki/Legendre_polynomials) using a 3-term recursion relation (Bonnet’s recursion formula).
+Compute [Legendre polynomials](https://en.wikipedia.org/wiki/Legendre_polynomials) and their derivatives using a 3-term recursion relation (Bonnet’s recursion formula).
 
 ```math
 P_\ell(x) = \left((2\ell-1) x P_{\ell-1}(x) - (\ell-1)P_{\ell - 2}(x)\right)/\ell
@@ -19,10 +19,12 @@ Currently this package evaluates the standard polynomials that satisfy ``P_\ell(
 \int_{-1}^1 P_m(x) P_n(x) dx = \frac{2}{2n+1} \delta_{mn}.
 ```
 
-There are two main functions: 
+There are four main functions: 
 
 * [`Pl(x,l)`](@ref Pl): this evaluates the Legendre polynomial for a given degree `l` at the argument `x`. The argument needs to satisfy `-1 <= x <= 1`.
 * [`collectPl(x; lmax)`](@ref collectPl): this evaluates all the polynomials for `l` lying in `0:lmax` at the argument `x`. As before the argument needs to lie in the domain of validity. Functionally this is equivalent to `Pl.(x, 0:lmax)`, except `collectPl` evaluates the result in one pass, and is therefore faster. There is also the in-place version [`collectPl!`](@ref) that uses a pre-allocated array.
+* [`dnPl(x, l, n)`](@ref dnPl): this evaluates the ``n``-th derivative of the Legendre polynomial ``P_\ell(x)`` at the argument ``x``. The argument needs to satisfy `-1 <= x <= 1`.
+* [`collectdnPl(x; n, lmax)`](@ref collectdnPl): this evaluates the ``n``-th derivative of all the Legendre polynomials for `l = 0:lmax`. There is also an in-place version [`collectdnPl!`](@ref) that uses a pre-allocated array.
 
 # Quick Start
 
@@ -33,7 +35,14 @@ julia> Pl(0.5, 3)
 -0.4375
 ```
 
-Evaluate all the polynomials for `l` in `0:lmax` as 
+Evaluate the `n`th derivative for one `l` as `dnPl(x, l, n)`:
+
+```jldoctest
+julia> dnPl(0.5, 3, 2)
+7.5
+```
+
+Evaluate all the polynomials for `l` in `0:lmax` as `collectPl(x; lmax)`
 
 ```jldoctest
 julia> collectPl(0.5, lmax = 3)
@@ -44,6 +53,19 @@ julia> collectPl(0.5, lmax = 3)
  -0.4375
 ```
 
+Evaluate all the `n`th derivatives as `collectdnPl(x; lmax, n)`:
+
+```jldoctest
+julia> collectdnPl(0.5, lmax = 5, n = 3)
+6-element OffsetArray(::Array{Float64,1}, 0:5) with eltype Float64 with indices 0:5:
+  0.0
+  0.0
+  0.0
+ 15.0
+ 52.5
+ 65.625
+```
+
 # Increase precision
 
 The precision of the result may be changed by using arbitrary-precision types such as `BigFloat`. For example, using `Float64` arguments we obtain
@@ -69,6 +91,16 @@ julia> setprecision(300) do
 0.33333333333333333333333333333333333333333333333333333333333333333333333333333333333333333317
 ```
 
+This is particularly important to avoid overflow while computing high-order derivatives. For example:
+
+```jldoctest
+julia> dnPl(0.5, 300, 200) # Float64
+NaN
+
+julia> dnPl(big(1)/2, 300, 200) # BigFloat
+1.738632750542319394663553898425873258768856732308227932150592526951212145232716e+499
+```
+
 # Reference
 
 ```@autodocs