使用 Python 的 FFT 和 IFFT 编码热方程解的问题

问题描述

我正在尝试使用 Python 的快速傅立叶变换 (FFT) 来实现热方程（一维）的解决方案。目标是一种有效的数值方案，用于说明在半无限域边界处指定的与时间相关的温度。这在传热文献中被称为 Duhamel 问题（参见，例如 Carslaw 和 Jaeger，1959，第 2.5 节）。对于指定为时间函数的边界处的强迫，已知精确解。在本案例中，我寻求对作为定期收集的离散数据给出的边界历史的响应。

对热方程进行傅里叶变换和适当的边界条件 ( ![equation]https://latex.codecogs.com/gif.download?T%280%2Ct%29%20%3D%20f%28t%29%3B%20%5Clim_%7Bx%20%5Crightarrow%20%5Cinfty%7D%20T%28x%2Ct%29%20%3D%200 ) 和初始条件 ( ![equation]https://latex.codecogs.com/gif.download?T%28x%2C0%29%20%3D%200 ) 相对于时间，问题的解由 ![equation]https://latex.codecogs.com/gif.download?%5Ctilde%7BT%7D%20%3D%20%5Ctilde%7Bf%7D%28%5Calpha%29%20%5Cexp%20%5B-%28i%20%5Calpha%20/%20%5Ckappa%29%5E%7B1/2%7D%20x%5D 给出，其中波浪号表示函数的傅立叶变换，alpha 是变换变量（频率），kappa 是热扩散系数。所以，方案是取边界数据的FFT，乘以 ![equation]https://latex.codecogs.com/gif.download?%5Cexp%20%5B%20-%20%28i%20%5Calpha%20/%20%5Ckappa%29%5E%7B1/2%7D%20x%5D ，并反转（使用 IFFT）以获得 T。

我构建了一个基准问题来验证我的编码。特别是，我考虑了强制函数 ![equation]https://latex.codecogs.com/gif.download?f%28t%29%20%3D%20a%20%5Cexp%20%5B%20-%20b%5E2%20%28t-t_0%29%5E2%5D ，其中 a、b 和 t_0 是固定常量。该函数可以代入经典的 Duhamel 解，结果可以通过数值积分来评估。此外，该函数的FT具有简单的封闭形式，可以代入上一段给出的解，并且可以通过直接数值积分对解的FT进行反演。这两个独立评估的结果是相同的，验证了上面给出的 FT 解决方案。

我已尝试为上述强制函数编写此问题，该函数由以固定时间间隔采样的离散数据给出。我调用Python的FFT例程来生成来自数据的强制函数的 FT。然后，我乘以 x 中的衰减指数，最后调用 IFFT 对给定的 x 固定值处的温度历史进行反演。该代码返回的值表现出一些正确的定性行为（例如，峰值温度的延迟到达；温度在后期逐渐下降），但在数量上（例如，温度过高）和质量上（例如，温度在早期不是渐近零）。我显然在对这个问题的编码中做错了什么，但我一直无法确定问题所在。我在衰减指数中指定频率的方式可能会发现一个可能的错误。

如果有人更熟悉 Python 的 FFT 和 IFFT 的使用，我将不胜感激。

我的代码如下：

def FTtrial3():
    import numpy as np
    import matplotlib.pyplot as plt
# nsam is the number of points given for the boundary data
    nsam = 32
    tmin = 0.
    tmax = 600.
    timedata = np.linspace(tmin,tmax,num=nsam)
# a,b,and t0 are fixed constants defining the boundary forcing
    a = 10.0
    b = 0.02
    time0 = 300.0 
    tempdata = a * np.exp(-b*b*(timedata - time0)**2) 
# take the FFT of the boundary data
    FTtemp = np.fft.fft(tempdata) 
# generate an array of the Fourier frequencies
    tint = (tmax-tmin)/(nsam-1)
    freq = np.fft.fftfreq(nsam,d=tint)
# specify the distance from the boundary for evaluation (in m)
    x = 0.1
# specify the thermal diffusivity (in m^2/s)
    D = 1.2E-05
# compute the FT of the solution
    FTsoln = FTtemp * np.exp(-np.sqrt( 1.j * freq / D) * x)
# invert to obtain the final solution
    soln = np.fft.ifft(FTsoln) 
    solnreal = np.real(soln)
# provide results of independent evaluation of Duhamel's solution for comparison
    pltdhml = np.array( [ [200.0,2.063E-05],[250.0,2.558E-03],[300.0,6.229E-02],[330.0,2.112E-01],[350.0,3.721E-01],[370.0,5.546E-01],[400.0,7.877E-01],[440.0,9.293E-01],[450.0,9.375E-01],[460.0,9.380E-01],[470.0,9.324E-01],[500.0,8.924E-01],[550.0,7.949E-01],[600.0,6.970E-01],[650.0,6.112E-01],[700.0,5.388E-01],[800.0,4.275E-01],[900.0,3.483E-01],[1000.0,2.902E-01] ] )
    tdhml,Tempdhml = pltdhml.T
# generate plot for comparison of FFT solution to Duhamel (exact) solution
    plt.scatter(tdhml,Tempdhml,color = 'green',marker = 'o',label = "Duhamel")
    plt.scatter(timedata,solnreal,color = 'red',marker = '^',label = "FT")
    plt.legend()
    plt.show()
FTtrial3()

解决方法

暂无找到可以解决该程序问题的有效方法，小编努力寻找整理中！

如果你已经找到好的解决方法，欢迎将解决方案带上本链接一起发送给小编。

小编邮箱:dio#foxmail.com (将#修改为@）

equation fft fft heat numpy python