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

Question

对于固定整数 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

等等

M(1) = 4. 

如果我要插n=50 ，比如说，你如何在 python 中建立这个联立方程组，以便 numpy.linalg.solve 可以解决它们？

更新 答案很好，但他们在方程组稀疏的地方使用密集求解器。发布跟进 Using scipy sparse matrices to solve system of equations .

Answer 1

更新:添加了使用 scipy.sparse 的实现

---

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

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

下面是示例案例 n=50 的快照，展示了如何推导系数矩阵和理解 block 结构。系数矩阵A为淡蓝色，常量右侧为橙色。寻求解决方案向量 x 在这里是浅绿色的，用于标记列。第一列显示上述给定方程中的哪一个。行 (= eq.) 已导出:

正如 Jaime 所建议的，乘以 n 可以改进代码。这没有反射(reflect)在上面的电子表格中，但已在以下代码中实现:

使用 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  ]]