如何在2D热方程求解中设置对角线热源及迭代终止条件？

一、定义对角线/曲线形态的热源 #

1. 核心思路 #

替代原代码中硬编码的矩形/直线热源，通过自定义判断函数识别每个网格点(i,j)是否属于目标热源区域，可灵活实现对角线、圆形、抛物线等任意形态的热源。

2. 具体实现 #

（1）对角线热源示例

定义一条从(10,10)到(40,40)的3像素宽粗对角线：

def is_diagonal_source(i, j):
    # abs(i-j)<=2控制对角线厚度，10<=i<=40限制对角线范围
    return abs(i - j) <= 2 and 10 <= i <= 40 and 10 <= j <= 40

（2）曲线热源示例（圆形）

定义圆心在(25,25)、半径10的圆形热源：

def is_circle_source(i, j):
    center_x, center_y = 25, 25
    radius = 10
    # 利用圆的方程判断点是否在圆内
    return (i - center_x)**2 + (j - center_y)**2 <= radius**2

（3）将热源应用到模拟中

替换原代码的边界条件设置逻辑，遍历所有网格点标记热源：

# 初始化初始温度场
u_initial = 0.0
source_temp = 100.0
u0 = np.full((plate_length, plate_length), u_initial)

# 标记热源区域
for i in range(plate_length):
    for j in range(plate_length):
        # 可同时启用多个热源
        if is_diagonal_source(i, j) or is_circle_source(i, j):
            u0[i, j] = source_temp

同时在迭代计算时，跳过热源点的温度更新（保持热源温度恒定）：

if is_diagonal_source(i,j) or is_circle_source(i,j):
    continue

二、修改迭代逻辑：温度达到最大值时自动停止 #

1. 核心思路 #

放弃固定迭代次数，改为动态循环：每次迭代后检查平板区域的最高温度是否接近热源温度（如100.0），达到阈值则停止；同时设置最大迭代次数作为兜底，避免无限循环。

2. 具体实现 #

（1）动态存储温度场

不再预分配固定大小的数组，改用列表存储每一步的温度场：

# 初始化温度场列表
u = []
u.append(u0)  # 添加初始状态

max_temp_target = source_temp
tolerance = 1e-3  # 浮点精度容忍度
max_iter_fallback = 500  # 最大迭代次数兜底
current_iter = 0

（2）迭代循环与停止条件

while current_iter < max_iter_fallback:
    u_prev = u[-1]
    u_next = u_prev.copy()
    
    # 更新非热源区域的温度
    for i in range(1, plate_length-1):
        for j in range(1, plate_length-1):
            if is_diagonal_source(i,j) or is_circle_source(i,j):
                continue
            u_next[i,j] = gamma * (u_prev[i+1][j] + u_prev[i-1][j] + u_prev[i][j+1] + u_prev[i][j-1] - 4*u_prev[i][j]) + u_prev[i][j]
    
    u.append(u_next)
    current_iter += 1
    
    # 检查停止条件
    current_max_temp = np.max(u_next)
    if abs(current_max_temp - max_temp_target) < tolerance:
        print(f"模拟自动停止：迭代次数={current_iter}，当前最高温度={current_max_temp:.3f}")
        break

（3）更新动画生成逻辑

动画帧数改为实际迭代次数（即列表u的长度）：

def animate(k):
    plotheatmap(u[k], k)

anim = animation.FuncAnimation(plt.figure(), animate, interval=1, frames=len(u), repeat=False)
anim.save("heat_equation_solution.gif")

完整修改后的代码 #

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

print("2D heat equation solver")

plate_length = 50
alpha = 2
delta_x = 1
delta_t = (delta_x ** 2)/(4 * alpha)
gamma = (alpha * delta_t) / (delta_x ** 2)

# 定义热源判断函数
def is_diagonal_source(i, j):
    return abs(i - j) <= 2 and 10 <= i <= 40 and 10 <= j <= 40

def is_circle_source(i, j):
    center_x, center_y = 25, 25
    radius = 10
    return (i - center_x)**2 + (j - center_y)**2 <= radius**2

# 初始化初始温度场
u_initial = 0.0
source_temp = 100.0
u0 = np.full((plate_length, plate_length), u_initial)

# 设置热源区域
for i in range(plate_length):
    for j in range(plate_length):
        if is_diagonal_source(i, j):
            u0[i, j] = source_temp

# 动态存储温度场
u = []
u.append(u0)

max_temp_target = source_temp
tolerance = 1e-3
max_iter_fallback = 500
current_iter = 0

# 迭代计算
while current_iter < max_iter_fallback:
    u_prev = u[-1]
    u_next = u_prev.copy()
    
    for i in range(1, plate_length-1):
        for j in range(1, plate_length-1):
            if is_diagonal_source(i,j):
                continue
            u_next[i,j] = gamma * (u_prev[i+1][j] + u_prev[i-1][j] + u_prev[i][j+1] + u_prev[i][j-1] - 4*u_prev[i][j]) + u_prev[i][j]
    
    u.append(u_next)
    current_iter += 1
    
    current_max = np.max(u_next)
    if abs(current_max - max_temp_target) < tolerance:
        print(f"Stopped at iteration {current_iter}, max temp: {current_max:.3f}")
        break

# 绘图函数
def plotheatmap(u_k, k):
    plt.clf()
    plt.title(f"Temperature at t = {k*delta_t:.3f} unit time")
    plt.xlabel("x")
    plt.ylabel("y")
    plt.pcolormesh(u_k, cmap=plt.cm.jet, vmin=0, vmax=100)
    plt.colorbar()
    return plt

# 生成动画
def animate(k):
    plotheatmap(u[k], k)

anim = animation.FuncAnimation(plt.figure(), animate, interval=1, frames=len(u), repeat=False)
anim.save("heat_equation_solution.gif")

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