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?
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.