OpenPFC / Operators and their solutions
A zero weight is a constraint
Solve the fifth positive Airy eigenmode without dividing by the zero coordinate at the centre of the grid.
Download runnable lesson (.zip)
Includes the model, inputs, required headers and a short run guide.
Program 22 of Trefethen (2000) solves a generalized eigenproblem:
\[ \begin{gathered} u''=\lambda xu, \\ u(-1)=u(1)=0. \end{gathered} \]
The weight \(x\) changes sign and vanishes at the centre. This is not a positive-definite mass matrix. Dividing every equation by \(x\) would be especially dangerous at the central Chebyshev point, whose floating-point coordinate is merely close to zero.
Eliminate the constraint, not a tiny denominator
The centre equation is \((D^2u)_c=0\). Write the other interior indices as O:
\[ \begin{gathered} u_c=-A_{cc}^{-1}A_{cO}u_O, \\ (A_{OO}-A_{Oc}A_{cc}^{-1}A_{cO})u_O=\lambda\,\mathrm{diag}(x_O)u_O. \end{gathered} \]
The remaining weights are nonzero. We solve this reduced finite eigenproblem, recover the central value, and normalize the fifth positive mode to \(\operatorname{Ai}(0)=1/[3^{2/3}\Gamma(2/3)]\).
for (int i = 0; i < n - 2; ++i)
for (int j = 0; j < n - 2; ++j)
a[i][j] =
(d2[free[i]][free[j]] -
d2[free[i]][center] * d2[center][free[j]] / d2[center][center]) /
x[free[i]];
const auto pairs = eigen_reference::solve(eigen_reference::complex_matrix(a));
std::vector<eigen_reference::Eigenpair> positive;
for (const auto &pair : pairs)
if (pair.value.real() > 0) positive.push_back(pair);
if (positive.size() < 5)
throw std::runtime_error("Fewer than five positive eigenvalues");
const auto &fifth = positive[4];
polynomial::Vector u(n + 1);
for (int j = 0; j < n - 2; ++j) u[free[j]] = fifth.vector[j].real();
for (int j : free) u[center] -= d2[center][j] * u[j] / d2[center][center];
const double ai0 = 1 / (std::pow(3.0, 2.0 / 3) * std::tgamma(2.0 / 3));
const double scaling = ai0 / u[center];
for (double &value : u) value *= scaling;This removes the infinite generalized eigenvalue associated with the zero weight; it does not discard a finite physical mode. Both the Schur reduction and the dense QR are example-owned mathematics. They use the same Chebyshev reference developed earlier, without claiming a public generalized-eigenvalue API in OpenPFC.
Separate eigenvalue convergence from residual size
openpfc run model.cpp input.json --source-root .. --launcher=noneThe defaults reproduce degrees 12,24,36,48. The coarse grid barely has five positive finite eigenvalues, so its fifth need not be an accurate continuum approximation. The CSV separates the selected eigenvalue and matrix residual from nodal and interpolated field values. A small residual is possible even on a poorly resolved grid.
The independent checks rebuild the full generalized equations, including the central constraint, and compare the refined eigenvalue with an Airy boundary determinant. The normalization at the centre fixes amplitude, not the number of oscillations or the eigenvalue. Change the resolution and predict which part of the shape should settle first.
The complete C++ model and starting input are available with this lesson. Keep collocation.hpp in the theme directory. This model also uses eigensystem.hpp in that directory.