如何在Python中建立和求解联立方程

对于固定整数n，我有一组2(n-1)联立方程如下。

M(p) = 1+((n-p-1)/n)*M(n-1) + (2/n)*N(p-1) + ((p-1)/n)*M(p-1)

N(p) = 1+((n-p-1)/n)*M(n-1) + (p/n)*N(p-1)

M(1) = 1+((n-2)/n)*M(n-1) + (2/n)*N(0)

N(0) = 1+((n-1)/n)*M(n-1)

M(p)定义为1 <= p <= n-1. N(p)定义为0 <= p <= n-2。还请注意p只是每个方程中的常整数，因此整个系统是线性的。

我一直在使用 Maple，但我现在想设置这些并在 python 中解决它们，也许使用numpy.linalg.solve（或任何其他更好的方法）。我实际上只想要的值M(n-1)。例如，当n=2答案应该是M(1) = 4， 我相信。这是因为方程变为

M(1) = 1+(2/2)*N(0)
N(0) = 1 + (1/2)*M(1)

所以

M(1)/2 = 1+1

and so

M(1) = 4. 

如果我想插入n=50比如说，如何在 python 中建立联立方程组以便 numpy.linalg.solve 可以求解它们？

Update答案很好，但他们使用密集的求解器，而方程组却稀疏。已发布跟进使用 scipy 稀疏矩阵求解方程组 https://stackoverflow.com/questions/14380899/using-scipy-sparse-matrices-to-solve-system-of-equations .

更新：使用 scipy.sparse 添加实现

这给出了按顺序的解决方案N_max,...,N_0,M_max,...,M_1.

要求解的线性系统的形状为A dot x == const 1-vector. x是寻求的解向量。
在这里我对方程进行了排序，使得x is N_max,...,N_0,M_max,...,M_1.
然后我建立了A-4 个分块矩阵的系数矩阵。

Here's a snapshot for the example case n=50 showing how you can derive the coefficient matrix and understand the block structure. The coefficient matrix A is light blue, the constant right side is orange. The sought after solution vector x is here light green and used to label the columns. The first column show from which of the above given eqs. the row (= eq.) has been derived:

正如 Jaime 建议的那样，乘以n改进代码。这没有反映在上面的电子表格中，但已在下面的代码中实现：

使用numpy实现：

import numpy as np
import numpy.linalg as linalg


def solve(n):
    # upper left block
    n_to_M = -2. * np.eye(n-1) 

    # lower left block
    n_to_N = (n * np.eye(n-1)) - np.diag(np.arange(n-2, 0, -1), 1)

    # upper right block
    m_to_M = n_to_N.copy()
    m_to_M[1:, 0] = -np.arange(1, n-1)

    # lower right block
    m_to_N = np.zeros((n-1, n-1))
    m_to_N[:,0] = -np.arange(1,n)

    # build A, combine all blocks
    coeff_mat = np.hstack(
                          (np.vstack((n_to_M, n_to_N)),
                           np.vstack((m_to_M, m_to_N))))

    # const vector, right side of eq.
    const = n * np.ones((2 * (n-1),1))

    return linalg.solve(coeff_mat, const)

使用 scipy.sparse 的解决方案：

from scipy.sparse import spdiags, lil_matrix, vstack, hstack
from scipy.sparse.linalg import spsolve
import numpy as np


def solve(n):
    nrange = np.arange(n)
    diag = np.ones(n-1)

    # upper left block
    n_to_M = spdiags(-2. * diag, 0, n-1, n-1)

    # lower left block
    n_to_N = spdiags([n * diag, -nrange[-1:0:-1]], [0, 1], n-1, n-1)

    # upper right block
    m_to_M = lil_matrix(n_to_N)
    m_to_M[1:, 0] = -nrange[1:-1].reshape((n-2, 1))

    # lower right block
    m_to_N = lil_matrix((n-1, n-1))
    m_to_N[:, 0] = -nrange[1:].reshape((n-1, 1))

    # build A, combine all blocks
    coeff_mat = hstack(
                       (vstack((n_to_M, n_to_N)),
                        vstack((m_to_M, m_to_N))))

    # const vector, right side of eq.
    const = n * np.ones((2 * (n-1),1))

    return spsolve(coeff_mat.tocsr(), const).reshape((-1,1))

示例n=4:

[[ 7.25      ]
 [ 7.76315789]
 [ 8.10526316]
 [ 9.47368421]   # <<< your result
 [ 9.69736842]
 [ 9.78947368]]

示例n=10:

[[ 24.778976  ]
 [ 25.85117842]
 [ 26.65015984]
 [ 27.26010007]
 [ 27.73593401]
 [ 28.11441922]
 [ 28.42073207]
 [ 28.67249606]
 [ 28.88229939]
 [ 30.98033266]  # <<< your result
 [ 31.28067182]
 [ 31.44628982]
 [ 31.53365219]
 [ 31.57506477]
 [ 31.58936225]
 [ 31.58770694]
 [ 31.57680467]
 [ 31.560726  ]]