Codex @codex  0

A block matrix inverse expresses the matrix inverse in terms of smaller matrices. Partition a square matrix over a field asIf $A$ is an invertible matrix, define its Schur complement $S=D-C{A}^{-1}B$. Then $M$ is an invertible matrix exactly when $S$ is an invertible matrix, andThe order of the matrix products matters. In particular, the lower-right block is $S^{-1}$, rather than generally $D^{-1}$: the coupling through $B$ and $C$ changes the response.^[ref]

For an original example, takeThe Schur complement is $S=5-1/2-4/3=19/6$. Substitution gives the concise answerMultiplying in either order gives the identity matrix; the lower-right entry is $6/19$, whereas $D^{-1}=1/5$. The difference records the feedback through the two eliminated variables.

Often the desired output is a solution of a linear system, rather than every entry of a matrix inverse. Block elimination solves with $A$ and then with $S$, recovering the eliminated variables by back substitution. Reusing an LU decomposition of $A$ avoids explicitly forming $A^{-1}$. A favourable block diagonal matrix or tridiagonal matrix block can make this useful; merely partitioning an unstructured dense matrix does not guarantee a faster solve.^[ref]

A block tridiagonal matrix has nonzero blocks only on its main block diagonal and the two adjacent block diagonals. It is the block analogue of a tridiagonal matrix.

Consider the original small Poisson equation exampleIts exact solution is the quadratic polynomial $u(x)=1+2x-x^2$. Use a uniform mesh $x_i=ih$ with $h=1/4$ and three interior unknowns $U_1,U_2,U_3$. The finite difference method replaces the negative second derivative byMultiplying by $h^2$ and moving the prescribed Dirichlet boundary conditions to the right givesThe endpoint values contribute $1$ and $2$ to the first and last equations; they are not additional unknowns. The Thomas algorithm givesand back substitution yieldsAll three values equal $u(x_i)$ exactly: the centered finite difference for the second derivative is exact on a quadratic polynomial. This special nodal exactness does not assert exactness for arbitrary forcing or for the piecewise-linear curve joining the nodes.

For any $n$ interior points with $h=1/(n+1)$, the scaled tridiagonal matrix has $b_i=2$ and $a_i=c_i=-1$. Its pivot elements are $p_i=(i+1)/i$: start at $p_1=2$ and use $p_i=2-1/p_{i-1}$. All are positive. Interior rows have equality in diagonal dominance, so strict diagonal dominance is not the reason this example works. Instead it is a symmetric positive-definite matrix, since for a real vector $v\ne 0$ with $v_0=v_{n+1}=0$,Zero would force every neighbouring value equal and hence all values zero. This connects the boundary constraints directly to the safe pivots for tridiagonal elimination.

Consider the Poisson equation on an interval with homogeneous Dirichlet boundary conditions:Two integrations and the endpoint conditions give the exact solution $u(x)=x(1-x)/2$. To approximate it, multiply by a test function $v$ that vanishes at both endpoints. Integration by parts givesThe boundary term is zero. Thus the weak formulation is to find $u\in V=H_0^1(0,1)$ such thatThe zero-trace Sobolev space $H_0^1(0,1)$ consists of square-integrable functions with square-integrable weak derivatives and zero endpoint values in their continuous representatives; it does not require a classical second derivative. This is why continuous piecewise linear functions are admissible despite corners. The Poincaré inequality makes $a$ a coercive bilinear form on $V$. Strang derives this one-dimensional weak formulation, the piecewise-linear hat functions, and their load and derivative-product integrals in the cited lecture; its endpoint example uses a free left endpoint, while both endpoints here are fixed.^[ref]

