使用solve_ivp求解二维热方程

问题描述 投票:0回答:1

我正在尝试使用有限差分法作为使用solve_ivp 方法的图像过滤技术来求解二维热方程。我的代码显示了一条错误消息“

y0
必须是一维的”,我不明白为什么会这样。

# Read an image 
im3 = imageio.imread('photo.jpg')
#find out the dimensions of the image
height3, width3, colours3 = im3.shape

#convert the image to a numpy array for easier mathematical use
nparray3 = np.asarray(im3)
#convert the image to its grayscle version
grayimg3 = color.rgb2gray(nparray3)

#add noise on the gray image: The result is the initial 
#condition of our system
noisygrayimg3 =     
grayimg3+1.5*np.random.rand(grayimg3.shape[0],grayimg3.shape[1])
#Reshape the initial conditions into a vector so it can be     
#passed into the ODE solver
U0 = np.reshape(noisygrayimg3, (width3*height3, 1)) 

#Define a grid of x and y variables and their derivatives
#x variable
x = np.linspace(0,1,width3) #patial vector
dx = x[1]-x[0] #spatial step size
onex = np.ones(width3)
Ix = np.eye(width3)
Dx = sp.spdiags([onex , -2*onex , onex], [-1, 0, 1], width3, 
width3) #x 1D matrix 
#y variable
y = np.linspace(0,1,height3) #patial vector
dy = y[1]-y[0] #spatial step size
oney = np.ones(height3)
Iy = np.eye(height3)
Dy = sp.spdiags([oney , -2*oney , oney], [-1, 0, 1], height, 
height3) #y 1D matrix 

#create the 2D Laplacian
A = sp.kron(Dx,Iy) + sp.kron(Ix,Dy) 

#Define the Diffusion coefficient
D = 0.005

#solve the diffusion equation
#STEP 1: define the function that returns du/dt which computes the 
right-hand side of the heat equation
def Diffusion_rhs(t,u):    
    #caluclate the right-hand side
    dudt = D * A* u
    return dudt

#STEP 2: define  the Time variable for the integration
tspan = np.array([0, 2])

#STEP 3: solve the system
sol = solve_ivp(Diffusion_rhs, tspan, U0)

如有任何帮助,我们将不胜感激。 谢谢。

python ode pde heat
1个回答
0
投票

为了简化实现,scipy 求解器仅允许平面数组作为接口及其内部过程中的状态。如果您的状态是结构化的,那么您需要在 ODE 函数中应用适当的转换,从开始时的数组到结构,以及最后返回对象的导数从结构到数组。 numpy 对象存在

reshape
方法。

对于

odeint
方法,存在
odeintw
包,它根据初始状态的结构生成 ODE 函数的包装器。

其他软件包允许实现基本向量运算的结构化对象。一个例子是 Julia ODE 求解器包,它也可以通过 python 连接。

最新问题
© www.soinside.com 2019 - 2025. All rights reserved.