PDE 离散化的矩阵与循环

Question

我正在尝试使用 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)

如果能就此问题提供任何指导，我将不胜感激。

Answer 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)

Answer 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)

PDE 离散化的矩阵与循环

Matrix vs. loop for discretization of PDE

python

julia

sparse-matrix