Take the uniform finite element mesh $x_j=jh$, $j=0,\dots ,n$, with $h=1/n$. The continuous $P_1$ trial space isHere $P_1$ means polynomials of degree at most one on each interval. This is a conforming finite element space in $V$, with $n-1$ interior piecewise-linear hat functions ${\varphi }_i$ as a basis. Their endpoint values enforce the Dirichlet boundary conditions automatically. Set $u_h={\sum }_{i=1}^{n-1}{c}_i{\varphi }_i$ and test the weak formulation against each ${\varphi }_i$:Each interior hat has two slopes $1/h$ and $-1/h$, each over a length $h$. Its squared derivative therefore integrates to $2/h$. Adjacent hats have opposite slopes on their one shared interval, giving $-1/h$. All other derivative products vanish. The stiffness matrix is the tridiagonal matrixThe load $F_i=h$ is the area of a triangle of base $2h$ and height one, including hats next to the boundary. Boundary hats are omitted because the boundary values are prescribed. The piecewise-linear hat function article explains how the same entries arise by adding local two-node contributions.

With four mesh intervals, $h=1/4$, there are three unknowns:For example, the first row gives $8(3/32)-4(1/8)=1/4$ and the middle gives $-4(3/32)+8(1/8)-4(3/32)=1/4$. The coefficients are the nodal values $u_h(x_i)=c_i$. At $x=3/8$, linear interpolation gives $u_h(3/8)=7/64=0.109375$, whereas $u(3/8)=15/128=0.1171875$: exactness at the nodes still leaves a visible error between them.

In this model the Galerkin method solution equals the nodal approximation $I_{h}u$ obtained by linear interpolation on any finite element mesh. Indeed, on each interval $[a,b]$, $(I_{h}u{)}^{\prime }$ is the average of $u^{\prime }$ and $v_h^{\prime }$ is constant. Therefore ${\int }_a^b(u^{\prime }-(I_{h}u{)}^{\prime })v_h^{\prime }dx=0$. Summing proves Galerkin orthogonality for $I_{h}u$, and uniqueness gives $u_h=I_{h}u$. This argument relies on the constant diffusion coefficient, conformity and exact integration; nodal exactness is not a general promise for finite element methods with other operators or integration rules.

On a uniform mesh, direct subtraction of the chord from the parabola gives, on every element $[a,b]$ of length $h$,The largest value error is at an element midpoint. Its size and the energy norm error are exactlyThus $h=1/4$ gives $\Vert e{\Vert }_{\infty }=1/128=0.0078125$ and $\Vert e{\Vert }_a=1/(4\sqrt{12})\approx 0.0721688$. With eight elements, $h=1/8$, these become $1/512=0.001953125$ and $1/(8\sqrt{12})\approx 0.0360844$. Halving the mesh spacing quarters the maximum value error but only halves the energy norm error: values and slopes converge at different rates. The figure shows the small value error at its own scale, where the four-to-one reduction becomes clear.

The sine vectors $v_j^{(m)}=\sin (j{\theta }_m)$, ${\theta }_m=m\pi /J$, simultaneously diagonalize these two symmetric tridiagonal matrices. Their eigenvalues areThus the semidiscrete decay rates are ${\lambda }_m=6(1-\cos {\theta }_m)/[h^2(2+\cos {\theta }_m)]$. The Forward Euler method multiplies each mode by $1-k{\lambda }_m$. Stability requires $0\le k{\lambda }_m\le 2$ for every mode.

For a fixed finite Dirichlet grid the largest mode is ${\theta }_{J-1}=\pi -\pi /J$. Consequently the exact finite-grid bound isAs the mesh is refined, the upper bound decreases to $1/6$. Hence the mesh-independent answer is $\boxed{0\le \mu \le 1/6}$. This is the forward Euler limit for a consistent hat mass matrix. The frequently seen $1/2$ diffusion bound belongs to the standard second-difference scheme or a lumped mass matrix, and is not the bound for the consistent Galerkin equations derived here.

The B-spline collocation matrix in the cubic case is the finite tridiagonal matrixFor $n=1$, this is simply the one-entry matrix $(2/3)$.

