问偏微分方程离散化的矩阵与循环
EN

Stack Overflow用户

提问于 2021-03-22 09:05:18

回答 2查看 56关注 0票数 1

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

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

python

julia

sparse-matrix

鹅厂大牛带你玩转AI+向量数据库

鹅厂证书、定制周边等你来拿

回答 2

Stack Overflow用户

回答已采纳

发布于 2021-03-27 11:20:48

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)

票数 0

Stack Overflow用户

发布于 2021-03-25 13:11:38

如果您没有进行就地修改，请使用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)