I'm trying to learn Python for basic work in linear algebra. I'm running into the following problem with a simple system of linear equations:
import scipy.linalg as la
import numpy as np
A = np.array([[186/450, 54/21, 30/60],
[12/450, 6/21 , 3/60],
[9/450, 6/21 , 15/60]])
l = np.array([18/450, 12/21, 30/36])
b = np.array([2, 0, 1/6])
y = np.array([180, 0, 30])
x = la.inv(np.eye(3) - A) @ y
lam = np.transpose(l) @ la.inv(np.eye(3) - A)
This returns
array([0.21212121, 2.12121212, 1.39393939])
which is incorrect. Performing the same operation in Julia,
A = [186/450 54/21 30/60;
12/450 6/21 3/60;
9/450 6/21 15/60]
l = [18/450, 12/21, 30/60]
b = [2, 0, 1/6]
y = [180, 0, 30]
λ = l' * inv(I - A)
yields the correct result, which is
1×3 adjoint(::Vector{Float64}) with eltype Float64:
0.181818 1.81818 0.909091
What am I missing here? I think I might be missing something in the opaque numpy array syntax.
CodePudding user response:
There is a typo in 'l' instantiation in your python code. (30/36 should be 30/60).
This code with the typo fixed produces the same result as in Julia.
import scipy.linalg as la
import numpy as np
A = np.array([[186/450, 54/21, 30/60],
[12/450, 6/21 , 3/60],
[9/450, 6/21 , 15/60]])
l = np.array([18/450, 12/21, 30/60]) #typo fixed here
b = np.array([2, 0, 1/6])
y = np.array([180, 0, 30])
x = la.inv(np.eye(3) - A) @ y
lam = np.transpose(l) @ la.inv(np.eye(3) - A)
Giving:
array([0.18181818, 1.81818182, 0.90909091])