c++ - 1D Heat Equation using DFT produces incorrect results (FFTW)
问题描述
I am trying to solve the 1D heat equation using a complex to complex IDFT. The problem is that the output after a single timestep does not seem to be correct. I have included a simple example below to illustrate the problem.
I initialize the temperature state as follows:
The initial modes in the frequency domain are:
k[ 0] = 12.5 + 0i
k[ 1] = 12.5 + 0i
k[ 2] = 12.5 + 0i
k[ 3] = 12.5 + 0i
k[ 4] = 12.5 + 0i
k[-3] = 12.5 + 0i
k[-2] = 12.5 + 0i
k[-1] = 12.5 + 0i
I then advance the state of the frequency domain to t=0.02
using the standard 1D heat equation:
double alpha = 0.2; // Thermal conductivity constant
double timestep = 0.02;
for (int i = 0; i < N; i++) {
int k = (i <= N / 2) ? i : i - N;
F[i][REAL] *= exp(-alpha * k * k * timestep); // Decay the real part
F[i][IMAG] *= exp(-alpha * k * k * timestep); // Decay the imaginary part
}
The frequency modes at t=0.02
become:
k[ 0] = 12.5 + 0i
k[ 1] = 12.45 + 0i
k[ 2] = 12.3 + 0i
k[ 3] = 12.06 + 0i
k[ 4] = 11.73 + 0i
k[-3] = 12.06 + 0i
k[-2] = 12.3 + 0i
k[-1] = 12.45 + 0i
After performing the IDFT to obtain the temperature domain state at t=0.02
I get:
The spatial and frequency domain both seem to be correctly periodic. However, the heat (values in spatial domain) does not seem to dissipate according to a Gaussian curve. Even more surprisingly, some temperatures dip below their initial value (they become negative!).
The conservation of energy does seem to hold correctly: adding all the temperatures together still nets 100.
This is my full heat equation code:
double alpha = 0.2; // Thermal conductivity constant
double timestep = 0.02; // Physical heat equation timestep
int N = 8; // Number of data points
fftw_complex* T = (fftw_complex*)fftw_alloc_complex(N); // Temperature domain
fftw_complex* F = (fftw_complex*)fftw_alloc_complex(N); // Frequency domain
fftw_plan plan = fftw_plan_dft_1d(N, F, T, FFTW_BACKWARD, FFTW_MEASURE); // IDFT from frequency to temperature domain
// Initialize all frequency modes such that there is a peak of 100 at x=0 in the temperature domain
// All other other points in the temperature domain are 0
for (int i = 0; i < N; i++) {
F[i][REAL] = 100.0 / N;
F[i][IMAG] = 0.0;
}
// Perform the IDFT to obtain the initial state in the temperature domain
fftw_execute(plan);
printTime1d(T, N);
printFrequencies1d(F, N);
// Perform a single timestep of the heat equation to obtain the frequency domain state at t=0.02
for (int i = 0; i < N; i++) {
int k = (i <= N / 2) ? i : i - N;
F[i][REAL] *= exp(-alpha * k * k * timestep); // Decay the real part
F[i][IMAG] *= exp(-alpha * k * k * timestep); // Decay the imaginary part
}
// Perform the IDFT to obtain the temperature domain state at t=0.02
fftw_execute(plan);
printTime1d(T, N);
printFrequencies1d(F, N);
The definition of printTime(...)
and printFrequencies(...)
is:
void printTime1d(fftw_complex* data, int N) {
int rounding_factor = pow(10, 2);
for (int i = 0; i < N; i++) {
std::cout << std::setw(8) << round(data[i][REAL] * rounding_factor) / rounding_factor;
}
std::cout << std::endl;
}
void printFrequencies1d(fftw_complex* data, int N) {
int rounding_factor = pow(10, 2);
for (int i = 0; i < N; i++) {
int k = (i <= N / 2) ? i : i - N;
double R = round(data[i][REAL] * rounding_factor) / rounding_factor;
double I = round(data[i][IMAG] * rounding_factor) / rounding_factor;
std::cout << "k[" << std::setw(2) << k << "]: " << std::setw(2) << R << ((I < 0) ? " - " : " + ") << std::setw(1) << abs(I) << "i" << std::endl;
}
std::cout << std::endl;
}
Perhaps good to note is that I have also conducted this experiment using a complex to real IDFT (with fftw's fftw_plan_dft_c2r_1d()
) and it gave the exact same results.
解决方案
Your problem is that you don't resolve the needed frequencies, getting instead the following Fourier image of the function after multiplication by the decay coefficients:
The above result is too far from what you should get instead—a Gaussian—at least something like this (using 80 points instead of 8):
Notice how the amplitudes in the first chart above don't even have any chance to come even close to zero, instead bumping into the Nyquist frequency. It's then obvious that you'll get artifacts resembling the Gibbs phenomenon: it's the usual behavior of Fourier partial sums.
The inverse Fourier transform of the 80-point data version is then as follows:
This result still does have negative components (since we use a finite number of harmonics), but these are much smaller in amplitude than what you got with only 8 harmonics.
Note that this does mean that, if you increase the value of time you're interested in, you could reduce the number of harmonics taken into account. This might be unexpected at first, but it's simply because the upper harmonics decay much faster than the lower ones, and they never ever increase back.
推荐阅读
- ssl - 端口 443 中的 tomcat ssl 配置
- javascript - 结合 CSS 和 SVG 过滤器来编辑和保存图像
- node.js - 如何在亚马逊 ses 上使用 sendBulkEmail 发送个性化电子邮件?
- matlab - MATLAB。位置 2 的索引超出数组边界(不得超过 1)
- javascript - [react]:useMemo 在传递数组时返回第一个元素
- json - Gson 反序列化工作异常
- spring-boot - @PreUpdate 在获取操作期间被调用
- html - 引导按钮未调整大小
- r - 模块化具有多个输出 ID 的闪亮代码
- java - Java 1.7.0_67 api 连接