我正在尝试使用复数到复数 IDFT 求解一维热方程。问题在于单个时间步后的输出似乎不正确。我在下面提供了一个简单的示例来说明该问题。
I initialize the temperature state as follows:
![Initial state of the temperature domain](https://i.stack.imgur.com/WIpqh.png)
频域中的初始模式为:
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
然后我将频域的状态推进到t=0.02
使用标准一维热方程:
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
}
频率模式为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:
![State of the spatial domain at t=0.02](https://i.stack.imgur.com/l0hZC.png)
空间域和频域似乎都是正确的周期性的。然而,热量(空间域中的值)似乎并未按照高斯曲线消散。更令人惊讶的是,一些温度低于其初始值(它们变为负值!)。
能量守恒似乎确实成立:将所有温度加在一起仍然是 100。
这是我的完整热方程代码:
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);
的定义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;
}
也许值得注意的是,我还使用复杂到真实的 IDFT 进行了这个实验(使用 fftw 的fftw_plan_dft_c2r_1d()
)并且给出了完全相同的结果。