Solving Linear Equation Systems with SciPy

Linear equation systems are collections of linear equations that share the same set of variables. They are fundamental in both mathematics and real-world applications, such as engineering, physics, and economics, where multiple relationships can be modeled simultaneously. For instance, when dealing with electrical circuits, mechanical systems, or economic models, you'll often encounter linear systems that need solutions.

In this lesson, we'll explore how to solve linear equation systems using Python, with a focus on leveraging the SciPy library.

Before we begin, let's recall some basics. You might remember that NumPy provides a way to work with arrays and matrices, which are essential for representing and solving linear systems.

SciPy builds on top of NumPy and provides a suite of functions specifically designed for scientific and technical computing.

Let's refresh it a bit. A linear equation in its general form can be written as:

where $a_1, a_2, \ldots, a_n$ are coefficients, $x_1, x_2, \ldots, x_n$ are variables, and $b$ is a constant term.

In matrix form, a system of linear equations can be represented as:

where:

• A is the coefficient matrix.
• x is the column vector of variables (unknowns).
• b is the constant terms vector.

Here's an example of formulating a 3x3 linear system:

Consider the following equations:

We can represent these equations in matrix form as:

where A contains coefficients for the variables x, y, z, and b contains free coefficients.

The solve function from the scipy.linalg module is used to find the solution vector x for the equation $Ax = b$. Let's walk through the steps involved.

First, we need to import the required modules. These are scipy.linalg for the solve function and numpy for handling arrays.

from scipy.linalg import solve
import numpy as np

We define the matrix A and vector b using NumPy arrays, which will represent our system of equations.

A = np.array([[3, 2, -1],
              [2, -2, 4],
              [-1, 0.5, -1]])

b = np.array([1, -2, 0])

We use the solve function to find the solution vector x.

x = solve(A, b)

The solve function works by calculating the determinant of matrix A. If det(A) ≠ 0, it finds the inverse of A and solves for x using the formula:

If det(A) = 0, the matrix A is singular, meaning it does not have an inverse, and the system may have no solutions or infinitely many solutions. In such cases, the solve function will raise a LinAlgError.

Finally, we print the solution to see the result of our computation.

print("Solution for vector x:", x)  # Solution for vector x: [1.  -2.  -2.]