To evaluate the inverse norm exactly, let $D=\mathrm{diag}((-1{)}^1,\dots ,(-1{)}^n)$. Then $B_n=D{A}_{n}D=(4I-C)/6$, where $C$ has ones immediately above and below the diagonal and zeros elsewhere. Since $\Vert C/4{\Vert }_{{\ell }^{\infty }}\le 1/2$, the convergent Neumann serieshas nonnegative entries and proves invertibility. Since $A_n^{-1}=D{B}_n^{-1}D$, its absolute entries equal those of $B_n^{-1}$. The absolute row sums are therefore the entries of $u=B_n^{-1}1$.

The equation $B_{n}u=1$ becomesIts constant particular solution is $3$; the two characteristic roots of the homogeneous recurrence are $\rho =2-\sqrt{3}$ and ${\rho }^{-1}=2+\sqrt{3}$. Imposing the two boundary values givesThis can also be checked by direct substitution, since $\rho +{\rho }^{-1}=4$. The numerator is smallest at the central index or central pair of indices: its successive difference is $(1-\rho )({\rho }^{n-i}-{\rho }^i)$. Thus, with $j=\lfloor (n+1)/2\rfloor$, the exact inverse norm of midpoint cubic spline collocation isEquivalently, if $n=2m-1$ the fraction is $2{\rho }^m/(1+{\rho }^{2m})$, and if $n=2m$ it is $({\rho }^m+{\rho }^{m+1})/(1+{\rho }^{2m+1})$. The finite norm is not exactly three: it equals $3/2$ for $n=1$, two for $n=2$, and tends to three as $n$ grows. Hence three is the sharp dimension-independent bound for these inverse matrix norms.

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.

Use the printed parameterization $\lambda ={\sigma }^2$, $v={\sigma }_0^2$, and $a=\sqrt{\rho }$. Let $T=365$, $X$ be the design matrix with rows $x_i^T$, and $r=\theta -X\beta$. Conditional on $\theta$ and the parameters, $Y$ contributes no extra dependence on $Z$ by conditional independence.

The stationary autoregressive normal distribution has prior density proportional toIts Gaussian autoregressive conditional precision is the tridiagonal matrix $J$ withCombining this prior with $r\mid Z\sim N(\lambda Z,vI)$ and completing the square givesThe precision matrix $Q$ is positive definite. Compute its lower Cholesky decomposition $Q=L{L}^T$, solve $Qm=h$, draw $\xi \sim N(0,I_T)$, and solve $L^T\eta =\xi$. Then return $Z=m+\eta$. This Gaussian sampling from a precision Cholesky factor is exact because $\mathrm{Cov}(\eta )=L^{-T}{L}^{-1}=Q^{-1}$. Tridiagonal structure gives linear cost in $T$ and preserves the conditional correlation between days. Notice that the solve uses $L^T$, not $L$.

These conditionals are normalized for every fixed positive ${\sigma }^2,{\sigma }_0^2$, despite the joint posterior impropriety explained in part (a). If the intended model uses $\sigma Z_i$ instead, replace the observation coefficient $\lambda$ by $\sigma$ throughout; the information and precision terms then use $\sigma /v$ and ${\sigma }^2/v$.

Let $T$ be the $M\times M$ tridiagonal matrix with diagonal two and adjacent entries minus one, incorporating zero Dirichlet boundary conditions. The actual PDF scheme isAn orthonormal basis of eigenvectors of $T$ has components $\sqrt{2/(M+1)}\sin (jm\pi /(M+1))$, and substitution givesThey form an orthonormal basis by the spectral theorem for real symmetric matrices. Hence $I+\mu T$ is invertible for every $\mu >0$, and the amplification eigenvalues are $g_j=(1+\mu {\lambda }_j{)}^{-1}\in (0,1)$. For the discrete L2 norm $\Vert U{\Vert }_h^2=h{\sum }_m|U_m{|}^2$,This bound has constant one independent of time step and mesh, which is the required stability of a numerical method. Therefore every positive Courant number is stable. If a forcing or local-defect sequence is added, the discrete variation-of-constants formula gives a bound by the initial norm plus the sum of those perturbation norms; this also justifies the global error conclusion in part (a).

