Alexander Godunov Personal website

Companion Website

Computational Physics

From Equations to Simulation

Chapter 8

Linear Systems and Eigenvalue Problems

These exercises supplement the book problems with checks on pivoting, sensitivity, iteration, and eigenvector interpretation. The projects use heat flow, electrostatics, and diffusion instead of the book’s spring, circuit, and quantum examples. Start with the stated inputs, report residuals and the meaning of each result, and then vary the parameters to see which conclusions persist.

Chapter 8 codes Chapter 8 resources Chapter 8 updates & errata Chapter 8 overview

Additional exercises

  1. A small pivot in an otherwise simple system

    Solve A*x = b with A = [[10−14, 1], [1, 1]] and b = [1 + 10−14, 2] in double precision. Compare elimination without pivoting, the companion Gaussian elimination code, and a library solve. Compute max|A*x − b| and max|x − [1,1]| using the original data. Why can the order of elimination matter even though the exact solution is [1,1]? Repeat with the top-left entry set to 10−8 and adjust b accordingly.

  2. A small residual and sensitive data

    Set δ = 10−8, A = [[1,1], [1,1+δ]], and b = [0,−δ]. Solve the system and check against x = [1,−1]. Then replace only the second component of b by −δ + 10−10 and solve again. The exact new solution is [0.99,−0.99]. Compare each solve’s residual against its own input data and the change in the solution. Repeat with δ = 10−6. Explain why a small residual alone does not measure uncertainty in the supplied physical data.

  3. Relaxation is a parameter, not a guarantee of speed

    Apply SOR to A = [[5,1,0], [1,5,1], [0,1,5]] and b = [7,10,−3]. Start from x = [0,0,0], use ω = 0.8, 1.0, 1.3, 1.8, and stop at max|A*x − b| ≤ 10−10 or 500 sweeps. Record the sweep counts and compare each result with [1,2,−1]. Which tested value is fastest for this particular system? Vary the off-diagonal entries while keeping the matrix strictly diagonally dominant.

  4. When the Thomas algorithm needs a pivot

    Consider the tridiagonal matrix A = [[0,1,0], [1,2,1], [0,1,2]] with b = [1,4,3]. Try the companion Thomas solver and record its failure indication. Then solve with a pivoted library routine and verify x = [1,1,1] to a residual below 10−12. Explain why tridiagonal structure alone does not make the no-pivot algorithm safe. Change the first diagonal entry to 2, recompute b = A*[1,1,1], and try again.

  5. What a converged power method has found

    Use A = diag(−5,2,1) and a residual tolerance ||A*y − λy||2 ≤ 10−8, with at most 500 iterations. First start from y = [1,1,1] and check that the eigenvalue with largest magnitude is negative. Then start from y = [0,1,1]. The latter run can converge with a small residual to λ = 2, because its initial vector has no component along the −5 eigenvector. Explain why a small residual ||A*y − λy||2 validates a computed eigenvalue and eigenvector but does not establish that the eigenvalue is dominant.

  6. Choosing a shift for one eigenvalue

    For A = [[2,1,0], [1,2,0], [0,0,6]], compute the full spectrum with Jacobi or a symmetric library solver. Then use shifted inverse iteration with y = [1,0,1], shift σ = 2.6, residual tolerance 10−10 for ||A*y − λy||2, and a limit of 100 iterations. It should select λ = 3 from the exact spectrum 1, 3, 6. Try σ = 3 and explain the singular shifted solve. At σ = 2, the eigenvalues 1 and 3 are equally close. With this starting vector their components have equal magnitude, so the iteration stalls on a mixture and reaches the 100-iteration limit without meeting the residual tolerance. Explain why a unique convergence target is no longer assured.

  7. Repeated eigenvalues and nonunique eigenvectors

    Find all eigenvalues and eigenvectors of A = [[3,1,0], [1,3,0], [0,0,2]] with the companion Jacobi program and a symmetric library solver. Use a squared off-diagonal tolerance of 10−18 and at most 100 Jacobi sweeps. Verify ||A*v − λv||2 < 10−8 and ||XTX − I||∞ < 10−8. The eigenvalues are 2, 2, 4. The two vectors belonging to 2 may differ between solvers: compare their two-dimensional span instead of matching their components one by one.

