热方程的 FTCS 算法答案

【问题标题】：FTCS Algorithm for the heat equation热方程的 FTCS 算法
【发布时间】：2017-11-06 02:56:59
【问题描述】：

我正在尝试在 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) 的合理值。也就是说，它们爆炸并给了我像 '4.11052068e+221' 这样的愚蠢值我参与的一些糟糕的编程实践正在破坏算法？还是我刚刚错误地实现了算法？

编辑：我想我已经得出这是因为当且仅当 r

【问题讨论】：

标签： python physics pde

【解决方案1】：

原始问题

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

小错误

dx = L/nx 和 dt = tmax/nt。（如果这让您感到困惑，想象一下nx = 2 和L = 1 的情况... 然后dx = 0.5)。
L 和 tmax 应该是浮点数。（即 - L = 1.0）。
您没有更新杆上倒数第二个位置的温度u。确保 for 循环位置延伸到 nx-1 而不是 nx-2。

需要考虑的事项

您可以考虑对时间、长度、热导率和初始温度使用实际值/单位。
当您应用边界条件时（在注释 #BC 下），只有时间 t=0 的值是非零的 u，因此您不需要遍历所有时间步.
虽然您需要遍历时间（因此是时间 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])

【讨论】：