Solving Linear Systems¶

This notebook demonstrates solving $Ax = b$ using the symbolic-numeric bridge: we set up systems of equations symbolically with Maxima's matrix algebra, then solve them numerically with np_solve (wrapping LAPACK). Along the way we explore residuals, condition numbers, ill-conditioning, and least-squares solutions.

Key functions: ndarray, np_solve, np_matmul, np_inv, np_norm, np_lstsq.

In [121]:
load("numerics")$
load("ax-plots")$

Circuit Analysis: Mesh Current Method [S+N]¶

Consider a resistor network with three meshes (loops). Using Kirchhoff's Voltage Law (KVL), each mesh gives one equation relating the mesh currents $I_1, I_2, I_3$ to the source voltages.

The circuit has:

  • $R_1 = 10\,\Omega$, $R_2 = 20\,\Omega$, $R_3 = 30\,\Omega$, $R_4 = 15\,\Omega$, $R_5 = 25\,\Omega$
  • $V_1 = 5\,\text{V}$ (mesh 1), $V_2 = 10\,\text{V}$ (mesh 3)

The mesh equations from KVL:

$$ (R_1 + R_2)\,I_1 - R_2\,I_2 = V_1 $$ $$ -R_2\,I_1 + (R_2 + R_3 + R_4)\,I_2 - R_4\,I_3 = 0 $$ $$ -R_4\,I_2 + (R_4 + R_5)\,I_3 = -V_2 $$

We set up the equations symbolically, extract the coefficient matrix and RHS vector, then solve numerically.

In [122]:
/* Define resistor and voltage values */
R1 : 10$ R2 : 20$ R3 : 30$ R4 : 15$ R5 : 25$
V1 : 5$ V2 : 10$

/* Mesh equations: LHS - RHS = 0 */
eq1 : (R1 + R2)*I1 - R2*I2 = V1;
eq2 : -R2*I1 + (R2 + R3 + R4)*I2 - R4*I3 = 0;
eq3 : -R4*I2 + (R4 + R5)*I3 = -V2;

print("Mesh equations:")$
print(eq1)$
print(eq2)$
print(eq3)$

30*I1-20*I2 = 5
-(15*I3)+65*I2-20*I1 = 0
40*I3-15*I2 = -10
Mesh equations:
30*I1-20*I2 = 5
-(15*I3)+65*I2-20*I1 = 0
40*I3-15*I2 = -10
In [123]:
/* Extract the coefficient matrix and RHS symbolically */
A : coefmatrix([eq1, eq2, eq3], [I1, I2, I3]);
b_rhs : makelist(rhs(eq), eq, [eq1, eq2, eq3])$
b_vec : transpose(matrix(b_rhs));

print("Coefficient matrix A:")$
A;
print("RHS vector b:")$
b_vec;

matrix([30,-20,0],[-20,65,-15],[0,-15,40])
matrix([5],[0],[-10])
Coefficient matrix A:
matrix([30,-20,0],[-20,65,-15],[0,-15,40])
RHS vector b:
Out[123]:
In [124]:
/* Convert to ndarrays and solve */
A_nd : ndarray(A)$
b_nd : ndarray(b_vec)$

x : np_solve(A_nd, b_nd)$

currents : np_to_list(x)$
print("Mesh currents (A):")$
print("  I1 =", currents[1])$
print("  I2 =", currents[2])$
print("  I3 =", currents[3])$

Mesh currents (A):
  I1 = 0.16063348416289594
  I2 = -0.009049773755656094
  I3 = -0.253393665158371
In [125]:
/* Verify: solve symbolically with Maxima and compare */
sym_sol : solve([eq1, eq2, eq3], [I1, I2, I3]);
print("Symbolic solution:")$
for s in sym_sol[1] do print("  ", s, " =", float(rhs(s)))$

[[I1 = 71/442,I2 = -(2/221),I3 = -(56/221)]]
Symbolic solution:
   I1 = 71/442  = 0.16063348416289594
   I2 = -(2/221)  = -0.00904977375565611
   I3 = -(56/221)  = -0.25339366515837103

Residual Analysis¶

After solving $Ax = b$, we check the quality of the solution by computing the residual $r = Ax - b$. In exact arithmetic $r = 0$; in floating-point we expect $\|r\|$ near machine epsilon times $\|A\| \|x\|$.

In [126]:
/* Compute residual r = A*x - b */
r : np_sub(np_matmul(A_nd, x), b_nd)$

print("Residual vector:", np_to_list(r))$
print("Residual norm ||r||_2 =", np_norm(r))$
print("Machine epsilon =", 1.1e-16)$
print("||A|| * ||x|| =", np_norm(A_nd) * np_norm(x))$

Residual vector: [0.0,-2.220446049250313e-16,1.7763568394002505e-15]
Residual norm ||r||_2 = 1.790180836524724e-15
Machine epsilon = 1.1e-16
||A|| * ||x|| = 26.804749193926753

