用于 PDE 离散化的矩阵与循环答案

【问题标题】：Matrix vs. loop for discretization of PDE用于 PDE 离散化的矩阵与循环
【发布时间】：2021-06-18 22:05:22
【问题描述】：

我正在尝试使用 diffeqpy 通过离散空间维度来求解 PDE，而我将时间维度视为一组常微分方程。我设法使用 for 循环解决了一个非常简单的问题。但是，当我尝试使用矩阵来扩大问题规模时，求解器提供的答案不正确。

以下代码有效：

from diffeqpy import de
import numpy as np

def f(du,u,p,t):
    #define shape of matrix
    s = (6,7)
    cc = np.matrix((np.zeros(s)))
              
    for j in range(0,6):
        for i in range(0,6):
            if (j == i):
                cc[j,i] = -1.0
                cc[j,i+1] = 1.0      
        
    for j in range(0,6):
        du[j] = cc[j,0]*u[0] + cc[j,1]*u[1] + cc[j,2]*u[2] + cc[j,3]*u[3] + cc[j,4]*u[4] + cc[j,5]*u[5] + cc[j,6]*u[6]

u0 = [0.1,0.0,0.0,0.0,0.0,0.0,1.0]

tspan = (0., 20.)
prob = de.ODEProblem(f, u0, tspan)
sol = de.solve(prob)

此代码类似于以下同样有效的代码：

from diffeqpy import de

def f(du,u,p,t):
    du[0] = -u[0]+u[1]
    du[1] = -u[1]+u[2]  
    du[2] = -u[2]+u[3]
    du[3] = -u[3]+u[4]
    du[4] = -u[4]+u[5]
    du[5] = -u[5]+u[6]
     
u0 = [0.1,0.0,0.0,0.0,0.0,0.0,1.0]

tspan = (0., 20.)
prob = de.ODEProblem(f, u0, tspan)
sol = de.solve(prob)

但是，当我尝试使用矩阵运算时，问题无法正确解决。我没有计算机科学背景。不过，我想了解更多。为什么以下代码不起作用？它与可变对象与不可变对象有关吗？如何利用矩阵将这个问题扩展到更大的离散化步骤？

from diffeqpy import de
import numpy as np

def f(du,u,p,t):
    #define shape of matrix
    
    s = (6,7)
    cc = np.matrix((np.zeros(s)))       
       
    for j in range(0,6):
        for i in range(0,6):
            if (j == i):
                cc[j,i] = -1.0
                cc[j,i+1] = 1.0     
   
    
    x = np.matrix(u).T
 
    du = (cc*x).T

u0 = [0.1,0.0,0.0,0.0,0.0,0.0,1.0]

tspan = (0., 20.)
prob = de.ODEProblem(f, u0, tspan)
sol = de.solve(prob)

我将不胜感激有关此问题的任何指导。

【问题讨论】：

不要只说“不起作用”。显示实际错误（带回溯）。
diffeqpy 是一个轻度使用（与numpy 相比）的包。它甚至没有标签。以及为什么使用sparse-matrix 标签。这些天不鼓励使用np.matrix；但我不知道这是否是你问题的一部分。如果我们不能运行这段代码（没有找到diffeqpy），我们也不能帮你调试。
代码没有错误。代码运行。但是，它不会改变问题的初始条件。因此，当我使用矩阵乘法时，“u”保持不变。
工作的f 不返回值，因此它们必须依赖于du 参数的就地修改。您的矩阵版本有一个du=...，它破坏了与参数du 的链接。

标签： python julia sparse-matrix

【解决方案1】：

如果您不进行就地修改，请使用 3 参数形式：

from diffeqpy import de
import numpy as np

def f(u,p,t):
    #define shape of matrix
    
    s = (6,7)
    cc = np.matrix((np.zeros(s)))       
       
    for j in range(0,6):
        for i in range(0,6):
            if (j == i):
                cc[j,i] = -1.0
                cc[j,i+1] = 1.0     
   
    
    x = np.matrix(u).T
 
    du = (cc*x).T

u0 = [0.1,0.0,0.0,0.0,0.0,0.0,1.0]

tspan = (0., 20.)
prob = de.ODEProblem(f, u0, tspan)
sol = de.solve(prob)

【讨论】：

谢谢你，克里斯，如果没有这个指导，我就无法让它工作。我只需要在函数底部添加一个“return du”，然后将 2D du 矩阵转换为 1D 矩阵。我已经发布了下面的代码。

【解决方案2】：

from diffeqpy import de
import numpy as np

def f(u,p,t):
    #define shape of matrix
    
    s = (6,6)
    cc = np.matrix((np.zeros(s)))       
       
    for j in range(0,6):
        for i in range(0,5):
            if (j == i):
                cc[j,i] = -1.0
                cc[j,i+1] = 1.0
    cc[-1,-1] = -1.0
                
    x = (np.matrix(u).T)
    
    du = (list(((cc*x).T).flat))
    
    return du

u0 = [0.1,0.0,0.0,0.0,0.0,0.1]

tspan = (0., 20.)
prob = de.ODEProblem(f, u0, tspan)
sol = de.solve(prob)

【讨论】：