python之二维传热之溢出错误

我尝试使用以下公式使用 python 模拟铝中的二维热传递:

dT/dt = K*(d^2T/d^2x + d^2T/d^2y)

来源: https://www.youtube.com/watch?v=V00p-TgQML0

Python 代码使用 dx = dy = 1 (mm) 但如果 dx 和 dy 变小，我会收到溢出错误，我不知道如何避免。

import numpy as np 
import matplotlib.pyplot as plt 
from copy import deepcopy 
from matplotlib.animation import FuncAnimation 
import time 
 
 
x = 11 
y = 11 
 
sd = 1 
nx = x*sd 
ny = y*sd 
 
dx = 1/float(sd) 
dy = dx 
 
#Initial Temperature 
T0 = 25 
 
# Source: https://en.wikipedia.org/wiki/Thermal_diffusivity   
# Aluminium Thermal diffusivity (mm**2/s)  
K = 97 
 
#Time 
t0 = 0 
te = 1 
st = 1000 
dt = 1/float(st) 
 
 
#Iterations 
N = (te - t0)*st 
 
T = np.zeros(shape=(nx, ny)) 
 
T[:,:] = T0 
 
# Dirichlet Condition 
T[nx/2,ny/2] = 1000 
 
 
MM = [] 
 
for n in range(N): 
 
    Tb = deepcopy(T) 
 
    for i in range(nx): 
        for j in range(ny): 
 
            #Neumann Boundary Conditions 
            #TOP 
            if i == 0: 
                T[i,j] = Tb[i+1,j] 
 
            #RIGHT 
            elif j == ny1-1 and i != 0 and i != nx1-1:  
                T[i,j] = Tb[i,j-1] 
 
            #BOTTOM 
            elif i == nx-1:   
                T[i,j] = Tb[i-1,j] 
 
            #LEFT 
            elif j==0 and i != 0 and i != nx-1: 
                T[i,j] = Tb[i,j+1] 
 
            else: 
 
                T[i,j] = Tb[i,j] + K*(dt/dx**2)*(Tb[i+1,j]+Tb[i-1,j]+Tb[i,j+1]+Tb[i,j-1]-4*Tb[i,j]) 
            T[nx/2,ny/2] = 200 
 
    MM.append(T.copy()) 
 
fig = plt.figure() 
pcm = plt.pcolormesh(MM[0]) 
plt.colorbar() 
 
# Function called to update the graphic 
def step(i): 
    if i >= len(MM): return 0 
    pcm.set_array(MM[i].ravel()) 
    plt.title("Time: {0} s\n".format(i*dt)) 
    plt.draw() 
 
anim = FuncAnimation(fig, step, interval=3) 
plt.show() 

为了重现溢出错误，我将 sd 的值更改为 10(每个 mm 将分为 10 个元素)。

请您参考如下方法:

基本上 T[i,j] 发散并达到非常高的 + 和 - 值，并且在某个点达到 double_scale 类型的极限。

您可以通过在 T[i,j] = Tb[i,j] + K*(dt/dx* *2)*(Tb[i+1,j]+Tb[i-1,j]+Tb[i,j+1]+Tb[i,j-1]-4*Tb[i,j]) 在 sd=10 和 st=1000 的情况下。

这不是 python 或 numpy 问题，而是在尝试使用您正在使用的方法对刚性微分方程进行数值求解时出现的数值数学问题。

当您决定以更高的空间分辨率求解时，您也应该以更高的时间分辨率求解。我已经测试了代码，它适用于:sd=2、st=5000 和 sd=4、st= 10000。您会看到这种模式。

或者为您的微分方程使用更好的数值解。就像向后微分公式 (BDF) 一样，您可以在其中采取更大的时间步长而不会导致数值求解器发散。在这里寻找灵感:

https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.integrate.ode.html