The residual norm is near machine epsilon, confirming that LAPACK's solver delivers a solution accurate to the limits of double-precision arithmetic.

Condition Number [S+N]¶

The condition number $\kappa(A) = \|A\| \cdot \|A^{-1}\|$ quantifies how sensitive the solution $x$ is to perturbations in $A$ or $b$.

  • $\kappa \approx 1$: well-conditioned (small errors stay small)
  • $\kappa \gg 1$: ill-conditioned (small errors are amplified)
  • $\kappa \sim 10^k$: roughly $k$ digits of accuracy are lost
In [127]:
/* Condition number of the circuit matrix */
A_inv : np_inv(A_nd)$
kappa : np_norm(A_nd) * np_norm(A_inv);

print("Condition number kappa(A) =", kappa)$
print("This is a well-conditioned system.")$

5.444965332150018
Condition number kappa(A) = 5.444965332150018
This is a well-conditioned system.

Ill-Conditioned System: Hilbert Matrix¶

The Hilbert matrix $H_{ij} = \frac{1}{i+j-1}$ is a classic example of ill-conditioning. Its condition number grows exponentially with size $n$, making numerical solution increasingly unreliable.

We build Hilbert matrices of increasing size, compute their condition numbers, and show how the solution error grows.

In [128]:
/* Build a Hilbert matrix of size n */
hilbert(n) := makelist(
  makelist(1.0/(i + j - 1), j, 1, n),
  i, 1, n
)$

/* Example: 5x5 Hilbert matrix */
H5 : apply(matrix, hilbert(5));
H5_nd : ndarray(H5)$

/* True solution: x = [1, 1, 1, 1, 1] */
x_true : ndarray([1, 1, 1, 1, 1], [5, 1])$
b5 : np_matmul(H5_nd, x_true)$

/* Solve and check error */
x_computed : np_solve(H5_nd, b5)$
error5 : np_sub(x_computed, x_true)$

kappa5 : np_norm(H5_nd) * np_norm(np_inv(H5_nd))$

print("Hilbert 5x5:")$
print("  Condition number =", kappa5)$
print("  Solution error ||x - x_true|| =", np_norm(error5))$
print("  Computed x =", np_to_list(x_computed))$

matrix([1.0,0.5,0.3333333333333333,0.25,0.2],
       [0.5,0.3333333333333333,0.25,0.2,0.16666666666666666],
       [0.3333333333333333,0.25,0.2,0.16666666666666666,0.14285714285714285],
       [0.25,0.2,0.16666666666666666,0.14285714285714285,0.125],
       [0.2,0.16666666666666666,0.14285714285714285,0.125,0.1111111111111111])
Hilbert 5x5:
  Condition number = 480849.11699433613
  Solution error ||x - x_true|| = 2.6493798901547713e-11
  Computed x =
              [1.0000000000001708,0.9999999999968402,1.0000000000135154,
               0.9999999999797115,1.00000000000988]
In [129]:
/* Condition number vs matrix size */
sizes : [2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]$
kappas : []$

for n in sizes do block(
  [Hn, Hn_nd, kn],
  Hn : apply(matrix, hilbert(n)),
  Hn_nd : ndarray(Hn),
  kn : np_norm(Hn_nd) * np_norm(np_inv(Hn_nd)),
  kappas : append(kappas, [kn])
)$

/* Display results */
for i : 1 thru length(sizes) do
  print(sconcat("  n=", sizes[i], ": kappa = ", kappas[i]))$

  n=2: kappa = 19.333333333333336
  n=3: kappa = 526.1588210797213
  n=4: kappa = 15613.793559642276
  n=5: kappa = 480849.11699433613
  n=6: kappa = 1.5118987126990862e7
  n=7: kappa = 4.817472525639491e8
  n=8: kappa = 1.5493617168394405e10
  n=9: kappa = 5.017293210077868e11
  n=10: kappa = 1.6331839891317984e13
  n=11: kappa = 5.326131742171254e14
  n=12: kappa = 1.6686299298869728e16
In [130]:
/* Plot: condition number grows exponentially with n */
cond_data : makelist([sizes[i], float(log(kappas[i])/log(10))], i, 1, length(sizes))$

ax_draw2d(
  color=blue, marker_size=6,
  points(cond_data),
  color=blue,
  lines(cond_data),
  title="Hilbert Matrix: Condition Number vs Size",
  xlabel="Matrix size n",
  ylabel="log10(kappa)",
  grid=true
)$
No description has been provided for this image

The condition number grows roughly as $\kappa \sim e^{3.5n}$. By $n = 12$ the condition number exceeds $10^{16}$, which means all digits of a double-precision solution may be corrupted. This is why algorithms like iterative refinement and higher-precision arithmetic exist.

Truss Analysis [S+N]¶

