如何在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可以解决它们？

更新 答案很好,但他们使用密集求解器,方程组很稀疏.发表后续使用scipy稀疏矩阵求解方程组.

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.
然后我A从4个块矩阵构建了- 系数矩阵.

这是示例案例的快照,n=50展示了如何推导系数矩阵并理解块结构.系数矩阵A为浅蓝色,恒定右侧为橙色.寻求的解决方案矢量x在这里是浅绿色并用于标记列.第一列显示了上面给出的eqs中的哪一个.行(= eq.)已派生:

正如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  ]]