Linear systems and numerical eigenvalue problems

State the linear problem and its units

A linear system relates unknown quantities through equations in which each unknown is multiplied by a supplied coefficient and the results are added. Collect two unknown displacement components into a column u = (u₁, u₂), expressed in metres in a fixed orthonormal two-dimensional basis. Orthonormal means the basis directions are perpendicular and each basis vector has unit length. Here the normalized basis directions are dimensionless; their component coefficients carry metres.

Write A = [[2, 1], [1, 3]] and b = (1, 2) m. The brackets list rows. The matrix equation Au = b means 2u₁ + u₂ = 1 m and u₁ + 3u₂ = 2 m. A is dimensionless, so each equation combines displacement quantities with the same unit. The supplied b is a column of linear combinations of displacement components; it is not a force column.

All matrices and quantities below are exact original teaching stipulations. They are not measured apparatus data or an established law of motion. In another problem, coefficients can carry units or components can describe different quantity kinds. Such a problem needs its own compatible units and scaling before a numerical norm is meaningful.

Eliminate with valid pivots

Elimination cancels an unknown from one equation using another equation. Apply the same operation to its right-hand side. For our system, subtract half the first equation from the second: 2.5u₂ = 1.5 m. Solve for u₂, then substitute it into the first equation to find u₁. These reversible equation operations preserve the solution set.

A pivot is a coefficient used for elimination or back substitution. Division by a zero pivot is undefined. A zero first pivot can sometimes be repaired by interchanging equations; it does not by itself show that the matrix is singular. Row interchange moves the entire equation, including its right-hand side.

For a two-by-two matrix [[a, b], [c, d]], the determinant is ad − bc. A nonzero determinant gives an inverse and a unique solution for each supplied right-hand side. Our A has determinant 5. Its inverse is A⁻¹ = [[3/5, −1/5], [−1/5, 2/5]]; multiplying A by A⁻¹ gives the identity matrix I = [[1, 0], [0, 1]].

A singular square matrix has no inverse. Its equations may conflict, or they may leave free choices. A repeated equation with a matching right-hand side does not conflict; a repeated left-hand side with a different right-hand side does. Inspect the transformed equations. Exact invertibility alone does not establish that a numerical implementation is stable.

Compare residual with solution error

Let u be the exact solution and û a candidate. Define the signed solution error e = û − u and the equation residual r = b − Aû. These are different displacement columns. Since Au = b, r = −Ae and e = −A⁻¹r when A is invertible. A residual checks how closely the candidate satisfies the equations; it is not itself the solution error.

For a displacement column z, define its infinity norm by ‖z‖∞ = max(|z₁|, |z₂|). It retains the component unit. For a dimensionless matrix M, define ‖M‖∞ as the largest sum of absolute entries in a row. The triangle inequality in each row gives ‖Mz‖∞ ≤ ‖M‖∞‖z‖∞. The units agree on both sides.

Combining that inequality with e = −A⁻¹r gives the exact bound ‖e‖∞ ≤ ‖A⁻¹‖∞‖r‖∞. For our A, ‖A⁻¹‖∞ = 4/5. A candidate residual of norm 0.02 m therefore bounds its solution error by 0.016 m. The bound need not be attained.

Aû = b − r also says that the candidate solves a system with the same A and an altered right-hand side. This is a narrow backward interpretation of this candidate. It does not prove that every output of an algorithm has small backward error. All bounds here use exact residual arithmetic; a computer's rounded residual requires its own error assessment.

Measure conditioning with a declared norm

Conditioning concerns how an exact problem responds to changes in its supplied data. It depends on the problem and the way size and relative change are measured. For an invertible dimensionless A, define κ∞(A) = ‖A‖∞‖A⁻¹‖∞. This dimensionless matrix condition number is useful for bounding sensitivity. It is not a physical uncertainty or a universal pass/fail score.

Keep A fixed and assume b ≠ 0. Then its solution u is also nonzero. From ‖b‖∞ = ‖Au‖∞ ≤ ‖A‖∞‖u‖∞ and the absolute error bound, obtain ‖e‖∞/‖u‖∞ ≤ κ∞(A)‖r‖∞/‖b‖∞. This is an exact inequality for the specified right-hand-side interpretation; changes in A require a separate analysis.

Our A has ‖A‖∞ = 4 and κ∞(A) = 16/5 = 3.2. A relative residual of 0.01 therefore gives a relative solution-error bound of 0.032, or 3.2%. A larger condition number can allow much greater amplification. It does not require every possible perturbation to attain that amplification.

The original sensitivity example below changes one right-hand-side value by 0.01 m while keeping its matrix exact. The solution changes by 1 m in the infinity norm. This happens in exact real arithmetic. Better arithmetic cannot remove the problem's input sensitivity; assessing an algorithm's added errors is a separate task.

