2

我有一组 6 个方程，我想让 numpy 为我求解。所以我构造了一个 6x6 的系数矩阵，并用各种值填充它。但是，我最终为此编写的代码非常难以辨认，并且几乎没有向我的代码读者传达我想要解决的方程式。

例如，填写系数矩阵如下所示：

# Coefficients matrix
# Order of variables: w, X, Y, Z, s, t
A = np.mat( np.zeros((6,6)) )

A[0:3,0] = cam_inv[...,2]
A[0:3,1:4] = -np.identity(3)
A[3:6,1:4] = np.identity(3)
A[3:,4] = -eigvecs[...,0]
A[3:,5] = -eigvecs[...,1]

# Constants matrix (RHS of equation)
b = np.mat( np.zeros((6,1)) )
b[0:3,0] = -cam_inv[...,0:2] * point
b[3:,] = mean.T

res = np.linalg.solve(A,b)

（其中 cam_inv、eigvecs、mean 和 point 是在别​​处计算的其他一些矩阵。）

显然上面的代码可以有更多的注释，但我觉得即使有一些注释，它仍然无法真正传达正在解决的基本方程。有没有更好的方法将方程输入求解器，从而产生更清晰的代码？

4

1 回答 1

1

The problem is that lines of A that represent the equalities don't have a one-to-one mapping to lines of code. What I do in my own work (Economics) is to have a function with a clear English name (or at least a single line of code, with no functional representation) for each of the rows of A. When necessary, I have a clear but slow or perhaps much longer version of the code that does the same thing as the code I ultimately use.

So for example (from Bretscher's Linear Algebra with Applications 1997, ex. 37 p. 29 a simple if unrealistic example), consider an economy with three industries, I1, I2, I3, each taking the other two industry outputs as inputs. What outputs should they produce to meet consumer and industrial demand?

A =np.zeros((3,3)) 
#Each unit of production by I1 requires 0.1 units of good 2 and .2 of good 3
A[:,0] = [0, 0.1, 0.2] 
#Each unit of production by I2 requires 0.2 units of good 1 and .5 of good 3
A[:,1] = [0.2, 0, 0.5] 
#Each unit of production by I3 requires 0.3 units of good 1 and .4 of good 2    
A[:,2] = [0.3, 0.4, 0] 
#The required production for consumers. 
b = np.array([320,150,90]).reshape(-1,1)
#The optimal production levels of x1, x2, and x3
res = np.linalg.solve(A,b)

It perhaps would be slower or less terse to do what I suggest, but it would be much clearer to read.

于 2013-04-05T13:06:14.360 回答