二维热传导PDE有限差分解法实现及边界响应异常问题

问题概述 #

已实现一维偏微分方程（PDE）的有限差分解法，扩展至二维热传导方程时，先遭遇IndexError: too many indices for array索引错误；调整代码改用odeint求解后，设置恒定均匀输入、初始与边界条件均为0的情况下，仅单一边界出现温度响应，其余区域无变化。

相关代码 #

一维实现代码 #

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
# number of discretization points
n = 100
# length of bar (1D object/surface)
L = 0.1
# make a vector containing location of each discretization point 
X = np.linspace(0, L, n)

h = L/(n - 1)

kappa = 2*10**(-9)
def odefunc(T, t):
    dTdt = np.zeros(X.shape)
    dTdt[0] = 0
    dTdt[-1] = 0
    for i in range(1, n-1):
        # 注：此处一维热传导方程多除了一个h²，正确应为 kappa/h²*(T[i+1]-2*T[i]+T[i-1])
        dTdt[i] = (kappa/h**2) * (T[i+1]- 2*T[i] + T[i -1])/h**2
    return dTdt
init = 0*np.ones(X.shape)
init[0] = 0.01
init[-1] = 0
tt = np.linspace(0.0, 1.2, 100).round(3)
sol = odeint(odefunc, init, tt)
sol[0] = 0

二维初始错误代码 #

# length of surface 
L = 0.1
# height of surface
ht = 0.1

nx,ny = (100,100)
x = np.linspace(0,L,nx)
y = np.linspace(0,ht,ny)
# make a grid containing location of each discretization point 
X,Y = np.meshgrid (x,y)

def odefunc(T, t):
 dTdt = np.zeros_like(X)

 dTdt[0][0] = 0
 dTdt[0][-1] = 0
 dTdt[-1][0] = 0
 dTdt[-1][-1] = 0
 Txx = np.zeros((100,100))
 Tyy = np.zeros((100,100))
 for i in range(1, nx-1):
     for j in range(1,ny-1):
         Txx = (T[i+1,j] - 2*T[i,j] + T[i-1,j])
         Tyy = (T[i,j+1]-2*T[i,j]+T[i,j-1])
         dTdt[i][j] = Txx + Tyy
 return dTdt

修正后代码（存在关键错误） #

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.integrate import odeint
# number of discretization points
n = 100
# length of surface (1D object/surface)
L = 0.1
# height of surface
ht = 0.1
#
nx,ny = (100,100)
dTdt = np.zeros((nx,ny))
# make a matrix containing location of each discretization point 
dx = L/nx
dy = ht/ny


def odefunc(T, t):
    magnitude = 2
    selection = 'uniform'
    nx,ny = (100,100)


    dTdt = np.zeros((nx,ny))
    # 边界条件被注释
    #dTdt[0,0] = 0
    #dTdt[0,-1] = 0
    #dTdt[-1,0] = 0
    #dTdt[-1,-1] = 0
    # 关键错误：将传入的状态T重置为全零数组
    T = np.zeros((100,100))
    input = np.zeros((100,100))
    for i in range(1, nx-1):
        for j in range(1,ny-1):
            Txx = (T[i+1,j] - 2*T[i,j] + T[i-1,j])/(dx**2)
            Tyy = (T[i,j+1]-2*T[i,j]+T[i,j-1])/(dy**2)
            input = 2
            dTdt[i,j] = Txx + Tyy +input
    result = dTdt.flatten() 

    return result

init = np.zeros((nx,ny))
inituse = init.flatten()

tt = np.linspace(0.0, 1.2, 100).round(3)
sol = odeint(odefunc, inituse, tt)

# reshape顺序错误，应为(time_steps, nx, ny)
discretesol = sol.reshape(100,100,100)

for k in range(0, len(tt), 2):
        x = np.linspace(0,L,nx)
        y = np.linspace(0,ht,ny)
        X,Y = np.meshgrid(x,y)
        # 索引顺序错误，应取discretesol[k,:,:]
        Z = discretesol[:,:,k]
        
        fig = plt.figure()
        ax = fig.gca(projection='3d')
        surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1)
        plt.show()