Physics projects

Use a library solve as an independent check when a project grows beyond the small examples in the companion codes. State the boundary conditions and units alongside every plotted result.

1. Temperature in a heated rod

Objective: compute the steady temperature of a rod with uniform internal heating and fixed temperatures at both ends.

Let −k T″(x) = q for 0 ≤ x ≤ L, with T(0) = T(L) = 300 K. Use L = 1 m, k = 200 W/(m K), and q = 10,000 W/m3. Divide the rod into 20 equal intervals, giving 19 interior unknowns and h = 0.05 m. The centered difference produces a tridiagonal system with diagonal 2, adjacent entries −1, and right-hand side q h2/k; add the boundary temperature to the first and last entries of the right-hand side. Solve with Thomas and with a library routine. Plot temperature against position and inspect its midpoint and symmetry.

Validation: the exact solution is T(x) = 300 K + (q/2k)x(L − x). It gives T(0.5 m) = 306.25 K. Centered differences reproduce this quadratic at every grid point to roundoff. Check the original discrete equation and compare each interior temperature with the exact profile.

Possible extension: set the right boundary to 320 K. Derive the new analytical solution and determine where the hottest point lies.

2. Screened potential in a grounded square

Objective: solve a two-dimensional boundary-value problem as a sparse linear system and see how screening changes its potential.

In a square 0 ≤ x,y ≤ 1 m, set V = 0 on all four boundaries and solve −∇2V + κ2V = s with κ = 2 m−1 and s = 1 V/m2. Use nine interior grid points along each axis, so h = 0.1 m and the matrix has 81 unknowns. The five-point stencil has diagonal 4 + (κh)2, entries −1 for each existing interior neighbor, and right-hand side h2s. Start SOR from zero, try ω = 1 and ω = 1.5, and require a maximum absolute discrete residual below 10−10 V within 5000 sweeps. Compare a contour plot with a sparse library solve.

Validation: the potential is symmetric about x = y = 0.5 m and is largest at the center. An independent direct solve of the stated 81-by-81 system gives V(0.5 m,0.5 m) ≈ 0.059673191 V. For the stated traversal order, zero initial guess, and residual threshold, ordinary Gauss–Seidel (ω = 1) takes about 159 sweeps; ω = 1.5 takes about 35. Verify the residual independently after either run.

Possible extension: vary κ and track the center potential. Stronger screening lowers it. At κ = 0, the discrete matrix loses strict diagonal dominance but remains symmetric positive definite. Gauss–Seidel and SOR with 0 < ω < 2 still converge for this finite grid. Compare their sweep counts with the screened case.

3. Diffusive relaxation modes

Objective: use eigenvalues to predict how nonuniform concentration decays in an insulated one-dimensional region.

For ∂u/∂t = D ∂2u/∂x2, use L = 1 m, D = 0.2 m2/s, six cell-centered values, and zero flux at both ends. With h = L/6, write du/dt = −K u, where K = (D/h2) times the symmetric tridiagonal matrix with entries −1 next to the diagonal, 2 on the interior diagonal, and 1 at both diagonal endpoints. Find all six eigenvalues and eigenvectors with Jacobi or a symmetric library solver. Use the first-cell pulse u(0) = [1,0,0,0,0,0] and project it onto the modes to predict its subsequent relaxation. Try the power method from this pulse for the fastest-decaying mode.

Validation: the eigenvalues are λj = 4D/h2 sin2(jπ/12), for j = 0,1,...,5. In particular, λ0 = 0, λ1 ≈ 1.929234 s−1, and λ5 ≈ 26.870766 s−1. The constant mode stays fixed, so the initial pulse tends toward u = [1/6,1/6,1/6,1/6,1/6,1/6] and its total concentration remains 1 in these discrete units. Check ||K*v − λv||2 for each mode and verify their orthogonality.

Possible extension: double D and verify that every nonzero decay rate doubles while the stationary mode remains unchanged.

← Back to Chapter 8