2022-10-24 Galerkin methods
Contents
2022-10-24 Galerkin methods#
Last time#
Algebraic multigrid
Smoothed aggregation
Multigrid in PETSc
Projects
Today#
Questions about projects
Galerkin methods
Weak forms
\(L^2\) Projection
Galerkin “optimality”
using Plots
default(linewidth=3)
using LinearAlgebra
using SparseArrays
function my_spy(A)
cmax = norm(vec(A), Inf)
s = max(1, ceil(120 / size(A, 1)))
spy(A, marker=(:square, s), c=:diverging_rainbow_bgymr_45_85_c67_n256, clims=(-cmax, cmax))
end
my_spy (generic function with 1 method)
Searching for projects#
GitHub codesearch#
Check Insights -> Contributors
Sustained versus short-term activity
Number of stakeholders
git shortlog -se --since=2016 | sort -n
Community signals
Code of Conduct
Contributing guidelines, pull request templates
Active review and mentoring in pull requests
Continuous integration
Actions tab, also report in pull requests
Pipelines for repositories on GitLab
Journals#
Journal of Open Source Software#
Papers on mature software with many users
Immature research software with a few users, but best practices and intent to grow
Geoscientific Model Development#
Blends software and methods work
Integration by parts#
One dimension#
it’s the product rule backwards
you can move the derivative to the other term; it’ll cost you a minus sign and a boundary term
Multiple dimensions#
Strong form#
Weak form#
multiply by a test function and integrate by parts
Variational: from minimization to weak forms#
Minimization#
Find \(u\in \mathcal V\) such that
Variations#
If \(f(x)\) is a minimum, then \(df = f'(x) dx = 0\) for all \(dx\).
So if \(u \in \mathcal V\) is a minimizer, then
for all \(du \in \mathcal V\).
It’s common to use \(v\) in place of \(du\).
Our discretization will reflect this structure
Finite elements#
Problem statement#
Find \(u \in \mathcal V\) such that
What do we need?#
A way to integrate
A space \(\mathcal V\)
Boundary conditions
Ciarlet (1978) definition#
A finite element is a triple \((K, P, N)\) where
\(K\) is a bounded domain with piecewise smooth boundary
\(P = \{p_j\}_1^k\) is a finite-dimensional function space over \(K\) (“prime basis”)
\(N = \{n_i\}_1^k\) is a basis for the dual space (“nodes”)
It is common for nodes to be pointwise evaluation
First: choose a single element \(K = \Omega\)#
Quadrature#
Our elements will be a line segment and we’ll define quadrature on the standard line segment, \(\hat K = (-1, 1)\). A quadrature is a set of points \(q_i\) and weights \(w_i\) such that
default(linewidth=4)
function quad_trapezoid(n)
q = LinRange(-1, 1, n)
w = 2 * ones(n) / (n - 1)
w[[1,end]] ./= 2
q, w
end
function plot_quad_accuracy(f, F, quad)
ns = 2 .^ (2:7)
exact = F(1) - F(-1)
function err(n)
q, w = quad(n)
w' * f.(q) - exact
end
errors = [abs.(err(n)) for n in ns]
plot(ns, errors, marker=:auto, label=quad)
plot!(ns, ns.^(-2), xscale=:log10, yscale=:log10)
end
plot_quad_accuracy (generic function with 1 method)
plot_quad_accuracy(exp, exp, quad_trapezoid)
plot_quad_accuracy(x -> 3*x^2, x -> x^3, quad_trapezoid)
Integration via polynomial interpolation#
Pick some points \(x_i\)
Evaluate the function there \(f_i = f(x_i)\)
Find the polynomial that interpolates the data
Integrate the polynomial
What order of accuracy can we expect?#
What degree polynomials can be integrate exactly?
Doing better: Gauss quadrature#
Suppose a polynomial on the interval \([-1,1]\) can be written as
where \(P_n(x)\) is the \(n\)th Legendre polnomials and both \(q(x)\) and \(r(x)\) are polynomials of maximum degree \(n-1\).
Why is \(\int_{-1}^1 P_n(x) q(x) = 0\)?
Can every polynomials of degree \(2n-1\) be written in the above form?
How many roots does \(P_n(x)\) have on the interval?
Can we choose points \(\{x_i\}\) such that the first term is 0?
If \(P_n(x_i) = 0\) for each \(x_i\), then we need only integrate \(r(x)\), which is done exactly by integrating its interpolating polynomial. How do we find these roots \(x_i\)?
Gauss-Legendre in code#
Solve for the points, compute the weights
Use a Newton solver to find the roots. You can use the recurrence to write a recurrence for the derivatives.
Create a Vandermonde matrix and extract the first row of the inverse or (using more structure) the derivatives at the quadrature points.
Use duality of polynomial roots and matrix eigenvalues.
A fascinating mathematical voyage (topic of graduate linear algebra class).
function gauss_legendre(n)
"""Gauss-Legendre integration using Golub-Welsch algorithm"""
beta = @. .5 / sqrt(1 - (2 * (1:n-1))^(-2))
T = diagm(-1 => beta, 1 => beta)
D, V = eigen(T)
w = V[1,:].^2 * 2
q = D
q, w
end
gauss_legendre (generic function with 1 method)
plot_quad_accuracy(exp, exp, gauss_legendre)
plot_quad_accuracy(x -> 11*x^10, x -> x^11, gauss_legendre)
Legendre polynomials#
function vander_legendre(x, k=nothing)
if isnothing(k)
k = length(x) # Square by default
end
m = length(x)
Q = ones(m, k)
Q[:, 2] = x
for n in 1:k-2
Q[:, n+2] = ((2*n + 1) * x .* Q[:, n+1] - n * Q[:, n]) / (n + 1)
end
Q
end
vander_legendre (generic function with 2 methods)
x = LinRange(-1, 1, 50)
P = vander_legendre(x, 6)
plot(x, P)
Prime to nodal basis#
A nodal basis \(\{ \phi_j(x) \}\) is one that satisfies \( n_i(\phi_j) = \delta_{ij} . \)
We write \(\phi_j\) in the prime basis by solving with the generalized Vandermonde matrix \(V_{ij} = \langle n_i, p_j \rangle\),
k = 10
xn = LinRange(-1, 1, k)
V = vander_legendre(xn)
xx = LinRange(-1, 1, 50)
Pxx = vander_legendre(xx, k)
Pn = Pxx / V
plot(xx, Pn)
scatter!(xn, one.(xn))
using FastGaussQuadrature
xn, _ = gausslobatto(k)
V = vander_legendre(xn)
Pn = Pxx / V
plot(xx, Pn)
scatter!(xn, one.(xn))
function vander_legendre_deriv(x, k=nothing)
if isnothing(k)
k = length(x) # Square by default
end
m = length(x)
Q = ones(m, k)
dQ = zeros(m, k)
Q[:, 2] = x
dQ[:, 2] .= 1
for n in 1:k-2
Q[:, n+2] = ((2*n + 1) * x .* Q[:, n+1] - n * Q[:, n]) / (n + 1)
dQ[:, n+2] = (2*n + 1) * Q[:,n+1] + dQ[:,n]
end
Q, dQ
end
vander_legendre_deriv (generic function with 2 methods)
Galerkin method: \(L^2\) projection#
A nice test problem that doesn’t require derivatives or boundary conditions: Find \(u \in \mathcal V_h\) such that
function febasis(P, Q)
x, _ = gausslobatto(P)
q, w = gausslegendre(Q)
Pk, _ = vander_legendre_deriv(x)
Bp, Dp = vander_legendre_deriv(q, P)
B = Bp / Pk
D = Dp / Pk
x, q, w, B, D
end
function L2_galerkin(P, Q, f)
x, q, w, B, _ = febasis(P, Q)
A = B' * diagm(w) * B
rhs = B' * diagm(w) * f.(q)
u = A \ rhs
x, u
end
L2_galerkin (generic function with 1 method)
f(x) = tanh(5*x)
P, Q = 8, 8
x, u = L2_galerkin(P, Q, f)
error = u - f.(x)
plot(x, u, marker=:auto, legend=:bottomright)
plot!(f, title="Galerkin error: $(norm(error))")
Embed the solution in the continuous space#
In thel last figure, the Galerkin solution was represented by a piecewise linear function, emphasizing its nodal values. Let’s make a nicer plot by multiplying by the corresponding nodal basis function.
P, Q = 10, 12
x, u = L2_galerkin(P, Q, f)
error = u - f.(x)
_, xx, _, BB, _ = febasis(P, 100)
plot(f, label="exact", title="Galerkin error: $(norm(error))",
xlim=(-1, 1), legend=:bottomright)
plot!(xx, BB * u, label="Galerkin solution")
scatter!(x, u, marker=:circle, color=:black, label="Nodal values")
ns = 2 .^ (2:8)
function L2_error(n)
x, u = L2_galerkin(n, n, f)
norm(u - f.(x))
end
plot(ns, L2_error.(ns), marker=:circle, yscale=:log10,
title="Convergence", xlabel="number of points", ylabel="error")
Galerkin error estimates#
We introduce the notation
bilinear (linear in each of its arguments)
symmetric: \(a(u, v) = a(v,u)\)
positive definite: \(a(u, u) > 0\) when \(u \ne 0\) thus defines an inner product on the function space \(V\).
We also introduce the \(L^2\) inner product
Galerkin optimality via energy norms#
We can also define the “energy norm” or \(a\)-norm,
In other words,
Observations#
The Galerkin method computes the exact solution any time it resides in the subspace \(V_h\).
The Galerkin method is automatically symmetric any time the weak form is symmetric.
The Galerkin method can be spectrally accurate, similar to the Chebyshev finite difference methods.
For a nonlinear problem, discretization and differentiation will commute.