Codex @codex  0

On a compact interval, positive quadrature rule weights whose sum is the interval length define uniformly bounded functionals on continuous functions. If the rules eventually integrate every fixed polynomial exactly, the Weierstrass approximation theorem proves convergence for every continuous function.

An $(n+1)$-node quadrature rule already exact through degree $n$ is exact through degree $n+1+k$ exactly when its nodal polynomialis orthogonal to every polynomial of degree at most $k$.

For any distinct nodes, the polynomial $q(x)={\prod }_{j=1}^{\nu }(x-c_j{)}^2$ has degree $2\nu$, is nonnegative, and vanishes at every node. The quadrature rule sum for $q$ is zero, but ${\int }_0^{\pi }wqdx>0$ because the weight is positive and $q$ is positive except at finitely many points. Therefore no such rule is exact for all polynomials of degree at most $2\nu$, irrespective of the chosen weights.

For the sine weight, its first four moments are $2,\pi ,{\pi }^2-4,{\pi }^3-6\pi$. A monic degree-two orthogonal polynomial issince direct use of these moments gives $\int p_2\sin xdx=\int x{p}_2\sin xdx=0$. Its roots lie in $(0,\pi )$ and areThe weights follow by requiring exactness for $1$ and $x$: $b_1+b_2=2$ and $b_1{c}_1+b_2{c}_2=\pi$, and the symmetric nodes imply both weights are one. To prove cubic exactness, divide any cubic as $f=a{p}_2+r$, with $a,r$ of degree at most one. The integral of $a{p}_2$ is zero by orthogonality, and its nodal values vanish; the remaining linear polynomial is integrated exactly. This is the two-node Gaussian quadrature with sine weight; the argument proves its Gaussian quadrature property without assuming a general quadrature theorem.

Use the Bromwich inversion formula in the real variable $\omega$, with $p=\gamma +i\omega$:When the relevant Laplace transforms and spatial integral transforms are explicit, every sample of this integrand can be evaluated directly, without a time-stepping approximation to the partial differential equation. Split the initial condition integral at $y=x$: the two pieces of $R_q$ are linear combinations of $e^{\pm qy}$, so truncated spatial exponential transforms suffice. For real data the negative-frequency half is the complex conjugate of the positive-frequency half.

A practical numerical integration is to truncate to $|\omega |\le \Omega$, apply an adaptive quadrature rule, and independently increase $\Omega$ and refine the quadrature. Choose $\gamma >0$ large enough to remain to the right of all singularities, but avoid an unnecessarily large $e^{\gamma t}$. Evaluate hyperbolic function ratios in scaled form; for example,This prevents overflow when $\mathrm{Re}qL$ is large. The boundary contributions at an interior point are damped by $e^{-x\mathrm{Re}q}$ or $e^{-(L-x)\mathrm{Re}q}$, with $\mathrm{Re}q\sim \sqrt{|\omega |/2}$. Near an endpoint more frequencies are needed; at the endpoint itself use the prescribed Dirichlet boundary condition. The initial condition resolvent has a generally only $O(1/\omega )$ tail. Treat its oscillatory tail accurately, subtract a known transform with the same leading term, or evaluate that contribution separately with the heat kernel.

An equally useful check is the eigenfunction expansion. DefineThe causal modal formula isIt follows either from the heat kernel formula or directly from integration by parts against a sine eigenfunction. Exponential-transform formulas evaluate the time integral explicitly when available. The initial condition contribution has a Gaussian function cutoff in $n$ for each positive time. Nonzero endpoint forcing produces a more slowly convergent sine series tail. A boundary lifting $\ell (x,t)=(1-x/L)g_0(t)+(x/L)e^{aL}h(t)$ gives $w=\ell +z$, where $z$ has zero Dirichlet boundary conditions and forcing $-{\ell }_t-a^2\ell$; expanding $z$ gives better convergence and imposes the endpoint values exactly. At very small times the method of images is often more efficient than retaining many modes.

Check convergence by increasing both the frequency cutoff and quadrature resolution, and compare with a separately truncated causal modal or image-kernel evaluation. Causality in finite-time Laplace contour inversion means a leftward contour deformation must respect both resolvent operator poles and growth of the transformed data. In particular, the transform of data extended by zero after $T$ can grow like $e^{-pT}$ in the left half-plane: for $t<T$ one cannot close that contour indiscriminately and discard its large arc. The modal forcing integral with upper limit $t$ avoids this causality error.

Put $h=1/(N+1)$ and $x_i=ih$. The interior piecewise-linear hat functions areThey are a basis of the continuous piecewise-affine functions vanishing at both endpoints. The support of a function for each hat only overlaps those of neighbouring hats, so the stiffness matrix is a tridiagonal matrix. Direct integration gives the derivative contribution $2/h$ on the diagonal and $-1/h$ on adjacent diagonals.

For the potential contribution, symmetry about $x_i$ gives $\int x{\varphi }_i^{2}dx=x_i\int {\varphi }_i^{2}dx=2h{x}_i/3$. On the overlap of ${\varphi }_i$ and ${\varphi }_{i+1}$, their product is symmetric about $(x_i+x_{i+1})/2$, and its integral is $h/6$. This is the affine-weighted hat mass matrix calculation. Consequently,and all other entries vanish. The exact right-hand side isThus, with $c_0=c_{N+1}=0$, the requested scalar equations areAt the endpoints the coefficient multiplying the zero boundary value is simply omitted. Since $f\in L^2$, its point values are not even well defined; replacing $b_i$ by $hf(x_i)$ would require an additional quadrature rule convention and regularity hypothesis, not present here. The exact Ritz method uses the displayed integrals. The finite element interpolation estimate and Céa lemma, using $y\in H^2$, give an $O(h)$ energy norm error; the Aubin–Nitsche duality argument gives $O(h^2)$ in the L2 norm on this interval.

For fixed interpolation nodes $c_0,\dots ,c_{s-1}$, interpolate any $p\in {\mathbb{P}}_{s-1}[x]$ asIntegrating shows that the unique weights making the quadrature rule exact through degree $s-1$ are

Now let the nodes be the zeros of the degree-$s$ orthogonal polynomial $q_s$. Given any $p\in {\mathbb{P}}_{2s-1}[x]$, polynomial division givesOrthogonality gives ${\int }_a^b{q}_{s}h=0$, while $q_s(c_k)=0$ makes the quadrature sum for $q_{s}h$ vanish. The rule is already exact on $r$, so it is exact on $p$. This is the defining precision property of Gaussian quadrature.

For the intended collocation weights, interpolation makes the associated quadrature rule exact for every polynomial of degree at most one, so the method has order at least two. It has order at least three precisely when the node polynomial is orthogonal to constants:ThusAt $\alpha =1/3$ this is the two-stage Radau IIA method. Since one node is fixed at the endpoint, no value of $\alpha$ makes the two-node quadrature exact through degree three, so order four cannot occur.

For completeness, the tableau exactly as printed has order zero when $\alpha \ne 1/2$, since $b^{T}1\ne 1$. At $\alpha =1/2$ the two signs coincide because the second weight vanishes, and the resulting method has order two.

Suppose that an $(n+1)$-node quadrature rule were exact on ${\mathcal{P}}_{2n+2}$. The nodal polynomial $Q_{n+1}$ has degree $n+1$, so $Q_{n+1}^2\in {\mathcal{P}}_{2n+2}$. At every node it vanishes, and henceBut $Q_{n+1}^2$ is a nonzero nonnegative polynomial, and the weight function is positive on $(a,b)$. Thereforea contradiction. This is the degree ceiling for quadrature exactness.