并行谱数值方法/使用傅里叶谱方法求导数

使用傅里叶谱方法求导数

谱方法是一类数值技术，通常利用FFT。谱方法可以在Matlab中轻松实现，但需要注意一些约定。首先要注意，Matlab的“fft”和“ifft”函数存储波数的顺序与迄今为止使用的顺序不同。Matlab和其他大多数FFT包中的波数排序为$0,1,...,{\frac {n}{2}},-{\frac {n}{2}}+1,-{\frac {n}{2}}+2,...,-1$。其次，Matlab没有充分利用真实输入数据。真实数据$u$的DFT满足对称性${\hat {u}}(-k)={\hat {u}}(k)$，因此只需要计算一半的波数。Matlab的“fft”命令没有充分利用此特性，并浪费内存存储正负波数。第三，对于平滑函数，谱精度（傅里叶系数幅值的指数衰减）更好，因此在可能的情况下，确保您的初始条件是平滑的 - **使用傅里叶谱方法时，这要求您的初始条件是周期性的**。

在傅里叶空间求导数

令$u_{j}$为函数$u(x)$的离散近似，该函数在$n$个离散点$x_{j}\in {h,2h,...,ih,..,2\pi -h,2\pi }$处采样，其中$h=2\pi /n$。现在对$u_{j}$进行傅里叶变换，使得${\text{FFT}}(u_{j})\equiv {\hat {u}}_{k}\qquad {\text{where}}\quad k\in {{\frac {-n}{2}}+1,...{\frac {n}{2}}}.$。然后可以从${\hat {u}}_{k}$计算${\frac {\partial ^{\nu }u_{j}}{\partial x^{\nu }}}$的傅里叶变换，如下所示

${\text{FFT}}({\frac {\partial ^{\nu }u_{j}}{\partial x^{\nu }}})\equiv (ik)^{\nu }{\hat {u}}_{k}$

( 1)

其中，当 $\nu$ 为奇数时，${\hat {u}}_{n/2}=0$ 。更多细节可以在 Trefethen 中找到。

因此，实空间中的微分变成了傅里叶空间中的乘法。然后，我们可以采用逆快速傅里叶变换 (IFFT) 来得到实空间的解。在下一节中，我们将使用这种技术来实现向前欧拉和向后欧拉时间步进方案，以计算几个 PDE 的解。

( A)

一个使用傅里叶变换来计算 $f(x)=\sin ^{2}\left(wx\right)$ 在 $(1,1+2\pi ]$ 上的前两个导数的 Matlab 程序。

% Approximates the derivative of a periodic function f(x) using Fourier
% transforms.  Requires a linear discretization and a domain (a,b] such
% that f(x) = f(x+b-a)

% domain
a = 1;
b = 1 + pi/2;
N = 100;
dx = (b-a)/N;
x = a + dx*(0:N-1);

% function
w = 2;
f = sin(w*x).^2;

% exact derivatives
dfdx = 2*w*sin(w*x).*cos(w*x);
d2fdx2 = 4*w^2*cos(w*x).^2 - 2*w^2;

% fourier derivatives
Nx = size(x,2);
k = 2*pi/(b-a)*[0:Nx/2-1 0 -Nx/2+1:-1];
dFdx = ifft(1i*k.*fft(f));
d2Fdx2 = ifft(-k.^2.*fft(f));

% graph result
figure;
plot(x,dfdx,'r-',x,d2fdx2,'g-',x,dFdx,'k:',x,d2Fdx2,'b:','LineWidth',2);
legend('df/dx','d^2f/dx^2','Fourier df/dx','Fourier d^2f/dx^2');

练习

1. 设 $u(x)=\sum _{k}{\hat {u}}_{k}\exp(ikx)$ 是函数 $u(x)$ 的傅里叶级数表示。解释为什么 ${\frac {\mathrm {d} ^{\nu }u}{\mathrm {d} x^{\nu }}}=\sum (ik)^{\nu }{\hat {u}}_{k}\exp(ikx),$ 前提是级数收敛。

2. 考虑线性 KdV 方程 $u_{t}+u_{xxx}=0$ ，它对 $x\in (0,2\pi ]$ 具有周期边界条件，初始数据为

   $u(x,0)=0\quad {\text{ if }}\quad 0<x\leq \pi$

   和

   $u(x,0)=1\quad {\text{ if }}\quad \pi <x\leq 2\pi$

3. 使用变量分离法，证明该方程的“解”为 $u(t,x)={\frac {1}{2}}-{\frac {2}{\pi }}\sum _{j=0}^{\infty }{\frac {\sin((2j+1)x-(2j+1)^{3}t)}{2j+1}}.$ 使用引号是因为这个解的表达式在时间上微分一次或空间上微分两次时不收敛。

4. 正如 Olver 所解释的，这个解在时间为 π 的无理数倍时具有分形结构，在时间为 π 的有理数倍时具有量子化结构。$\pi$ Listing B 中的 Matlab 程序使用快速傅里叶变换来找到线性化 KdV 方程的解。解释该程序是如何找到线性化 KdV 方程的解。

5. 将 Matlab 程序产生的数值解与解析解进行比较。尝试确定哪一个更准确，并尝试找到证据或解释来支持你的建议。

( B)

Matlab 程序 代码下载

% This program computes the solution to the linearly dispersive
% wave equation using the Fast Fourier Transform

N = 512; % Number of grid points.
h = 2*pi/N; % Size of each grid.
x = h*(1:N); % Variable x as an array.
t = .05*pi; % Time to plot solution at
dt = .001; % Appropriate time step.
u0 = zeros(1,N); % Array to hold initial data
u0(N/2+1:N)= ones(1,N/2); % Defining the initial data
k=(1i*[0:N/2-1 0 -N/2+1:-1]); % Fourier wavenumbers
k3=k.^3;
u=ifft(exp(k3*t).*fft(u0)); % Calculate the solution
plot(x,u,'r-'); % Plot the solution
xlabel x; ylabel u; % Label the axes of the graphs
title(['Time ',num2str(t/(2*pi)),' \pi']);

参考文献

Olver, P.J. (2010). "Dispersive Quantization". American Mathematical Monthly. 117: 599–610. {{cite journal}}: Cite has empty unknown parameter: |coauthors= (help)

Trefethen, L.N. (2000). Spectral Methods in Matlab. SIAM.