jagot/AnalyticBSplines.jl

AnalyticBSplines.jl

Simple library for defining and evaluating B-splines analytically. Companion package to BSplinesQuasi.jl.

Usage

We first define a knot set (using rationals such that polynomial integrals/derivatives will be exact).

using AnalyticBSplines
using BSplinesQuasi # For knot-sets

t = LinearKnotSet(3, 0:1//2:2)

LinearKnotSet{Rational{Int64}} of order k=3 on [0//1,2//1] (4 intervals)

We then generate the basis

B = gen_basis(t)
B[end]

6-element Array{PPoly{Rational{Int64}},1}:
 Piecewise polynomial{Rational{Int64}}
 Poly(1//1 - 4//1*x + 4//1*x^2) on Interval[0//1,1//2]                                                                                                                 
 Piecewise polynomial{Rational{Int64}}
 Poly(4//1*x - 6//1*x^2) on Interval[0//1,1//2]
 Poly(2//1 - 4//1*x + 2//1*x^2) on Interval[1//2,1//1]                                                                 
 Piecewise polynomial{Rational{Int64}}
 Poly(2//1*x^2) on Interval[0//1,1//2]
 Poly(-3//2 + 6//1*x - 4//1*x^2) on Interval[1//2,1//1]
 Poly(9//2 - 6//1*x + 2//1*x^2) on Interval[1//1,3//2]                  
 Piecewise polynomial{Rational{Int64}}
 Poly(1//2 - 2//1*x + 2//1*x^2) on Interval[1//2,1//1]
 Poly(-11//2 + 10//1*x - 4//1*x^2) on Interval[1//1,3//2]
 Poly(8//1 - 8//1*x + 2//1*x^2) on Interval[3//2,2//1]
 Piecewise polynomial{Rational{Int64}}
 Poly(2//1 - 4//1*x + 2//1*x^2) on Interval[1//1,3//2]
 Poly(-16//1 + 20//1*x - 6//1*x^2) on Interval[3//2,2//1]                                                       
 Piecewise polynomial{Rational{Int64}}
 Poly(9//1 - 12//1*x + 4//1*x^2) on Interval[3//2,2//1]                                                                                                                

We can then easily generate the overlap matrix (or any other polynomial scalar operator):

scalar_op(B, t)

6×6 Array{Rational{Int64},2}:
 1//10   7//120   1//120   0//1    0//1    0//1  
 7//120  1//6     5//48    1//240  0//1    0//1  
 1//120  5//48   11//40   13//120  1//240  0//1  
 0//1    1//240  13//120  11//40   5//48   1//120
 0//1    0//1     1//240   5//48   1//6    7//120
 0//1    0//1     0//1     1//120  7//120  1//10 

derop(B, t, 1)

6×6 Array{Any,2}:
 -1//2    5//12   1//12    0       0      0   
 -5//12   0//1    3//8    1//24    0      0   
 -1//12  -3//8    0//1    5//12   1//24   0   
   0     -1//24  -5//12   0//1    3//8   1//12
   0       0     -1//24  -3//8    0//1   5//12
   0       0       0     -1//12  -5//12  1//2 

derop(B, t, 2)

6×6 Array{Any,2}:
 4//3  -2//1   2//3    0      0     0  
 2//1  -8//3   1//3   1//3    0     0  
 2//3   1//3  -2//1   2//3   1//3   0  
  0     1//3   2//3  -2//1   1//3  2//3
  0      0     1//3   1//3  -8//3  2//1
  0      0      0     2//3  -2//1  4//3

Note that the second-derivative operator is not symmetric towards the end-points of the grid; it will only be so if

Bᵢ(b)Bⱼ′(b) - Bᵢ(a)Bᵢ′(a) - Bᵢ′(b)Bⱼ(b) + Bᵢ′(a)Bⱼ(a) = 0,

where [a,b] are the end-points. These boundary conditions can be satisfied by e.g. removing the first and last BSpline (this amounts to Dirichlet0 conditions). We can accomplish this thusly:

B′ = gen_basis(t, 0, 0)
derop(B′, t, 2)

6×6 Array{Any,2}:
 0    0      0      0      0    0
 0  -8//3   1//3   1//3    0    0
 0   1//3  -2//1   2//3   1//3  0
 0   1//3   2//3  -2//1   1//3  0
 0    0     1//3   1//3  -8//3  0
 0    0      0      0      0    0

Finally, the k:th derivative operator of a basis of order k, is zero:

derop(B, t, 3)

6×6 Array{Any,2}:
 0  0  0  0  0  0
 0  0  0  0  0  0
 0  0  0  0  0  0
 0  0  0  0  0  0
 0  0  0  0  0  0
 0  0  0  0  0  0