You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert
Legendre polynomials and their derivatives using automatic differentiation
7
+
Iterative computation of Legendre Polynomials
8
8
9
9
# Getting Started
10
10
11
-
## Prerequisites
12
-
13
-
The package depends on `HyperDualNumbers` and `OffsetArrays`.
14
-
15
11
## Installing
16
12
17
13
To install the package, run
@@ -20,134 +16,29 @@ To install the package, run
20
16
] add LegendrePolynomials
21
17
```
22
18
23
-
## Using
24
-
25
-
This package computes Legendre polynomials and their first and second derivatives in one pass using hyperdual numbers. All the functions require an argument `x` that satisfies `-1<=x<=1`, and either an integer `l` that indicates the degree of the polynomial, or a possibly optional keyword argument `lmax` that indicates the maximum degree of the Legendre polynomials that will be computed. The polynomials are computed iteratively using a three-term recurrence relation, so all the values in the range will have to be computed.
26
-
27
-
There are two classes of functions, one that returns the value of the Legendre polynomial and its derivatives for a single `l`, and one that returns the values for all `l` in `0:lmax` for a specified `lmax` as an array. The arrays returned will have the Legendre polynomials and their derivatives stored along columns.
28
-
29
-
There are three quantities that can be computed: `Pl`, `dPl` and `d2Pl`. It's possible to compute them individually or as combinations of two or three at once. There are allocating as well as non-allocating methods that do the same.
30
-
31
-
### Individual degrees
19
+
## Quick start
32
20
33
-
The way to compute Legendre polynomials of degree `l` and argument `x`is through the function `Pl(x,l)`. This is the general signature followed by the other functions as well that compute the derivatives for a single degree. The functions are relatively accurate and performant till degrees of millions, although a thorough test of accuracy has not been carried out (verified roughly with results from Mathematica for some degrees).
21
+
To compute the Legendre polynomials for a given argument `x`and a degree `l`, use `Pl(x,l)`:
34
22
35
23
```julia
36
-
julia>Pl(0.5,10)
37
-
-0.18822860717773438
38
-
39
-
julia>@btimePl(0.5,1_000_000)
40
-
13.022 ms (0 allocations:0 bytes)
41
-
-0.0006062610545162491
24
+
julia>Pl(0.5, 20)
25
+
-0.04835838106737356
42
26
```
43
27
44
-
The accuracy can be increased by using Arbitrary Precision Arithmetic, through the use of
45
-
`BigInt` and `BigFloat`. This comes at the expense of performance though.
28
+
To compute all the polynomials for `0 <= l <= lmax`, use
The way to compute the derivatives is through `dPl(x,l)` and `d2Pl(x,l)`, for example
58
-
59
-
```julia
60
-
julia>dPl(0.5,5)
61
-
-2.2265625000000004
62
-
63
-
julia>d2Pl(0.5,5)
64
-
-6.5625
65
-
```
66
-
67
-
Combinations of these three can be computed as well, for example
68
-
69
-
```julia
70
-
julia>Pl_dPl(0.5,5)
71
-
(0.08984375, -2.2265625000000004)
72
-
73
-
julia>Pl_d2Pl(0.5,5)
74
-
(0.08984375, -6.5625)
75
-
76
-
julia>dPl_d2Pl(0.5,5)
77
-
(-2.2265625000000004, -6.5625)
78
-
79
-
julia>Pl_dPl_d2Pl(0.5,5)
80
-
(0.08984375, -2.2265625000000004, -6.5625)
81
-
```
82
-
83
-
### All degrees up to a cutoff
84
-
85
-
The second class of methods return all the polynomials up to a cutoff degree `lmax`. They are returned as `OffsetArrays` that have 0-based indexing, keeping in mind that the polynomials start from `l=0`.
86
-
87
-
The polynomials and their derivatives can be computed in general by calling the function `P(x;lmax)`, where `P` has to be chosen appropriately as necessary. The keyword argument has to be specified in the allocating functions, whereas it may be inferred from the array in the non-allocating versions.
88
-
89
-
An example of the allocating functions are
90
-
91
-
```julia
92
-
julia>Pl(0.5,lmax=3)
93
-
4-element OffsetArray(::Array{Float64,1}, 0:3) with eltype Float64 with indices 0:3:
94
-
1.0
95
-
0.5
96
-
-0.125
31
+
julia>collectPl(0.5, lmax =5)
32
+
6-element OffsetArray(::Array{Float64,1}, 0:5) with eltype Float64 with indices 0:5:
33
+
1.0
34
+
0.5
35
+
-0.125
97
36
-0.4375
98
-
99
-
julia>dPl(0.5,lmax=3)
100
-
4-element OffsetArray(::Array{Float64,1}, 0:3) with eltype Float64 with indices 0:3:
101
-
0.0
102
-
1.0
103
-
1.5
104
-
0.3750000000000001
105
-
106
-
julia>d2Pl(0.5,lmax=3)
107
-
4-element OffsetArray(::Array{Float64,1}, 0:3) with eltype Float64 with indices 0:3:
108
-
0.0
109
-
0.0
110
-
3.0
111
-
7.5
112
-
```
113
-
114
-
It's possible to compute combinations analogous to the ones seen earlier, for example
This returns a 3-tuple of `OffsetArrays``Pl`, `dPl` and `d2Pl`.
122
-
123
-
There are non-allocating functions as well that can be called using a pre-allocated array. As an example
124
-
125
-
```julia
126
-
julia> P=zeros(0:3,0:2)
127
-
4×3OffsetArray(::Array{Float64,2}, 0:3, 0:2) with eltype Float64 with indices 0:3×0:2:
128
-
0.00.00.0
129
-
0.00.00.0
130
-
0.00.00.0
131
-
0.00.00.0
132
-
133
-
julia>Pl_dPl_d2Pl!(P,0.5)
134
-
4×3OffsetArray(::Array{Float64,2}, 0:3, 0:2) with eltype Float64 with indices 0:3×0:2:
135
-
1.00.00.0
136
-
0.51.00.0
137
-
-0.1251.53.0
138
-
-0.43750.3757.5
139
-
140
-
julia>fill!(P,0);
141
-
142
-
julia>dPl!(P,0.5)
143
-
4×3OffsetArray(::Array{Float64,2}, 0:3, 0:2) with eltype Float64 with indices 0:3×0:2:
144
-
0.00.00.0
145
-
0.01.00.0
146
-
0.01.50.0
147
-
0.00.3750.0
37
+
-0.2890625
38
+
0.08984375
148
39
```
149
40
150
-
Note that the column number that will be populated depends on the order of the derivative, assuming 0-based indexing. Therefore `Pl!` will fill column 0, `dPl!` will fill column 1 and `d2Pl!` will fill column 2. Combinations of these will fill multiple columns as expected.
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.
4
+
We demonstrate one example of this using the package [`HyperDualNumbers.jl`](https://github.com/JuliaDiff/HyperDualNumbers.jl) v4:
5
+
6
+
```@meta
7
+
DocTestSetup = quote
8
+
using LegendrePolynomials
9
+
end
10
+
```
11
+
12
+
```jldoctest hyperdual
13
+
julia> x = 0.5;
14
+
15
+
julia> Pl(x, 3)
16
+
-0.4375
17
+
18
+
julia> using HyperDualNumbers
19
+
20
+
julia> xh = Hyper(x, one(x), one(x), zero(x));
21
+
22
+
julia> p = Pl(xh, 3)
23
+
-0.4375 + 0.375ε₁ + 0.375ε₂ + 7.5ε₁ε₂
24
+
```
25
+
26
+
The Legendre polynomial ``P_\ell(x)`` may be obtained using
27
+
28
+
```jldoctest hyperdual
29
+
julia> realpart(p)
30
+
-0.4375
31
+
```
32
+
33
+
The first derivative ``dP_\ell(x)/dx`` may be obtained as
34
+
35
+
```jldoctest hyperdual
36
+
julia> ε₁part(p)
37
+
0.375
38
+
```
39
+
40
+
The second derivative ``d^2P_\ell(x)/dx^2`` may be obtained using
41
+
42
+
```jldoctest hyperdual
43
+
julia> ε₁ε₂part(p)
44
+
7.5
45
+
```
46
+
47
+
Something similar may also be evaluated for all `l` iteratively. For example, the function `collectPl` evaluates Legendre polynomials for the `HyperDualNumber` argument for a range of `l`:
48
+
49
+
```jldoctest hyperdual
50
+
julia> collectPl(xh, lmax=4)
51
+
5-element OffsetArray(::Array{Hyper{Float64},1}, 0:4) with eltype Hyper{Float64} with indices 0:4:
52
+
1.0 + 0.0ε₁ + 0.0ε₂ + 0.0ε₁ε₂
53
+
0.5 + 1.0ε₁ + 1.0ε₂ + 0.0ε₁ε₂
54
+
-0.125 + 1.5ε₁ + 1.5ε₂ + 3.0ε₁ε₂
55
+
-0.4375 + 0.375ε₁ + 0.375ε₂ + 7.5ε₁ε₂
56
+
-0.2890625 - 1.5625ε₁ - 1.5625ε₂ + 5.625ε₁ε₂
57
+
```
58
+
59
+
We may extract the first derivatives by broadcasting the function `ε₁part` on the array as:
60
+
61
+
```jldoctest hyperdual
62
+
julia> ε₁part.(collectPl(xh, lmax=4))
63
+
5-element OffsetArray(::Array{Float64,1}, 0:4) with eltype Float64 with indices 0:4:
64
+
0.0
65
+
1.0
66
+
1.5
67
+
0.375
68
+
-1.5625
69
+
```
70
+
71
+
Similarly the function `ε₁ε₂part` may be used to obtain the second derivatives.
72
+
73
+
Several convenience functions to compute the derivatives of Legendre polynomials were available in `LegendrePolynomials` v0.2, but have been removed in v0.3. The users are encouraged to implement convenience functions to extract the derivatives as necessary. As an exmaple, we may compute the polynomials and their first and second derivatives together as
*[`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`.
25
+
*[`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.
26
+
27
+
# Quick Start
28
+
29
+
Evaluate the Legendre polynomial for one `l` at an argument`x` as `Pl(x, l)`:
30
+
31
+
```jldoctest
32
+
julia> Pl(0.5, 3)
33
+
-0.4375
34
+
```
35
+
36
+
Evaluate all the polynomials for `l` in `0:lmax` as
37
+
38
+
```jldoctest
39
+
julia> collectPl(0.5, lmax = 3)
40
+
4-element OffsetArray(::Array{Float64,1}, 0:3) with eltype Float64 with indices 0:3:
41
+
1.0
42
+
0.5
43
+
-0.125
44
+
-0.4375
45
+
```
46
+
47
+
# Increase precision
48
+
49
+
The precision of the result may be changed by using arbitrary-precision types such as `BigFloat`. For example, using `Float64` arguments we obtain
0 commit comments