Menu
Ran a simple script in Numpy to solve a system of 3 equations:
from numpy import *
a = matrix('1 4 1; 4 13 7; 7 22 13')
b = matrix('0;0;1')
print linalg.solve(a,b)
But when I ran it through the command prompt, I somehow got:
[[ 3.46430741e+15]
[ -6.92861481e+14]
[ -6.92861481e+14]]
although Wolfram Alpha says there's no solution.
Does anyone know why there seems to be this difference between the Numpy/WRA answers?
670
From the previous section, we know that to solve a system of linear equations, we need to perform two operations: matrix inversion and a matrix dot product. The Numpy library from Python supports both the operations. If you have not already installed the Numpy library, you can do with the following pip command:
Thankfully, in Python, we have NumPy library to solve this matrix and vector calculation. By using NumPy library, we can make that 3-lines of code become 1-line of code like this one below, And if we compare the time, for this case, the conventional way will take around 0.000224 and the NumPy method is just 0.000076.
Let's first create the matrix A in Python. To create a matrix, the array method of the Numpy module can be used. A matrix can be considered as a list of lists where each list represents a row. In the following script we create a list named m_list, which further contains two lists: [4,3] and [-5,9]. These lists are the two rows in the matrix A.
There are multiple ways to solve such a system, such as Elimination of Variables, Cramer's Rule, Row Reduction Technique, and the Matrix Solution. In this article we will cover the matrix solution.
If you take pen and paper (or use numpy.linalg.matrix_rank) and calculate the ranks of the coefficient and the augmented matrix you'll see that, according to Rouché–Capelli theorem, there is no solution. So Wolfram is right.
NumPy searches for numerical solution of equation using LU decomposition. Without going into details LU decomposition involves divisions, divisions in FP arithmetic can lead to significant errors.
If you check:
a = np.matrix('1 4 1; 4 13 7; 7 22 13')
b = np.matrix('0;0;1')
c = np.linalg.solve(a,b)
np.linalg.norm(a * c - b)
## 2.0039024427351748
you'll see that NumPy solution is far from being correct.
95
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With
