热方程的FTCS算法

Question

我试图用Python实现1维热方程的FTCS算法。

import numpy as np

L = 1 #Length of rod in x direction
k = 0.3 #Thermal conductivity of rod
tmax = 5 #how many seconds
nx = 100 #number of spacial steps
nt = 100; #number of time steps
xi = np.linspace(0,L,nx) 
ti = np.linspace(0,tmax,nt)
dx = L/(nx-1)
dt = tmax/(nt-1)
r = k*dt/(dx)**2
r2 = 1-2*r

u=np.zeros((nt,nx))
#IC
phi = 100;
for x in range(0,nx):    
    u[0][x] = phi

#BC
for t in range(0,nt):    
    u[t][0] = 0;
    u[t][nx-1] = 0
#FTCS Algorithm

for t in range(0,nt-1): #timestep
    for x in range(1,nx-2):
        u[t+1][x] = r*(u[t][x-1]+ u[t][x+1]) + r2*(u[t][x])

然而，考虑到我的初始条件为100，对于ut = u(x，t)，我并没有得到可信的值。也就是说，它们会爆炸，给我一些愚蠢的值，比如'4.11052068e+221‘--我参与的一些糟糕的编程实践正在破坏算法吗？还是我刚刚实现的算法不正确？

编辑:我认为这是因为算法是稳定的当且仅当r< 1/2。因为我的r是2.5左右或者其他什么，那么这个数字就会爆炸，但是如果有人看到任何其他错误，请告诉我！！

Answer 1

原始问题

FTCS格式的稳定性取决于常数r的大小。如果是r<1/2，那么每一步引入的舍入误差都会呈指数衰减。如果是r>1/2，那么这些舍入误差就会成倍增加。(正如你在编辑中提到的那样)。

小错误

1. dx = L/nx和dt = tmax/nt.(如果这让你感到困惑，想象一下nx = 2和L = 1的情况.然后是dx = 0.5)。
2. L和tmax应该是浮动的。(ie - L = 1.0)。
3. 在杆上第二到最后的位置，你不能更新温度u。确保for循环超过位置扩展到nx-1，而不是nx-2。

要考虑的事情

1. 您可以考虑使用实数/单位的时间，长度，导热系数和初始温度。
2. 当您应用边界条件(在注释#BC下)时，只有时间t=0的u值为非零，因此不需要迭代所有时间步骤。
3. 虽然迭代时间(因此是time for循环)是必要的，但实际上您可以将空间导数向量化。可以用以下内容替换for循环：for t in range(0, nt-1): u[t+1, 1:nx-1] = r*(u[t,0:nx-2] + u[t,2:]) + r2*(u[t,1:nx-1])