在scipy中用全矩阵模拟常微分方程（对象太深，无法满足数组要求）

Question

我有一个包含微分方程的系统，使用一个 m*m 的矩阵 S。矩阵 S 中的每个元素 S[i,j] 代表某种物种的浓度，这个浓度会受到 S[i-1,j] 和 S[i,j-1] 的影响。

我可以在每一步计算出每个元素的变化率 dx/dt（这个是通过 update_matrix 函数得到的），但接下来我需要对这些变化率进行积分，以更新我的初始浓度（x 和 x_counts 是同一个意思）。不过，scipy 的 integrate.odeint 函数不支持矩阵作为输入（或者返回值），而且会报“对象太深”的错误。

有没有办法让我调整一下，能够在这个矩阵上使用积分器呢？我提供了一段代码：diffeq 函数返回每个元素的 dx/dt，update_matrix 函数返回一个矩阵 B，其中 B[i,j] = dx[i,j]/dt。

def diffeq(x,i,j,coeffs):
    if (i==1 and j==0):
        return (-(coeffs.M(1,0)-coeffs.L(1,0)))+coeffs.d(1,0)-coeffs.get_phi()*x[1][0]
    elif j==0:
        return (-(coeffs.M(i,0)+coeffs.L(i,0)))*x[i][0]+coeffs.L(i-1,0)*x[i-1][0]+coeffs.d(i,j)-coeffs.get_phi()*x[i][0]
    elif (i>1 and j>0):
        return (-(coeffs.M(i,j)+coeffs.L(i,j)))*x[i][j]+coeffs.M(i,j-1)*x[i][j-1]+coeffs.L(i-1,j)*x[i-1][j]+coeffs.d(i,j)-coeffs.get_phi()*x[i][j]
    elif i==1 and j>0:
        return (-(coeffs.M(1,j)+coeffs.L(1,j)))*x[1][j]+coeffs.M(1,j-1)*x[1][j-1]+coeffs.d(1,j)-coeffs.get_phi()*x[1][j]
    elif i==0 and j==1:
        return -x[0][1]+coeffs.d(0,1)-coeffs.get_phi()*x[0][1]
    elif i==0 and j>1:
        return -j*x[0][j]+(j-1)*x[0][j-1]+coeffs.d(0,j)-coeffs.get_phi()*x[0][j]


def update_matrix(x,coeffs,m):
    update_matrix=numpy.zeros((m,m))
    for i in range(m+1):
        for j in range(m+1-i):
            update_matrix[m][m]=diffeq(x,i,j,coeffs)
    return update_matrix


def run_simulation_R2(a,q,m):

    x_counts=numpy.zeros((m,m))
    x_counts[1][0]=1
    x_counts[0][1]=1
    coeffs=R2(a,q,m,x_counts)
    t=range(0,100)
    output = integrate.odeint(update_matrix, x_counts, t, args=(coeffs, m))

Answer 1

如果 odeint 需要的是一个向量，而不是一个矩阵，那你就得给它一个向量。如果你不想大幅度改变你的代码逻辑，可以在函数外部把 x 设成一个形状为 (m**2,) 的向量，但在函数内部仍然保持它是一个 (m, m) 的矩阵。你可以在需要的地方灵活使用 .reshape(-1) 来实现这个转换。你没有提供足够的信息让我们完全测试，但类似下面的代码可能会有效：

def update_matrix(x,coeffs,m):
    x = x.reshape(m, m)
    update_matrix=numpy.zeros((m,m))
    for i in range(m+1):
        for j in range(m+1-i):
            update_matrix[m][m]=diffeq(x,i,j,coeffs)
    return update_matrix.reshape(-1)

def run_simulation_R2(a,q,m):
    x_counts=numpy.zeros((m,m))
    x_counts[1][0]=1
    x_counts[0][1]=1
    coeffs=R2(a,q,m,x_counts)
    t=range(0,100)
    output = integrate.odeint(update_matrix, x_counts.reshape(-1), t,
                              args=(coeffs, m))
    return output.reshape(m, m)