並行譜數值方法/使用傅立葉譜方法求導數

使用傅立葉譜方法求導數

譜方法是一類數值技術，通常利用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.