A simple 2D truss with 3 bars and 2 free nodes. We set up the equilibrium equations from statics: at each node, the sum of forces in $x$ and $y$ must be zero.

Geometry:

  • Node A at $(0, 0)$ — pinned support
  • Node B at $(4, 0)$ — roller support (vertical reaction only)
  • Node C at $(2, 3)$ — free node, external load $P = 100\,\text{N}$ downward

Members: AC, BC, AB

We solve for the member forces $F_{AC}$, $F_{BC}$, and $F_{AB}$ using the method of joints at node C.

In [131]:
/* Truss geometry */
/* Angles of members from node C */
theta_AC : atan2(-3, -2)$  /* C to A */
theta_BC : atan2(-3, 2)$   /* C to B */

/* Equilibrium at node C: sum Fx = 0, sum Fy = 0 */
/* F_AC along CA direction, F_BC along CB direction */
/* External load: P = 100 N downward */
P : 100$

eq_x : F_AC * cos(theta_AC) + F_BC * cos(theta_BC) = 0;
eq_y : F_AC * sin(theta_AC) + F_BC * sin(theta_BC) = P;

print("Equilibrium equations at node C:")$
print(eq_x)$
print(eq_y)$

/* Solve symbolically */
truss_sol : solve([eq_x, eq_y], [F_AC, F_BC]);
print("Member forces (symbolic):")$
print(truss_sol)$

/* Numeric values */
print("F_AC =", float(rhs(truss_sol[1][1])), "N")$
print("F_BC =", float(rhs(truss_sol[1][2])), "N")$

(2*F_BC)/sqrt(13)-(2*F_AC)/sqrt(13) = 0
-((3*F_BC)/sqrt(13))-(3*F_AC)/sqrt(13) = 100
Equilibrium equations at node C:
(2*F_BC)/sqrt(13)-(2*F_AC)/sqrt(13) = 0
-((3*F_BC)/sqrt(13))-(3*F_AC)/sqrt(13) = 100
[[F_AC = -((50*sqrt(13))/3),F_BC = -((50*sqrt(13))/3)]]
Member forces (symbolic):
[[F_AC = -((50*sqrt(13))/3),F_BC = -((50*sqrt(13))/3)]]
F_AC = -60.092521257733154 N
F_BC = -60.092521257733154 N
In [132]:
/* Verify numerically with np_solve */
A_truss : float(coefmatrix([eq_x, eq_y], [F_AC, F_BC]))$
b_truss : float(transpose(matrix(makelist(rhs(eq), eq, [eq_x, eq_y]))))$

A_truss_nd : ndarray(A_truss)$
b_truss_nd : ndarray(b_truss)$

forces : np_solve(A_truss_nd, b_truss_nd)$
print("Numeric solution:", np_to_list(forces))$

Numeric solution: [-60.092521257733154,-60.092521257733154]

Over-Determined System: Least Squares¶

When there are more equations than unknowns ($m > n$), the system $Ax = b$ is typically inconsistent. The least-squares solution minimises $\|Ax - b\|_2$ and is computed via SVD by np_lstsq.

Example: fitting a line $y = c_0 + c_1 x$ to 5 data points.

In [133]:
/* 5 data points: (x, y) */
data_x : [1.0, 2.0, 3.0, 4.0, 5.0]$
data_y : [2.1, 3.9, 6.2, 7.8, 10.1]$

/* Design matrix: [1, x_i] for each data point */
A_ls : apply(matrix,
  makelist([1.0, data_x[i]], i, 1, length(data_x))
)$
b_ls : transpose(matrix(data_y))$

print("Design matrix A (5x2):")$
A_ls;
print("Observations b:")$
b_ls;

Design matrix A (5x2):
matrix([1.0,1.0],[1.0,2.0],[1.0,3.0],[1.0,4.0],[1.0,5.0])
Observations b:
Out[133]:
In [134]:
/* Least-squares solve */
A_ls_nd : ndarray(A_ls)$
b_ls_nd : ndarray(b_ls)$

[x_ls, residuals, rank_ls, S_ls] : np_lstsq(A_ls_nd, b_ls_nd)$

coeffs : np_to_list(x_ls)$
print("Least-squares coefficients:")$
print("  intercept c0 =", coeffs[1])$
print("  slope     c1 =", coeffs[2])$
print("  rank =", rank_ls)$

Least-squares coefficients:
  intercept c0 = 0.04999999999999677
  slope     c1 = 1.9900000000000002
  rank = 2
In [135]:
/* Plot data and fitted line */
c0 : coeffs[1]$
c1 : coeffs[2]$
fitted_line : c0 + c1*x$

ax_draw2d(
  color=red, marker_size=8,
  name="Data",
  points(data_x, data_y),
  color=blue, name="Least-squares fit",
  explicit(fitted_line, x, 0, 6),
  title="Over-Determined System: Line Fit via np_lstsq",
  xlabel="x", ylabel="y",
  grid=true
)$
No description has been provided for this image