Specify an eigenvalue problem

An eigenpair satisfies Cv = λv with v ≠ 0. The nonzero vector v is an eigenvector and the scalar λ is its eigenvalue. Multiplying v by any nonzero scalar gives the same eigenvalue and the same line of directions. Normalization selects a representative rather than a physical amplitude. The zero vector is excluded because it would satisfy the equation for every λ and select no eigenvalue.

Use an original dimensionless matrix C = [[3, 1], [1, 3]] in our fixed orthonormal basis. Take dimensionless direction columns v₊ = (1, 1) and v₋ = (1, −1). Direct multiplication gives Cv₊ = 4v₊ and Cv₋ = 2v₋. The eigenvalues are dimensionless here because C is dimensionless. This is not a rule that all physical eigenvalues have no units.

These two directions form a basis: every real two-component column can be written as a sum of multiples of v₊ and v₋. C is symmetric because its entries match across the diagonal. Real symmetric matrices have an orthonormal eigenvector basis, a property used later to interpret a residual. No motion, physical mode or stability law follows from the matrix alone.

Iterate with a stated starting direction

A power iteration repeatedly multiplies by a matrix and rescales. For our C, start x₀ = (1, 0). At each step calculate y = Cxₖ and set xₖ₊₁ = y/‖y‖∞, provided y ≠ 0. The positive norm rescales the column so its largest absolute component is 1. We use this norm normalization, rather than dividing by a signed component.

Here x₀ = (v₊ + v₋)/2. After k unscaled multiplications, Cᵏx₀ = (4ᵏv₊ + 2ᵏv₋)/2. The relative coefficient of v₋ is (2/4)ᵏ. After our normalization, xₖ = (1, (2ᵏ − 1)/(2ᵏ + 1)) for integer k ≥ 0. Thus these exact iterates approach (1, 1). The ratio 1/2 describes the relative eigencomponent decay; it is not an exact halving of every possible error measure per step.

A standard sufficient setting is a real matrix with a basis of real eigenvectors, a single eigenvalue strictly largest in absolute value, and a starting column with a nonzero component along that eigenvector. Exact normalized iteration approaches its eigenvector line. If that dominant eigenvalue is negative, positive norm normalization can alternate the representative's sign. Our example has a positive dominant eigenvalue. These conditions are narrower than claiming success from every nonzero start.

Starting our C from (1, −1) stays on that other eigenvector line in exact arithmetic. Finite iterations are generally approximations; a small change between two iterates is not by itself an eigenpair accuracy bound. Actual software must also handle zero normalization, finite range, rounding and a declared stopping rule.

Check an eigenpair residual with symmetry

For real two-component columns, the dot product is xᵀz = x₁z₁ + x₂z₂ and the Euclidean norm is ‖x‖₂ = √(x₁² + x₂²). For x ≠ 0, define the Rayleigh quotient ρ = (xᵀCx)/(xᵀx). It is a candidate eigenvalue, not an automatic exact one. Define the eigenpair residual s = Cx − ρx. All these quantities are dimensionless in our stated direction example.

The normalized residual ‖s‖₂/‖x‖₂ is unchanged when x is multiplied by a nonzero real scalar. For a real symmetric matrix, it bounds the distance from ρ to at least one exact eigenvalue. To see why, express x in an orthonormal eigenvector basis: ‖s‖₂²/‖x‖₂² is a weighted average of (λⱼ − ρ)² with nonnegative weights summing to 1. That average cannot be smaller than the smallest squared distance. The orthonormal-basis property is essential to this argument. The same argument works for any real candidate α with residual Cx − αx.

For our C and x = (1, 3/5), the quotient is ρ = 66/17. The residual is s = (−24/85, 8/17), and ‖s‖₂/‖x‖₂ = 8/17. The exact eigenvalue 4 lies 2/17 from ρ, within that bound. The nonzero residual also confirms that this finite iterate is not an exact eigenvector.

The bound concerns at least one eigenvalue; it does not on its own identify a particular eigenvector or bound its direction error. It cannot be transferred unchanged to every nonsymmetric matrix. The final exercise gives a nonsymmetric counterexample. A physical interpretation requires an independently specified physical model as well.

Record the numerical and physical limits

Record the equations or matrix, basis, quantity kinds, component units, supplied data, precision, method and pivot or normalization rule. State the residual sign convention, norm and stopping criterion. When a bound is used, retain its assumptions and distinguish an upper bound from the actual error.

Our worked values are exact reference calculations. A rounded matrix-vector product or solve can add errors, including errors in the reported residual itself. Small residuals, visually stable digits and successful exact hand examples do not establish a general stability theorem or validate an implementation. A condition estimate also depends on the norm and data scaling.