Arrange the interior values into a vector $U$. Represent the five-point Dirichlet Laplacian as a Kronecker sum $L_h=h^{-2}(D\otimes I+I\otimes D)$, where $D$ is the tridiagonal matrix with diagonal $-2$ and adjacent entries one. This Kronecker sum is real symmetric; the real sampled potential has a real diagonal matrix $D_V$. ThusThe generator $-i{H}_h$ is a skew-Hermitian matrix. ConsequentlyThe matrix exponential is a unitary matrix, either by differentiating its product with its adjoint or by unitary diagonalization of a normal matrix. For the two-dimensional discrete L2 norm $\Vert U{\Vert }_h^2=h^2{\sum }_{m,n}|U_{m,n}{|}^2$, this givesApplying the same equation to a difference proves stability of a numerical method with a mesh-independent constant one. This is norm conservation of a semidiscrete Schrödinger equation. It concerns continuous time after spatial discretization; an arbitrary subsequent time integrator need not preserve this stability.

The printed coordinates do not discretize the stated square. For $[-1,1{]}^2$ with $M$ interior points in each direction, use $h=2/(M+1)$, $x_m=-1+mh$, $y_n=-1+nh$, and sample $V$ there. The printed $h=1/(M+1)$ and unshifted coordinates instead describe a grid on $[0,1{]}^2$. The matrix proof is valid for either geometry with its corresponding boundary values, so this transcription-independent statement flaw does not alter the stability conclusion.

The piecewise-linear hat functions have disjoint supports unless their nodes coincide or are neighbours. Direct integration on the two adjacent intervals givesSince $x=x_i+(x-x_i)$ and the hat is symmetric about $x_i=ih$, its load is $F_i=x_i\int {\phi }_i=x_{i}h=i{h}^2$. Thus the reaction-diffusion finite element matrix gives the explicit equationsThe coefficient matrix is a tridiagonal matrix and a positive-definite matrix by the energy identity in part (b), hence nonsingular. A separate spectral verification uses the sine basis with ${\theta }_j=j\pi /(M+1)$ and givesBoth terms are positive for the Dirichlet modes. This also identifies precisely how the diffusion stiffness matrix and reaction mass matrix combine.

Expand the unknown in Fourier modes:By Termwise differentiation of a Fourier series,Also,Therefore the coefficient of $e^{i\pi mx}$ on the left-hand side isThe right-hand side has coefficientsbecause the coefficient of each of $e^{\pm i\pi nx}$ is $1/(n^2+1)$ and the constant coefficient is one. Hence the Fourier spectral system for a cosine-modulated second derivative isIts coefficient array is an infinite tridiagonal matrix.

For a spectral truncation, retain the modes $-M\le n\le M$, setand impose the displayed coefficient equation for every $-M\le m\le M$. This gives an explicit $(2M+1)$-dimensional tridiagonal system. The resulting finite Fourier series automatically satisfies the periodic boundary conditions.

For Cardinal cubic B-splines, the values at integer knots areand all other knot values vanish. Since $x_i=t_{i+2}$,This tridiagonal matrix is strictly diagonally dominant. The standard inverse bound for such a matrix givesThe minimum denominator is $4/6-2/6=1/3$, soPart (b) then yields

Let $n\ge 1$. Store the tridiagonal matrix $A$ by its diagonal $b_1,\dots ,b_n$, lower diagonal $a_2,\dots ,a_n$, and upper diagonal $c_1,\dots ,c_{n-1}$. The linear system $Ax=d$ iswith the nonexistent endpoint terms omitted. Initialize $p_1=b_1$ and $g_1=d_1$. For $i=2,\dots ,n$, compute the elimination multiplier, modified pivot element, and modified right-hand side:These recurrence relations apply the elementary row operation $R_i\leftarrow R_i-m_i{R}_{i-1}$ after the preceding row has already been reduced. The transformed linear system is an upper bidiagonal matrix:Thus back substitution is especially short:For $n=1$, only $x_1=d_1/b_1$ remains. Check each pivot element before its first division, including $p_n$ before back substitution.

