热方程的FTCS算法

Question

我试图在Python中实现一维热方程的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，我没有得到u[t][x]=u（x，t）的合理值。i、 它们会爆炸，给我一些愚蠢的值，比如‘4.11052068e+221’，我参与的一些糟糕的编程实践会破坏算法吗？或者我只是错误地实现了算法？在

编辑：我认为这是因为算法是稳定的当且仅当r<；1/2。因为我的r大约是2.5或者其他的，所以数字会爆炸，但是如果有人看到任何其他错误，请告诉我！！在

Answer 1

原始问题

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

小的ish错误

1. dx = L/nx和{}。（如果您对此感到困惑，想象一下nx = 2和{}的情况。。。 然后dx = 0.5）。在
2. L和{}应该是浮动的。（即-L = 1.0）。在
3. 你没有更新棒上倒数第二个位置的温度u。确保for循环覆盖位置扩展到nx-1，而不是nx-2。在

需要考虑的事项

1. 您可以考虑使用时间、长度、热导率和初始温度的实际值/单位。

2. 当您应用边界条件时（在注释BC下），只有时间t=0的值为u，这些值是非零的，因此您不需要迭代所有的时间步。

3. 虽然有必要遍历时间（因此产生了timefor循环），但实际上可以将空间导数矢量化。可以将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])