我正在尝试使用有限差分法作为使用solve_ivp 方法的图像过滤技术来求解二维热方程。我的代码显示了一条错误消息“
y0
必须是一维的”，我不明白为什么会这样。
import numpy as np
import imageio.v2 as imageio
from skimage import color
from skimage import io
import scipy.sparse as sp
from scipy.integrate import solve_ivp
im3 = imageio.imread('photo.jpg')
#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))
#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):
# Read an image
im3 = imageio.imread('photo.jpg')
#find out the dimsnions of the image
height3, width3, colours3 = im3.shape
#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
500x500
#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], height3, height3) #y 1D matrix
400x400
#create the 2D Laplacian
A = sp.kron(Dx,Iy) + sp.kron(Ix,Dy)
#Define the Diffusion coefficient
D = 0.005
#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)
任何帮助将不胜感激。 谢谢。
错误信息“
y0
必须是一维的”来自
solve_ivp
，因为它被设计为处理一阶常微分方程（ODE）系统，其中未知函数是一维向量。在你的情况下，你试图求解一个偏微分方程（PDE），即二维热方程。
虽然你可以使用
solve_ivp
求解 PDE，但你需要将二维空间问题离散化为一维 ODE 系统。这可以使用有限差分法（你已经提到的）完成。这个想法是用有限差分逼近空间导数，这会产生一个关于每个网格点上未知函数值的 ODE 系统。
下面是如何使用
solve_ivp
通过有限差分法求解二维热方程的更详细的细分：
- 空间离散化：
-
创建一个网格，覆盖你的图像域。你已经用
x和y完成了此操作。 -
用中心差分逼近拉普拉斯算子：
u_xx[i, j] ≈ (u[i+1, j] - 2*u[i, j] + u[i-1, j]) / dx**2 u_yy[i, j] ≈ (u[i, j+1] - 2*u[i, j] + u[i, j-1]) / dy**2其中u[i, j]表示网格点(x[i], y[j])处的温度。 -
ODE 系统：
-
将每个网格点上的偏微分方程转换为一个常微分方程：
du[i, j]/dt = D * (u_xx[i, j] + u_yy[i, j]) -
这将为你提供
width3 * height3个 ODE 的一个系统。 -
向量化：
-
将二维数组
u[i, j]转换为solve_ivp可以处理的一维向量u_vec。你可以通过堆叠列来做到这一点。 -
调整拉普拉斯逼近，以便它在
u_vec上运行。这可以使用矩阵乘法来完成，其中矩阵表示了有限差分方案。 -
求解 ODE 系统：
-
定义一个函数
diffusion_rhs(t, u_vec)，该函数计算 ODE 系统的右侧。 -
使用
solve_ivp求解系统，并将初始条件作为展平的噪声图像传递。 -
重塑：
-
solve_ivp返回的解也将是一维向量。你需要将其重塑回二维数组，以获得每个网格点和时间步长的温度。
以下是用
solve_ivp
求解二维热方程的 Python 代码：
import numpy as np
import imageio.v2 as imageio
from skimage import color
from scipy.integrate import solve_ivp
# Load and preprocess the image
im3 = imageio.imread('photo.jpg')
grayimg3 = color.rgb2gray(im3)
noisygrayimg3 = grayimg3 + 1.5 * np.random.rand(grayimg3.shape[0], grayimg3.shape[1])
# Spatial discretization
height3, width3 = grayimg3.shape
x = np.linspace(0, 1, width3)
y = np.linspace(0, 1, height3)
dx = x[1] - x[0]
dy = y[1] - y[0]
# Construct the Laplacian matrix
D = 0.005 # Diffusion coefficient
A = (np.eye(width3 * height3, k=-1) + np.eye(width3 * height3, k=1) - 4 * np.eye(width3 * height3)) / dx**2
A += (np.eye(width3 * height3, k=-width3) + np.eye(width3 * height3, k=width3)) / dy**2
A = D * A
# Define the RHS of the ODE system
def diffusion_rhs(t, u_vec):
return A.dot(u_vec)
# Set up initial condition and time span
u0_vec = noisygrayimg3.flatten()
tspan = (0, 2)
# Solve the ODE system
sol = solve_ivp(diffusion_rhs, tspan, u0_vec, t_eval=np.linspace(tspan[0], tspan[1], 50))
# Reshape the solution back to a 2D array
u_sol = np.reshape(sol.y, (height3, width3, -1))
# Now u_sol contains the temperature at each grid point and time step
这段代码首先创建了表示二维拉普拉斯算子的稀疏矩阵
A
。然后，它定义了一个函数
diffusion_rhs
，该函数将展平的温度向量作为输入，并返回相同长度的导数向量。最后，它使用
solve_ivp
求解 ODE 系统，并将初始条件作为展平的噪声图像传递。该解决方案然后被重塑回二维数组，以获得每个网格点和时间步长的温度。