使用solve_ivp求解2D热方程做图像滤波时遇维度错误求助

错误原因 #

solve_ivp要求初始条件y0必须是一维NumPy数组，你当前通过np.reshape(noisygrayimg3, (width3*height3, 1))得到的是二维列向量（形状为(N,1)），不符合函数输入要求，因此触发报错。

另外代码还有两处细节错误需要修正：

• 构造Dy矩阵时，参数里的height是未定义变量，应该改为height3
• 有限差分的二阶导数矩阵需要除以空间步长的平方（dx²和dy²），否则拉普拉斯算子的尺度不正确，会导致扩散效果偏离预期

修改方案 #

1. 将初始条件U0转换为一维数组，可通过以下任意一种方式实现：
   • U0 = noisygrayimg3.flatten()
   • U0 = np.reshape(noisygrayimg3, width3*height3)
   • U0 = noisygrayimg3.ravel()
2. 修正Dy矩阵的构造参数，将未定义的height替换为height3
3. 给Dx和Dy分别除以dx**2和dy**2，还原正确的拉普拉斯算子尺度

修改后的完整代码 #

import imageio
import numpy as np
from skimage import color
import scipy.sparse as sp
from scipy.integrate import solve_ivp

# 读取图像
im3 = imageio.imread('photo.jpg')
# 获取图像尺寸
height3, width3, colours3 = im3.shape

# 转换为灰度图
grayimg3 = color.rgb2gray(im3)

# 添加噪声作为初始条件
noisygrayimg3 = grayimg3 + 1.5 * np.random.rand(*grayimg3.shape)
# 将初始条件转为一维数组（符合solve_ivp要求）
U0 = noisygrayimg3.flatten()

# 构造x方向的二阶差分矩阵
x = np.linspace(0, 1, width3)
dx = x[1] - x[0]
onex = np.ones(width3)
Ix = np.eye(width3)
# 除以dx²，还原正确的二阶导数尺度
Dx = sp.spdiags([onex, -2*onex, onex], [-1, 0, 1], width3, width3) / (dx**2)

# 构造y方向的二阶差分矩阵
y = np.linspace(0, 1, height3)
dy = y[1] - y[0]
oney = np.ones(height3)
Iy = np.eye(height3)
# 修正height为height3，同时除以dy²
Dy = sp.spdiags([oney, -2*oney, oney], [-1, 0, 1], height3, height3) / (dy**2)

# 构造二维拉普拉斯算子矩阵
A = sp.kron(Dx, Iy) + sp.kron(Ix, Dy)

# 扩散系数
D = 0.005

# 定义热方程右端项
def Diffusion_rhs(t, u):
    dudt = D * A @ u  # 使用矩阵乘法运算符@更清晰
    return dudt

# 时间区间
tspan = [0, 2]

# 求解ODE
sol = solve_ivp(Diffusion_rhs, tspan, U0)

# 可选：将最终结果转换回图像形状
filtered_img = sol.y[:, -1].reshape(height3, width3)

内容的提问来源于stack exchange，提问作者User123