Hands-On Practical: Solving Systems with NumPy

Theory provides the map, but practical application is where you learn the terrain. A system of linear equations is represented as $Ax = b$, and the matrix inverse, $A^{-1}$, helps find the solution vector $x$. You will now translate this mathematical understanding into code using NumPy, the standard library for numerical computing in Python.

While computing the inverse is a valid method, NumPy offers a more direct, efficient, and numerically stable function to solve for $x$ directly. Let's work through a practical example.

Setting Up the Problem

Imagine you're buying fruit. On the first day, you buy 2 apples and 3 bananas, and your total is $8. On the second day, you buy 4 apples and 1 banana for a total of$9. If the prices remain the same, what is the price of a single apple and a single banana?

We can set this up as a system of two linear equations with two unknowns: Let $x_1$ be the price of an apple. Let $x_2$ be the price of a banana.

The system is:

$$2x_1 + 3x_2 = 8$$ $$4x_1 + 1x_2 = 9$$

In matrix form, $Ax = b$, this becomes:

$$\begin{pmatrix} 2 & 3 \\ 4 & 1 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \end{pmatrix} = \begin{pmatrix} 8 \\ 9 \end{pmatrix}$$

Our goal is to find the vector $x$, which contains the prices of an apple and a banana.

Solving the System with NumPy

First, let's import NumPy and create our matrix A and vector b.

import numpy as np

# The matrix of coefficients (A)
A = np.array([
    [2, 3],
    [4, 1]
])

# The vector of outcomes (b)
b = np.array([8, 9])

print("Matrix A:\n", A)
print("\nVector b:\n", b)

Now, we can use the numpy.linalg.solve() function. This function is designed specifically to solve equations of the form $Ax = b$. It takes the matrix A and the vector b as arguments and returns the solution vector x.

# Solve for x
x = np.linalg.solve(A, b)

print("\nSolution vector x:\n", x)

The output shows that $x_1 = 1.9$ and $x_2 = 1.4$. In the context of our problem, this means an apple costs $1.90, and a banana costs$1.40.

Verifying the Solution

A good practice is to verify that our solution is correct. If our x is right, then multiplying it by our original matrix A should give us our original outcome vector b.

Let's calculate $Ax$ and see if it equals $b$.

# Verify the solution by calculating A @ x
# The @ operator performs matrix multiplication
verification = A @ x

print("Verification (A @ x):\n", verification)
print("\nOriginal vector b:\n", b)

# Check if the verification is close to b
# This is useful for handling tiny floating-point errors
is_close = np.allclose(verification, b)
print("\nIs the solution correct?", is_close)

The result of the matrix multiplication A @ x is [8. 9.], which matches our original b vector exactly. The function np.allclose() is a reliable way to compare floating-point arrays, returning True if they are element-wise equal within a small tolerance.

The Inverse Method vs. solve()

In the previous sections, we discussed solving for $x$ by finding the inverse of $A$ and computing $x = A^{-1}b$. You can do this in NumPy as well.

# Find the inverse of A
A_inv = np.linalg.inv(A)

# Calculate x using the inverse
x_from_inverse = A_inv @ b

print("Solution from inverse:\n", x_from_inverse)

You get the same result. So why does np.linalg.solve() exist? For large matrices, np.linalg.solve() is both faster and more numerically stable than computing the inverse and then performing matrix multiplication. The solve() function uses more advanced decomposition techniques that avoid some of the challenges of direct inversion. As a general rule, you should prefer np.linalg.solve() over np.linalg.inv().

Visualizing the Solution

For a 2D system like ours, the solution represents the point where two lines intersect. We can plot the two equations to see this visually.

The first equation, $2x_1 + 3x_2 = 8$, can be rewritten as $x_2 = (8 - 2x_1) / 3$. The second equation, $4x_1 + x_2 = 9$, can be rewritten as $x_2 = 9 - 4x_1$.

The plot below shows these two lines and their intersection point, which is our solution.

The solution to the system is the unique point where the two lines cross. This geometric interpretation is a useful way to think about what solving a system of equations accomplishes.

Handling Non-Invertible (Singular) Systems

What happens if a system has no unique solution? This occurs when the matrix A is singular (non-invertible), which you learned means its determinant is zero. Let's see how NumPy handles this.

Consider this system:

$$2x_1 + 3x_2 = 8$$ $$4x_1 + 6x_2 = 16$$

The second equation is just two times the first one. Geometrically, these are the same line, so there are infinite solutions. The corresponding matrix is singular.

# A singular matrix (second row is 2 * first row)
A_singular = np.array([
    [2, 3],
    [4, 6]
])

b_singular = np.array([8, 16])

# Let's try to solve it
try:
    x_singular = np.linalg.solve(A_singular, b_singular)
except np.linalg.LinAlgError as e:
    print("Error:", e)

When you run this code, NumPy correctly identifies the problem and raises a LinAlgError with the message "Singular matrix". This is NumPy's way of telling you that the system you've defined cannot be solved to find a single, unique solution vector x. This practical feedback is extremely useful, preventing your programs from producing nonsensical results when faced with an unsolvable system.