Commit acba7fd1 authored by Juan Ignacio Polanco's avatar Juan Ignacio Polanco

Compute Galerkin matrix ⟨bi, bj⟩

parent ad0821fc
......@@ -5,6 +5,7 @@ version = "0.1.0"
[deps]
BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
......
......@@ -20,8 +20,10 @@ export knots, order, coefficients
export augment_knots
export evaluate_bspline, evaluate_bspline!
export integral
export galerkin_matrix, galerkin_matrix!
using BandedMatrices
using FastGaussQuadrature
using Reexport
using LinearAlgebra: Symmetric
using StaticArrays: MVector
......
......@@ -35,7 +35,12 @@ function galerkin_matrix!(S::Symmetric, B::BSplineBasis)
A = parent(S)
k = order(B)
h = (k + 1) >> 2 # k/2 if k is even
t = knots(B)
h = (k + 1) ÷ 2 # k/2 if k is even
T = eltype(S)
# Quadrature information (weights, nodes).
quad = _quadrature_prod(k)
# Upper part: j >= i
for j = 1:N
......@@ -43,10 +48,53 @@ function galerkin_matrix!(S::Symmetric, B::BSplineBasis)
# In other words, we know that S[i, j] = 0 outside the chosen interval.
istart = clamp(j - h, 1, N)
bj = BSpline(B, j)
tj = j:(j + k) # support of b_j (knot indices)
for i = istart:j
ti = i:(i + k) # support of b_i
t_inds = intersect(ti, tj) # common support of b_i and b_j
@assert !isempty(t_inds) # there is a common support (the B-splines see each other)
@assert length(t_inds) == k + 1 - (j - i)
bi = BSpline(B, i)
A[i, j] = _integrate_prod(bi, bj, t, t_inds, quad)
end
end
S
end
# Generate quadrature information for B-spline product.
# Returns weights and nodes for integration in [-1, 1].
#
# See https://en.wikipedia.org/wiki/Gaussian_quadrature.
#
# Some notes:
#
# - On each interval between two neighbouring knots, each B-spline is a
# polynomial of degree (k - 1). Hence, the product of two B-splines has degree
# (2k - 2).
#
# - The Gauss--Legendre quadrature rule is exact for polynomials of degree
# <= (2n - 1), where n is the number of nodes and weights.
#
# - Conclusion: on each knot interval, `k` nodes should be enough to get an
# exact integral. (I verified that results don't change when using more than
# `k` nodes.)
_quadrature_prod(k) = gausslegendre(k)
# Integrate product of functions over the subintervals t[inds].
function _integrate_prod(f, g, t, inds, (x, w))
int = 0.0 # compute stuff in Float64, regardless of type wanted by the caller
N = length(w) # number of weights / nodes
for i in inds[2:end]
# Integrate in [t[i - 1], t[i]].
# See https://en.wikipedia.org/wiki/Gaussian_quadrature#Change_of_interval
a, b = t[i - 1], t[i]
α = (b - a) / 2
β = (a + b) / 2
for n = 1:N
y = α * x[n] + β
int += α * w[n] * f(y) * g(y)
end
end
int
end
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment