Boundary Value Problem Methods

On a uniform grid \(x_i=a+ih\), \(h=(b-a)/N\), replace derivatives by centred differences:

This is the three-point stencil \(\big[\,{-}\dfrac{p}{h^2}-\dfrac{b}{2h},\ \dfrac{2p}{h^2}+c,\ {-}\dfrac{p}{h^2}+\dfrac{b}{2h}\,\big]\) at indices \(i{-}1,i,i{+}1\), giving a tridiagonal system solved by the Thomas algorithm. Dirichlet values are substituted directly; Neumann uses a second-order ghost node so the whole scheme stays \(O(h^2)\).

Postprocessing. The centred scheme has a regular error expansion in even powers of \(h\), so Richardson extrapolation across a grid sequence (\(h,\ h/2,\dots\)) cancels the leading \(h^2\) term and lifts the order to \(O(h^4)\) and beyond - the same idea as Romberg. Deferred correction gives a comparable gain on a single grid.

Continuous piecewise-polynomial Lagrange elements of degree \(k\) (selectable P\(_1\) linear or P\(_2\) quadratic). On a grid \(x_i=a+ih\) the nodal basis functions are defined piecewise, element by element.

P\(_1\) (linear). The nodal basis is the hat function

P\(_2\) (quadratic). Each element \([x_i,x_{i+1}]\) carries three shape functions in the local coordinate \(\xi=\dfrac{x-x_i}{h}\in[0,1]\) - two vertex functions and one interior "bubble" at the midpoint \(x_{i+1/2}\):

Glued across the shared node, the global quadratic vertex basis \(\varphi_i\) and midpoint basis \(\varphi_{i+1/2}\) are

Each element \([x_i,x_{i+1}]\) contributes a local \((k{+}1)\times(k{+}1)\) matrix - for P\(_1\),

assembled into the global stiffness matrix via the local-to-global index map; the load \(\int f\varphi_i\) and all element integrals are evaluated by Gauss quadrature. Neumann data enters naturally through the boundary term; Dirichlet is imposed by elimination. Increasing the degree raises the order to \(O(h^{k+1})\): P\(_1\) is second order, P\(_2\) third order. For pure diffusion the P\(_1\) nodal values are exact - a classical Galerkin superconvergence result.

Postprocessing. FEM is the classical setting for superconvergence recovery. The nodal values are superconvergent (exact here for pure diffusion), and the derivative \(u_h'\) is superconvergent at the element Gauss points: averaging/recovering it there - gradient recovery, Zienkiewicz–Zhu / superconvergent patch recovery - yields a flux one order more accurate than the raw \(u_h'\). Functional outputs \(J(u)\) can be postprocessed to roughly double order through an adjoint (Babuška–Miller).

Sample \(u\) at the Chebyshev–Gauss–Lobatto nodes \(x_j=\cos\!\big(\dfrac{j\pi}{N}\big),\ j=0,\dots,N\), and enforce the strong equation there using the Chebyshev differentiation matrix \(D\), whose entries are

with \(c_0=c_N=2\) and \(c_1=\dots=c_{N-1}=1\). It is a dense \((N{+}1)\times(N{+}1)\) matrix

for example at \(N=2\) (nodes \(x_0=1,\ x_1=0,\ x_2=-1\)),

The strong equation is collocated as

a dense system. Boundary rows are overwritten to enforce the conditions (an identity row for Dirichlet, a row of \(D\) for Neumann). For smooth solutions the error decays faster than any power of \(h\) - a handful of points reaches machine precision.

Postprocessing. On the smooth solutions here there is little to recover - the method is already spectrally accurate. Postprocessing matters for non-smooth data, where spectral filtering or Gegenbauer reprojection restores high-order accuracy away from a discontinuity and tames the Gibbs oscillations that otherwise pollute the whole domain. The spectral coefficients also give an accurate derivative directly, by applying \(D\) to the computed \(\mathbf{u}\).

Piecewise polynomials of degree \(k\), discontinuous across element interfaces. The symmetric interior-penalty (SIPG) bilinear form couples neighbours through averages \(\{\cdot\}\), jumps \([\cdot]\), and a penalty \(\sigma\):

The penalty must be large enough for coercivity; too small and the method loses stability (more so at higher \(k\)), too large and it is over-stiff. The matrix is block tridiagonal. Dirichlet is imposed weakly (Nitsche).

Postprocessing. DG is the poster child for postprocessing. A SIAC filter (smoothness-increasing accuracy-conserving) convolves the degree-\(k\) solution against a B-spline kernel and lifts the order from \(k+1\) to \(2k+1\) (linear problems), cashing in DG's superconvergence in the negative-order norm - and also smooths the inter-element jumps. It uses the same B-spline machinery as the spline method above.

A Galerkin method whose trial space is the degree-\(d\) spline space \(S_d^{\,d-1}\) on a clamped uniform knot vector - globally \(C^{d-1}\), and interpolatory at the endpoints so Dirichlet data pins the first and last coefficients directly. On each knot span (element) \(I=[x_\ell,x_{\ell+1}]\) the polynomial pieces are written in Bernstein–Bézier form (Lai & Schumaker, Spline Functions on Triangulations). Using the barycentric coordinates of \(I\),

the degree-\(d\) Bernstein basis polynomials and the B-form of a piece \(s|_I\) are

with B-coefficients (Bézier ordinates) \(c_{ij}\) carried at the domain points \(\xi_{ij}=\dfrac{i\,x_\ell+j\,x_{\ell+1}}{d}\). Global \(C^{d-1}\) smoothness is enforced through the Lai–Schumaker smoothness conditions tying the \(c_{ij}\) of adjacent spans, which reduces the per-span B-form to the clamped B-spline basis. The weak form is assembled by Gauss quadrature over each span, giving a banded matrix of bandwidth \(2d+1\); the extra smoothness delivers order \(d+1\) with far fewer degrees of freedom than P\(_1\) - the isogeometric idea behind spline finite elements.

Postprocessing. The globally \(C^{d-1}\) basis makes the derivatives \(s',s'',\dots\) directly available and continuous - no gradient recovery needed - and they are read off cheaply from the B-coefficients by de Boor's derivative formula. Isogeometric superconvergence and k-refinement give further accuracy from the same spline space.

Recast the BVP as a first-order system \(y_1=u,\ y_2=u'\),

and integrate it as an initial value problem with RK4, treating the unknown boundary datum \(s\) (the missing \(u'(a)\) or \(u(a)\)) as a parameter. The residual \(F(s)=u(b;s)-\beta\) is driven to zero by the secant method. Because the ODE is linear, \(F\) is affine in \(s\) and two shots locate the root exactly. Accuracy is that of the integrator, \(O(h^4)\).

Postprocessing. Two handles come from the integrator. A continuous (dense-output) extension of RK4 gives a smooth \(u(x)\) between steps - values and derivatives anywhere, not just at the grid nodes. And because RK4 has a regular \(h^4\) error expansion, Richardson extrapolation over halved step sizes raises the global order further.