异常原因分析 #

1. 核心状态覆盖错误：修正后代码的odefunc中，将odeint传入的当前温度状态数组T直接重置为全零数组，导致所有计算都基于无演化的全零温度场，完全忽略了时间步的状态传递。
2. 边界条件未正确处理：代码注释了边界点的dTdt设置，且循环仅处理内部点，边界点的温度变化率始终为0，未匹配热传导方程的物理边界规则。
3. 状态数组reshape与切片错误：sol的形状为(时间步数, 总离散点数)，每个行对应一个时间步的扁平温度场，原代码reshape顺序错误，后续切片取的不是对应时间步的温度分布。
4. 热扩散系数遗漏：二维热传导方程的标准形式包含热扩散系数，原代码遗漏该参数，导致温度演化速率不符合物理规律。

修复方案 #

关键修改点 #

1. 保留并恢复传入的状态数组：删除T = np.zeros((100,100))，将扁平的输入状态转为二维数组：T = T.reshape(nx, ny)。
2. 明确边界条件：以恒温边界（温度固定为0）为例，设置边界点的温度变化率为0。
3. 修正reshape与切片逻辑：将solreshape为(时间步数, nx, ny)，取时间步时使用discretesol[k,:,:]。
4. 补充热扩散系数：按照热传导方程标准形式加入热扩散系数kappa，保证物理逻辑正确。

修复后完整代码 #

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# 离散参数
nx, ny = 100, 100
L, ht = 0.1, 0.1
dx = L / nx
dy = ht / ny
# 热扩散系数（示例值，可根据实际调整）
kappa = 2e-9
# 均匀热源强度
q = 2

def odefunc(T_flat, t):
    # 将扁平的状态数组转为二维温度场
    T = T_flat.reshape(nx, ny)
    dTdt = np.zeros_like(T)
    
    # 处理内部点的热传导方程
    for i in range(1, nx-1):
        for j in range(1, ny-1):
            Txx = (T[i+1,j] - 2*T[i,j] + T[i-1,j]) / dx**2
            Tyy = (T[i,j+1] - 2*T[i,j] + T[i,j-1]) / dy**2
            dTdt[i,j] = kappa * (Txx + Tyy) + q
    
    # Dirichlet边界条件：边界温度固定为0，温度变化率为0
    dTdt[0, :] = 0
    dTdt[-1, :] = 0
    dTdt[:, 0] = 0
    dTdt[:, -1] = 0
    
    # 转为扁平数组返回
    return dTdt.flatten()

# 初始条件：全区域温度为0
init = np.zeros((nx, ny))
init_flat = init.flatten()

# 时间点
tt = np.linspace(0.0, 1.2, 100)
# 求解ODE
sol = odeint(odefunc, init_flat, tt)

# reshape为(时间步, nx, ny)
discretesol = sol.reshape(len(tt), nx, ny)

# 可视化部分
for k in range(0, len(tt), 10):
    x = np.linspace(0, L, nx)
    y = np.linspace(0, ht, ny)
    X, Y = np.meshgrid(x, y)
    Z = discretesol[k, :, :]
    
    fig = plt.figure(figsize=(8,6))
    ax = fig.add_subplot(projection='3d')
    surf = ax.plot_surface(X, Y, Z, cmap='viridis', rstride=1, cstride=1)
    ax.set_xlabel('X')
    ax.set_ylabel('Y')
    ax.set_zlabel('Temperature')
    ax.set_title(f'Temperature Field at t={tt[k]:.2f}')
    fig.colorbar(surf)
    plt.show()

说明 #

修复后的代码会正确模拟二维热传导过程：在均匀热源作用下，内部区域温度逐渐升高，边界因恒温条件保持为0，温度场从边界向中心逐步上升，符合热传导的物理规律。

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