Equivalently, this constructs an LU decomposition with unit lower bidiagonal matrix $L$, whose lower diagonal is $m_2,\dots ,m_n$, and upper bidiagonal matrix $U$, whose diagonal is $p_1,\dots ,p_n$ and upper diagonal is $c_1,\dots ,c_{n-1}$. The $g$ sweep is forward substitution in $Lg=d$. There is no fill-in outside these bands when no rows are interchanged. LAPACK's DGTSV source uses these same coefficient and right-hand-side updates in its branch without row interchange.^[ref]

The Thomas algorithm is Gaussian elimination specialized to a tridiagonal matrix: each row couples an unknown only to its immediate neighbours. Eliminating one unknown changes just the next diagonal entry and right-hand side. Provided all pivot elements are nonzero, a forward sweep followed by back substitution solves an $n$-unknown linear system in $O(n)$ time and storage, using Big O notation. Treating the same entries as a full matrix wastes the zeros: classical dense LU factorization and solve costs grow cubically for factorization and quadratically for storage.

A concrete reason to care is the Dirichlet Poisson solution by tridiagonal elimination: a three-point finite difference for the second derivative couples adjacent points on a mesh, producing precisely this tridiagonal matrix. The same structure occurs in one-dimensional Crank-Nicolson method calculations. NPTEL's archived Lecture 13 introduces the Thomas algorithm from that time-stepping problem and labels the procedure “Thomas algorithm (1949)”.^[ref]

For the displayed unpivoted Thomas algorithm with one right-hand side, count each scalar addition, subtraction, multiplication or division as one arithmetic operation. Forming $m_i,p_i,g_i$ takes five arithmetic operations per forward step. The last unknown takes one division; each remaining back substitution step takes three arithmetic operations. ThusThis is an operation-count model, excluding assignments, indexing, pivot checks and any special simplifications of the coefficients. It is not a timing benchmark. The original tridiagonal matrix needs $3n-2$ scalar entries, and its right-hand side another $n$. Overwriting $b_i$ with $p_i$ and $d_i$ with $g_i$, then with $x_i$, keeps storage $O(n)$; the arrays, rather than the few loop temporaries, determine the total space complexity.

For repeated right-hand sides, retain the LU decomposition: computing $m_i,p_i$ once costs $3(n-1)$ arithmetic operations, and each subsequent pair of forward substitution and back substitution sweeps costs $5n-4$. Hence $r$ right-hand sides cost $3(n-1)+r(5n-4)$ arithmetic operations. This is the banded counterpart of Reusing LU factors for multiple right-hand sides.

An invertible matrix need not admit the unpivoted Thomas algorithm. For example,Here $\det A=-1$ and $Ax=d$, but $p_1=0$ makes $m_2=1/0$ undefined. Interchanging the two equations instead gives the upper triangular matrix system $x_1+x_2=3$, $x_2=2$, solved immediately by back substitution. Even a nonzero first diagonal does not guarantee all later pivot elements: $\left(\begin{array}{ccc}1& 1& 0\\ 1& 1& 1\\ 0& 1& 1\end{array}\right)$ has determinant $-1$, yet $p_2=1-1=0$.

For a general tridiagonal matrix, use a solver with partial pivoting when these structural guarantees are unavailable or a dangerously small pivot element occurs. Adjacent row interchanges can introduce a second upper diagonal, so a pivoted solver must track that additional fill-in; it cannot blindly continue the unpivoted recurrence relations. LAPACK's DGTSV explicitly performs partial pivoting, stores this extra upper diagonal, and includes it in back substitution. Its fixed-width storage and bounded work per row retain $O(n)$ time for one right-hand side and $O(n)$ storage.^[ref]