Lesson 10 - From Balance Laws to Linear Systems¶
Learning goals. By the end of this lecture, students should be able to:
translate a small steady balance model into simultaneous linear equations,
identify the physical meaning of the unknown vector, coefficient matrix, and right-hand side,
assemble a matrix equation of the form \(\mathbf{A}\mathbf{x}=\mathbf{b}\),
use a black-box numerical solve and check the residual,
interpret an inverse conceptually as a source-to-solution response matrix,
explain why a coupled system must be solved as a system rather than one row at a time.
Key terms. balance law, unknown vector, coefficient matrix, right-hand side, residual.
Why linear systems appear in engineering models¶
Many engineering calculations begin with a balance. Over a control volume, cell, energy group, or state, the balance has the schematic form
If the model is steady, the accumulation term is zero or prescribed. If the rates are proportional to the quantities being modeled, the result is a linear equation.
One balance is useful. Several balances coupled together are more useful. If region 1 exchanges material, energy, or particles with region 2, the unknown \(x_1\) in region 1 appears in the balance for region 2, and the unknown in region 2 \(x_2\) appears in the balance for region 1. Systems for which this coupling is linear are the focus of linear algebra, the language of which includes vectors and matrices.
The basic object of this chapter is therefore
The vector \(\mathbf{x}\) contains all the unknown quantities. Each row of \(\mathbf{A}\) represents one balance equation. The right-hand side vector \(\mathbf{b}\) contains sources, requirements, boundary data, or other known forcing terms.
A fixed-step walk with absorption¶
To make the idea concrete, consider a one-dimensional line of four cells. A source injects particles into the first cell. During each step, a fraction \(c\) of the population survives and scatters to neighboring cells. The remaining fraction \(1-c\) is absorbed or lost. At the ends of the domain, particles that step out of the system leak away.
Let
Let the source be
For this system, the source vector is the right-hand side, so later we will write \(\mathbf{b}=\mathbf{s}\). In this toy model the entries are source-normalized particle counts: \(\mathbf{s}\) represents one injected source particle, \(x_i\) is the expected number of visits or population units per source particle, and \(c\) is dimensionless.
The first cell receives the external source and particles scattered from cell 2. Since only half of the surviving particles in cell 2 move left, the balance is
The second cell receives particles scattered from cells 1 and 3:
Likewise,
These four equations are coupled. For example, \(x_2\) depends on \(x_1\) and \(x_3\), but \(x_1\) and \(x_3\) also depend on \(x_2\). Guessing one unknown at a time is not the natural approach; solving all four equations together is.
A boundary-value variant¶
The same algebraic balance idea also appears with dimensionless probabilities rather than dimensional populations. Suppose a particle moves on the integer sites \(0,1,\ldots,N\). Let
The \(u_i\) are dimensionless probabilities or source-normalized responses.
The boundary values are known:
If an interior collision scatters left or right with equal probability and the scattering ratio is \(c\), then the interior balance is
Moving unknowns to the left gives
For interior nodes away from the boundary, this produces the same tridiagonal pattern seen below: a diagonal coefficient and two nearest-neighbor coefficients. At the right boundary, however, the known value \(u_N=1\) is moved to the right-hand side:
Thus boundary data can appear in \(\mathbf{b}\) even when the original interior balance has no external source term.
Takeaway.
Boundary conditions are known values. When a known boundary value appears in an equation for an unknown interior value, move that known contribution to the right-hand side.
Worked example: boundary data in the right-hand side.
Problem. Take \(N=4\), so the unknown responses are \(u_1,u_2,u_3\), with
\[u_0=0,\qquad u_4=1.\]Use \(c=0.8\). Assemble the matrix equation and identify where the boundary data enters the right-hand side.
Solution. Since \(c/2=0.4\), the interior equations are
\[\begin{split}\begin{alignedat}{3} u_1 &= 0.4u_0+0.4u_2, &\qquad& u_1-0.4u_2 &= 0,\\ u_2 &= 0.4u_1+0.4u_3, &\qquad& -0.4u_1+u_2-0.4u_3 &= 0,\\ u_3 &= 0.4u_2+0.4u_4, &\qquad& -0.4u_2+u_3 &= 0.4. \end{alignedat}\end{split}\]The right boundary value \(u_4=1\) is known, so its contribution becomes the nonzero entry in the right-hand side:
\[\begin{split}\underbrace{\begin{bmatrix} 1 & -0.4 & 0\\ -0.4 & 1 & -0.4\\ 0 & -0.4 & 1 \end{bmatrix}}_{\mathbf{A}} \underbrace{\begin{bmatrix}u_1\\u_2\\u_3\end{bmatrix}}_{\mathbf{u}} = \underbrace{\begin{bmatrix}0\\0\\0.4\end{bmatrix}}_{\mathbf{b}}.\end{split}\]Solving gives
\[\begin{split}\mathbf{u}\approx \begin{bmatrix} 0.0941\\ 0.2353\\ 0.4941 \end{bmatrix}.\end{split}\]The response increases as the starting site moves toward the right boundary. The values are below the pure random-walk probabilities \(1/4,1/2,3/4\) because \(c<1\) leaves room for absorption or loss before either boundary is reached.
The units also check out. The matrix entries \(1\) and \(-0.4\) are dimensionless probabilities or survival fractions, the boundary responses \(u_0\) and \(u_4\) are dimensionless, and every entry of \(\mathbf{b}\) has the same dimensionless response units as \(\mathbf{u}\).
Putting the balances into matrix form¶
Move all unknown terms to the left-hand side. Then rewrite each equation so that every unknown appears in the same order, even when its coefficient is zero:
Now the coefficient matrix is visible:
This matrix has a clear physical interpretation. The diagonal entries multiply the population being balanced in each cell. The off-diagonal entries describe coupling between neighboring cells. The zeros say that non-neighboring cells do not exchange particles directly in this model. In this source-normalized example, \(\mathbf{A}\) is dimensionless and \(\mathbf{x}\), \(\mathbf{b}\), and the residual \(\mathbf{r}=\mathbf{b}-\mathbf{A}\mathbf{x}\) all use the same population-per-source-particle units. In a dimensional reactor balance matrix, entries may instead contain terms such as \(V_i\Sigma_{a,i}\), leakage coefficients, or scattering/removal rates with units chosen so each row matches \(\mathbf{b}\).
If \(c=0.8\), then
Solving first, opening the box later¶
For now, we use np.linalg.solve as a “black box” tool - we give it an input, it gives us an output, and we ignore the details. Later chapters explain the operations and factorizations behind this call.
With \(c=0.8\), the example above can be defined and solved via
[1]:
import numpy as np
c = 0.8
A = np.array([[1.0, -c/2, 0.0, 0.0],
[-c/2, 1.0, -c/2, 0.0],
[0.0, -c/2, 1.0, -c/2],
[0.0, 0.0, -c/2, 1.0]])
b = np.array([1.0, 0.0, 0.0, 0.0])
x = np.linalg.solve(A, b)
r = b - A @ x
print(x)
print(np.linalg.norm(r))
[1.24633431 0.61583578 0.29325513 0.11730205]
0.0
The computed solution is approximately
The values decrease as we move away from the source, which matches the physical picture. The source enters cell 1 directly. Particles reach cells 2, 3, and 4 only after one or more scattering events, and absorption reduces the expected population along the way.
Residuals as a first quality check¶
The computed \(\mathbf{x}\) is only useful if it nearly satisfies the original equations. The residual is
If \(\mathbf{r}=\mathbf{0}\), the equations are satisfied exactly. In floating-point arithmetic, we usually ask whether \(\mathbf{r}\) is small. Each component of \(\mathbf{r}\) is a defect in the balance equations. The residual has the same units as \(\mathbf{b}\). In this random-walk example it is source-normalized; in a dimensional source balance it would have source or rate units.
For the random-walk example, the first residual component is
Thus \(r_1\) directly checks the first cell balance. The other components check the other cells. This is the first important diagnostic for numerical linear algebra: after computing a candidate solution, put it back into the model.
Takeaway.
The residual \(\mathbf{r}=\mathbf{b}-\mathbf{A}\mathbf{x}\) measures how much the computed solution fails to satisfy the original equations. Its components are defects in the balance equations.
A second view: source-to-solution response¶
Now we can introduce the inverse idea with minimal extra machinery. Keep \(\mathbf{A}\) fixed. The solver takes a source pattern \(\mathbf{b}\) and returns the corresponding solution \(\mathbf{x}\):
For this lecture, call that mapping the response map.
Returning to the source-normalized population model, the easiest way to see the response map is to use unit source patterns. Define
Let \(\mathbf{m}_j\) be the solution caused by the unit source \(\mathbf{e}_j\):
Thus \(\mathbf{m}_1\) is the response to a unit source in cell 1, \(\mathbf{m}_2\) is the response to a unit source in cell 2, and so on.
For \(c=0.8\), these four response vectors are approximately the columns of
The first column is exactly the solution we computed for the original source \(\mathbf{b}=\mathbf{e}_1\).
Why does this help with a general source? Because the equations are linear. A source vector such as
can be built from unit sources. The corresponding solution is the same combination of unit-source responses:
The compact notation for this weighted combination is
The coefficient matrix now has a physical role. Given a candidate state \(\mathbf{x}\), the product \(\mathbf{A}\mathbf{x}\) evaluates all of the balances and tells us what source or requirement vector would be needed to support that state:
But the usual engineering question runs the other way. We are given the source, forcing, or boundary data \(\mathbf{b}\), and we want the state \(\mathbf{x}\) that the system produces. We can represent that opposite direction by a response map:
In words,
When a square system has one solution for every right-hand side, this response map is the inverse of the balance operator:
The inverse is therefore best understood first as a map from inputs to responses. The units of a response-map entry are the units of solution divided by the units of source; in this normalized example the entries are dimensionless. It is useful notation, but it is usually not the object we form in a computation. In practice, we solve \(\mathbf{A}\mathbf{x}=\mathbf{b}\) directly and check the residual \(\mathbf{b}-\mathbf{A}\mathbf{x}\).
[2]:
responses = []
for j in range(4):
e = np.zeros(4)
e[j] = 1.0
responses.append(np.linalg.solve(A, e))
M = np.column_stack(responses)
print(M)
# Original b = e_1, so its solution is column 0.
print(M[:, 0])
[[1.24633431 0.61583578 0.29325513 0.11730205]
[0.61583578 1.53958944 0.73313783 0.29325513]
[0.29325513 0.73313783 1.53958944 0.61583578]
[0.11730205 0.29325513 0.61583578 1.24633431]]
[1.24633431 0.61583578 0.29325513 0.11730205]
Takeaway.
A good first reading is: \(\mathbf{A}\) maps a candidate solution to the source pattern required by the balances, while \(\mathbf{M}\) maps a source pattern to the solution. The inverse is the response map, not the recommended way to compute the solution.
Exercises¶
For the four-cell fixed-step walk with \(c=0.8\), write the matrix \(\mathbf{A}\) for the source-normalized population model.
In the boundary-value variant, let \(u_0=0\), \(u_4=1\), and \(c=0.8\). Explain why the right boundary contributes \(0.4\) to the third entry of \(\mathbf{b}\).
For
\[\begin{split}\mathbf{A}=\begin{bmatrix} 1 & -0.4 & 0\\ -0.4 & 1 & -0.4\\ 0 & -0.4 & 1 \end{bmatrix}, \qquad \mathbf{x}=\begin{bmatrix}1\\1\\1\end{bmatrix}, \qquad \mathbf{b}=\begin{bmatrix}0.6\\0.2\\0.6\end{bmatrix},\end{split}\]compute the residual \(\mathbf{r}=\mathbf{b}-\mathbf{A}\mathbf{x}\).
If \(\mathbf{M}=\mathbf{A}^{-1}\), what does column 2 of \(\mathbf{M}\) represent in a source-response interpretation?
Answers.
\(\mathbf{A}\) has \(1\) on the diagonal and \(-0.4\) on the neighboring off-diagonals:
\[\begin{split}\begin{bmatrix} 1&-0.4&0&0\\ -0.4&1&-0.4&0\\ 0&-0.4&1&-0.4\\ 0&0&-0.4&1 \end{bmatrix}.\end{split}\]The right boundary term enters the balance for cell 3 as \((c/2)u_4=(0.8/2)(1)=0.4\). It is known boundary data, so it moves to the right-hand side.
\(\mathbf{A}\mathbf{x}=[0.6,0.2,0.6]^{\mathsf{T}}\), so \(\mathbf{r}=\mathbf{0}\).
Column 2 is the solution produced by a unit source in entry 2, with all other source entries zero. It is a response column, not a probability column.
End-of-lecture summary¶
Takeaway.
Steady linear balance laws and boundary data lead naturally to systems of equations.
The unknown vector stores the quantities being solved for, one per cell, region, group, or state.
Each row of the matrix represents one balance equation.
np.linalg.solvegives a practical first solve, but the residual \(\mathbf{b}-\mathbf{A}\mathbf{x}\) should always be checked as a defect in the balance equations.The inverse can be read conceptually as a response matrix whose columns are unit-source solutions.
Lesson 11 opens the black box by explaining the dot products and matrix-vector products behind \(\mathbf{A}\mathbf{x}\).