Changing a basis, using different component units or rescaling equations can change numerical diagnostics. Document such choices. Solving the supplied equations accurately is separate from showing that their coefficients and assumptions describe an experiment. No measurement uncertainty, probability distribution or real-world mode is assigned to these teaching values.

The next map lesson, Interpolation, extrapolation, and approximation, asks how a chosen function relates supplied values and what can be justified between or beyond them. Here the task is to solve and check a stated linear problem before giving its numerical output a physical interpretation.

Worked examples

Try it yourself

1. Solve Bq = c for B = [[0, 1], [1, 1]] and c = (2, 3) m. Does the zero top-left entry make the system singular?

   Show the answer and reasoning

   No. The solution is q = (1, 2) m and the determinant is −1.

   The first equation directly gives q₂ = 2 m; the second then gives q₁ = 1 m. Alternatively, interchange the two entire equations to obtain a nonzero first pivot. A zero entry at the chosen pivot position is not the same as a zero determinant.

2. For S = [[1, 2], [2, 4]], compare the right-hand sides (3, 6) m and (3, 7) m. Is there a unique solution in either case?

   Show the answer and reasoning

   Neither is unique. The first has infinitely many solutions; the second has none.

   The determinant is 0. Subtracting twice the first equation from the second gives 0 = 0 m for (3, 6) m, leaving q₁ + 2q₂ = 3 m with a free real choice of q₂. For (3, 7) m it gives the contradiction 0 = 1 m. Singularity permits these different outcomes.

3. For the original A and b, assess û = (0.19, 0.61) m. Find r, e and the absolute infinity-norm bound.

   Show the answer and reasoning

   r = (0.01, −0.02) m, e = (−0.01, 0.01) m, and ‖e‖∞ ≤ 0.016 m.

   Aû = (0.99, 2.02) m. Subtract it from b for r; subtract the exact (0.2, 0.6) m from û for e. Their norms are 0.02 m and 0.01 m. The inverse norm remains 4/5, giving the bound (4/5)(0.02 m) = 0.016 m. Signed columns differ even when their norm values agree with the worked example.

4. Multiply both A and b in the first worked system by the dimensionless number 100, retaining its same candidate û. What happens to the exact solution, residual, relative residual and matrix condition number?

   Show the answer and reasoning

   The solution is unchanged; the residual is multiplied by 100; the relative residual and κ∞ stay unchanged.

   The equation 100Au = 100b is equivalent to Au = b. Its residual is 100r and its right-hand-side norm is 100‖b‖∞, so their ratio is unchanged. Meanwhile ‖100A‖∞ = 100‖A‖∞ and (100A)⁻¹ = A⁻¹/100, making the norm product invariant. An absolute residual depends on equation scaling. This is an algebraic rewriting, not a new physical measurement.

5. For the original C, start normalized power iteration from x₀ = (1, −1). Will the exact iteration approach the eigenvector line for eigenvalue 4?

   Show the answer and reasoning

   No. It remains xₖ = (1, −1), the eigenvector representative for eigenvalue 2.

   C(1, −1) = (2, −2). Dividing by its positive infinity norm 2 returns the starting column at every step. This column has no component along (1, 1), so the sufficient dominant-start condition fails. Do not rely on unspecified rounding to repair that exact starting choice.

6. Take the nonsymmetric dimensionless N = [[1, 100], [0, 1]], candidate eigenvalue α = 2 and v = (1, 0.01). Calculate Nv − αv and compare its normalized Euclidean norm with the distance from α to the exact eigenvalues.

   Show the answer and reasoning

   The residual is (0, −0.01), its normalized norm is less than 0.01, and the only exact eigenvalue is 1, at distance 1 from α.

   Nv = (2, 0.01), while αv = (2, 0.02). Thus the residual norm is 0.01 and ‖v‖₂ = √1.0001 > 1. The characteristic determinant is det(N − λI) = (1 − λ)², so both eigenvalues counted with multiplicity equal 1. This original counterexample shows why the symmetric nearest-eigenvalue bound cannot be asserted for every real matrix.

What to carry forward

• Specify quantity kinds, units and a basis before interpreting a matrix equation numerically.
• A zero chosen pivot may be repairable by row interchange; a singular system may be inconsistent or leave free choices.
• Equation residual and solution error are different, with a fixed invertible matrix relating them.
• A condition number and a declared norm help bound sensitivity; a bound is not necessarily the attained error or measurement uncertainty.
• Power iteration needs stated spectral and starting-direction conditions, and its finite iterates need residual checks.
• A real symmetric eigenpair residual can bound proximity to an eigenvalue; retain symmetry and distinguish that statement from direction accuracy and physical validity.