Matlab 和 Python 中的示例

現在我們想使用傅立葉譜方法找到近似的數值解。在本節中，我們主要關注具有周期邊界條件的熱方程，其中 $x\in [0,2\pi )$ 。這裡使用的許多技術也適用於更復雜的偏微分方程，這些方程無法直接使用變數分離來求解。

一維熱方程

{\frac {\partial u}{\partial t}}=\alpha {\frac {\partial ^{2}u}{\partial x^{2}}}

( 1)

是一個眾所周知的二階偏微分方程，可以使用變數分離法找到其精確的級數解。它出現在幾個情況下，例如預測細長均勻橫截面杆中的溫度。該方程及其推導可以在偏微分方程和微積分的入門書籍中找到，例如和。常數 $\alpha$ 是熱擴散率， $u(x,t)$ 是溫度。我們已經描述瞭如何使用變數分離法求解熱方程。首先讓我們離散化 $x$ ，使得 $x_{j}$ ，其中 $j=0,1,2,...,n$ 。 $x_{j}$ 在 $[0,2\pi )$ 中均勻分佈。現在讓我們對一維熱方程的兩邊進行 FFT，得到

${\widehat {\frac {\partial u}{\partial t}}}=\alpha {\widehat {\frac {\partial ^{2}u}{\partial x^{2}}}}.$

然後我們使用方程 ${\text{FFT}}({\frac {\partial ^{\nu }u_{j}}{\partial x^{\nu }}})\equiv (ik)^{\nu }{\hat {u}}_{k}.$ 重新寫出空間導數。

下標 $k$ 表示 $k^{th}$ 傅立葉模的係數。

${\frac {\partial {\hat {u}}_{k}}{\partial t}}=\alpha (ik)^{2}{\hat {u}}_{k},$

因此，偏微分方程現在變成了獨立的常微分方程的集合。雖然我們可以精確地求解這些常微分方程的時間解，但我們將使用一些技術，這些技術也將使我們能夠獲得對我們無法精確求解的偏微分方程的近似解。我們將討論兩種求解這些常微分方程的方法：前向尤拉和後向尤拉。

前向尤拉

使用時間上的前向尤拉方法，我們得到

${\frac {{\hat {u}}_{k}^{n+1}-{\hat {u}}_{k}^{n}}{h}}=\alpha (ik)^{2}{\hat {u}}_{k}^{n}$

${\hat {u}}_{k}^{n+1}={\hat {u}}_{k}^{n}+\alpha h(ik)^{2}{\hat {u}}_{k}^{n}$

剩下的就是對所有時間步長完成後計算出的解進行逆FFT變換，將其轉換回實空間。這是一個線性偏微分方程，所以最後只需要一個逆FFT變換。我們稍後將看到，對於非線性偏微分方程，情況並非如此。該方法的 Matlab 實現見列表 A 。

( A)

使用前向尤拉時間步長求解熱方程的 Matlab 程式程式碼下載

%Solving Heat Equation using pseudo-spectral and Forward Euler
%u_t= \alpha*u_xx
%BC= u(0)=0, u(2*pi)=0
%IC=sin(x)
clear all; clc;

%Grid
N = 64; %Number of steps
h = 2*pi/N; %step size
x = h*(1:N); %discretize x-direction

alpha = .5; %Thermal Diffusivity constant
t = 0;
dt = .001;

%Initial conditions
v = sin(x);
k=(1i*[0:N/2-1 0 -N/2+1:-1]);
k2=k.^2;

%Setting up Plot
tmax = 5; tplot = .1;
plotgap= round(tplot/dt);
nplots = round(tmax/tplot);
data = [v; zeros(nplots,N)]; tdata = t;

  
for i = 1:nplots
    v_hat = fft(v); %Fourier Space
    for n = 1:plotgap
        v_hat = v_hat+dt*alpha*k2.*v_hat; %FE timestepping
    end
    v = real(ifft(v_hat)); %Back to real space
    data(i+1,:) = v;
    t=t+plotgap*dt;
    tdata = [tdata; t]; %Time vector
end

%Plot using mesh
mesh(x,tdata,data), grid on,
view(-60,55), xlabel x, ylabel t, zlabel u, zlabel u

( Ap)

使用前向尤拉時間步長求解熱方程的 Python 程式。程式碼下載

#!/usr/bin/env python
"""
Solving Heat Equation using pseudo-spectral and Forward Euler
u_t= \alpha*u_xx
BC= u(0)=0, u(2*pi)=0
IC=sin(x)
"""

import math
import numpy
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
from matplotlib.ticker import LinearLocator

# Grid
N = 64 # Number of steps
h = 2*math.pi/N # step size
x = h*numpy.arange(0,N) # discretize x-direction

alpha = 0.5 # Thermal Diffusivity constant
t = 0
dt = .001

# Initial conditions
v = numpy.sin(x)
I = complex(0,1)
k = numpy.array([I*y for y in range(0,N/2) + [0] + range(-N/2+1,0)])
k2=k**2;

# Setting up Plot
tmax = 5; tplot = .1;
plotgap = int(round(tplot/dt))
nplots = int(round(tmax/tplot))

data = numpy.zeros((nplots+1,N))
data[0,:] = v
tdata = [t]

for i in xrange(nplots):
    v_hat = numpy.fft.fft(v)

    for n in xrange(plotgap):
        v_hat = v_hat+dt*alpha*k2*v_hat # FE timestepping

    v = numpy.real(numpy.fft.ifft(v_hat)) # back to real space
    data[i+1,:] = v

    # real time vector
    t = t+plotgap*dt
    tdata.append(t)

# Plot using mesh
xx,tt = (numpy.mat(A) for A in (numpy.meshgrid(x,tdata)))

fig = plt.figure()
ax = fig.gca(projection='3d')
surf = ax.plot_surface(xx, tt, data,rstride=1, cstride=1, cmap=cm.jet,
        linewidth=0, antialiased=False)
fig.colorbar(surf, shrink=0.5, aspect=5)
plt.xlabel('x')
plt.ylabel('t')
plt.show()

後向尤拉

為了推匯出這種方法，我們首先應用 FFT 變換，然後使用後向尤拉進行時間步長。然後，我們將隱式形式改寫成給出下一迭代值的形式，

${\frac {\partial {\hat {u}}_{k}}{\partial t}}=\alpha (ik)^{2}{\hat {u}}_{k}$

${\frac {{\hat {u}}_{k}^{n+1}-{\hat {u}}_{k}^{n}}{h}}=\alpha (ik)^{2}{\hat {u}}_{k}^{n+1}$

${\hat {u}}_{k}^{n+1}(1-\alpha h(ik)^{2})={\hat {u}}_{k}^{n}$

${\hat {u}}_{k}^{n+1}={\frac {{\hat {u}}_{k}^{n}}{(1-\alpha h(ik)^{2})}}.$ 圖 i 展示了熱方程的一個數值解（獲得精確解的方法可以在 Boyce 和 DiPrima 等文獻中找到^[1]）。其中 $n=64$ 是使用列表 B 中的 Matlab 程式獲得的。

( B)

使用後向尤拉時間步進法求解熱方程的 Matlab 程式程式碼下載

%Solving Heat Equation using pseudospectral methods with Backwards Euler:
%u_t= \alpha*u_xx
%BC = u(0)=0 and u(2*pi)=0 (Periodic)
%IC=sin(x)
clear all; clc;

%Grid
N = 64; h = 2*pi/N; x = h*(1:N);

% Initial conditions
v = sin(x);
alpha = .5;
t = 0;
dt = .001; %Timestep size

%(ik)^2 Vector
k=(1i*[0:N/2-1 0 -N/2+1:-1]);
k2=k.^2;

%Setting up Plot
tmax = 5; tplot = .1;
plotgap= round(tplot/dt);
nplots = round(tmax/tplot);
data = [v; zeros(nplots,N)]; tdata = t;

  
for i = 1:nplots
    v_hat = fft(v); %Converts to fourier space
    for n = 1:plotgap
        v_hat = v_hat./(1-dt*alpha*k2); %Backwards Euler timestepping
    end
    v = ifft(v_hat); %Converts back to real Space
    data(i+1,:) = real(v); %Records data
    t=t+plotgap*dt; %Records time
    tdata = [tdata; t];
end

%Plot using mesh
mesh(x,tdata,data), grid on, %axis([-1 2*pi 0 tmax -1 1]),
view(-60,55), xlabel x, ylabel t, zlabel u, zlabel u,

( Bp)

使用後向尤拉時間步進法求解熱方程的 Python 程式。程式碼下載

#!/usr/bin/env python
"""
Solving Heat Equation using pseudospectral methods with Backwards Euler:
u_t= \alpha*u_xx
BC = u(0)=0 and u(2*pi)=0 (Periodic)
IC=sin(x)
"""

import math
import numpy
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
from matplotlib.ticker import LinearLocator

# Grid
N = 64; h = 2*math.pi/N; x = [h*i for i in xrange(1,N+1)]

# Initial conditions
v = [math.sin(y) for y in x]
alpha = 0.5
t = 0
dt = .001 #Timestep size

# (ik)^2 Vector
I = complex(0,1)
k = numpy.array([I*n for n in range(0,N/2) + [0] + range(-N/2+1,0)])
k2=k**2;

# Setting up Plot
tmax = 5.0; tplot = 0.1
plotgap= int(round(tplot/dt))
nplots = int(round(tmax/tplot))
data = numpy.zeros((nplots+1,N))
data[0,:] = v
tdata = [t]

for i in xrange(nplots):
    v_hat = numpy.fft.fft(v) # convert to fourier space
    for n in xrange(plotgap):
        v_hat = v_hat / (1-dt*alpha*k2) # backward Euler timestepping

    v = numpy.fft.ifft(v_hat) # convert back to real space
    data[i+1,:] = numpy.real(v) # records data

    t = t+plotgap*dt # records real time
    tdata.append(t)

# Plot using mesh
xx,tt = (numpy.mat(A) for A in (numpy.meshgrid(x,tdata)))

fig = plt.figure()
ax = fig.gca(projection='3d')
surf = ax.plot_surface(xx, tt, data,rstride=1, cstride=1, cmap=cm.jet,
        linewidth=0, antialiased=False)
fig.colorbar(surf, shrink=0.5, aspect=5)
plt.xlabel('x')
plt.ylabel('t')
plt.show()

練習

編寫程式使用 Crank-Nicolson 方法求解熱方程。
求解對流方程 $u_{t}=u_{x}$ 對於 $x\in [0,2\pi )$ ，初始條件為
1. $u(t=0,x)=\cos(x)$
2. $u(t=0,x)=0$ $u(t=0,x)=0$ 對於 $x<\pi$ 以及 $u(t=0,x)=1$ 對於 $x\geq \pi$
直到時間 $T=1$ 。你可以透過使用變數分離或假設解的形式為 $u(x,t)=f(x+t)$ 來實現這一點，並推匯出 $f$ 是什麼來滿足初始條件。在這兩種情況下，請使用前向尤拉、後向尤拉和 Crank-Nicolson 時間步進方案。在計算每種情況下的精確解之後，檢查這三種方法中，最終時間最大誤差如何依賴於時間步長。

非線性方程

一維 Allen-Cahn 方程

到目前為止，我們只處理線性方程。現在是時候討論非線性偏微分方程了。Allen-Cahn 方程模擬了材料中相分離的過程。它是由 S. Allen 和 J. W. Cahn 提出的^[2]，其表示式為

{\frac {\partial u}{\partial t}}=\epsilon {\frac {\partial ^{2}u}{\partial x^{2}}}+u-u^{3},

( 2)

其中 $\epsilon$ 是一個很小的正常數。數值求解該方程的方法與熱方程的求解方法類似，但也存在一些顯著區別。最大的區別在於 FFT( $u^{3}$ ) $\neq ({\text{FFT}}(u))^{3}$ ，因此必須先計算 $u^{3}$ ，然後再進行 FFT 變換。FFT 是一個線性運算，而立方運算是非線性運算，因此運算順序很重要

{\frac {\partial {\hat {u}}_{k}}{\partial t}}=\epsilon {\frac {\partial ^{2}{\hat {u}}_{k}}{\partial x^{2}}}+{\hat {u}}_{k}-{\widehat {u^{3}}}_{k}.

( 3)

接下來，像我們之前在熱方程中所做的那樣，對等式右邊第一項進行重寫

${\frac {\partial {\hat {u}}_{k}}{\partial t}}=\epsilon (ik)^{2}{\hat {u}}_{k}+{\hat {u}}_{k}-{\widehat {u^{3}}}_{k}.$

為了數值求解該方程，我們將使用隱式（後向尤拉）和顯式（前向尤拉）方法的組合。我們將跳過前向尤拉方法，因為它不穩定。

隱式-顯式方法

您可能已經注意到，由於非線性項的存在，後向尤拉方法目前無法用於 Allen-Cahn 方程。如果您嘗試實現後向尤拉方法，您會發現無法提取所有 ${\hat {u}}_{k}^{n+1}$ 。幸運的是，有一個簡單直觀的解決方法，不會對解的精度造成負面影響。除了非線性項外，將所有項都隱式地寫出（後向尤拉）。將此應用於 Allen-Cahn 方程，我們發現 ^[3]

${\frac {{\hat {u}}_{k}^{n+1}-{\hat {u}}_{k}^{n}}{h}}=\epsilon (ik)^{2}{\hat {u}}_{k}^{n+1}+{\hat {u}}_{k}^{n}-{\widehat {(u^{n})^{3}}}_{k}$

${\hat {u}}_{k}^{n+1}\left(-\epsilon (ik)^{2}+{\frac {1}{h}}\right)={\frac {1}{h}}{\hat {u}}_{k}^{n}+{\hat {u}}_{k}^{n}-{\widehat {(u^{n})^{3}}}_{k}$

${\hat {u}}_{k}^{n+1}={\frac {{\hat {u}}_{k}^{n}({\frac {1}{h}}+1)-{\widehat {(u^{n})^{3}}}_{k}}{\left(-\epsilon (ik)^{2}+{\frac {1}{h}}\right)}}.$

現在我們有了可以處理的形式。我們可以將初始條件設定為 $u(x,0)={\frac {1}{4}}\sin(x)$ ，並繪製由列表 C 中的 Matlab 程式碼計算的空間時間演化。計算結果如圖 ii 所示。

( C)

使用隱式顯式時間步進求解一維 Allen-Cahn 方程的 Matlab 程式程式碼下載

%Solving 1D Allen-Cahn Eq using pseudo-spectral and Implicit/Explicit method
%u_t=u_{xx} + u - u^3
%where u-u^3 is treated explicitly and u_{xx} is treated implicitly
%BC = u(0)=0, u(2*pi)=0 (Periodic)
%IC=.25*sin(x);
clear all; clc;

%Grid and Initial Data
N = 8000; h = 2*pi/N; x = h*(1:N); t = 0;

dt = .001; %timestep size
epsilon= .001;

%initial conditions
v = .25*sin(x);

%(ik) and (ik)^2 vectors
k=(1i*[0:N/2-1 0 -N/2+1:-1]);
k2=k.^2;

%setting up plot
tmax = 5; tplot = .2;
plotgap= round(tplot/dt);
nplots = round(tmax/tplot);
data = [v; zeros(nplots,N)]; tdata = t;
  
for i = 1:nplots
    for n = 1:plotgap
        v_hat = fft(v); %converts to Fourier space
        vv = v.^3; %computes nonlinear term in real space
        vv = fft(vv); %converts nonlinear term to Fourier space
        v_hat = (v_hat*(1/dt+1) - vv)./(1/dt-k2*epsilon); %Implicit/Explicit
        v = ifft(v_hat); %Solution back to real space
    end
    data(i+1,:) = real(v); %Records data each "plotgap"
    t=t+plotgap*dt; %Real time
    tdata = [tdata; t];
end

mesh(x,tdata,data), grid on, axis([-1 2*pi 0 tmax -1 1]),
view(-60,55), xlabel x, ylabel t, zlabel u

二維 Allen-Cahn 方程

現在我們將看一下 Allen-Cahn 方程的二維形式，其中 $u(x,y,t)$ 滿足

{\frac {\partial u}{\partial t}}=\epsilon \left({\frac {\partial ^{2}u}{\partial x^{2}}}+{\frac {\partial ^{2}u}{\partial y^{2}}}\right)+u-u^{3}.

( 4)

透過對兩邊進行 FFT 將其轉換為傅立葉空間

${\frac {\partial {\hat {u}}_{k}}{\partial t}}=\epsilon \left({\frac {\partial ^{2}{\hat {u}}_{k}}{\partial x^{2}}}+{\frac {\partial ^{2}{\hat {u}}_{k}}{\partial y^{2}}}\right)+{\hat {u}}_{k}-{\widehat {u^{3}}}_{k}$

{\frac {\partial {\hat {u}}_{k}}{\partial t}}=\epsilon \left((ik_{x})^{2}{\hat {u}}_{k}+(ik_{y})^{2}{\hat {u}}_{k}\right)+{\hat {u}}_{k}-{\widehat {(u^{3})}}_{k}

( 5)

其中 $k_{x}$ 和 $k_{y}$ 用於提醒我們在相應方向進行 FFT 變換。我們還定義

f(u)\equiv u-u^{3}

( 6)

處理右側前兩項的方法是先對 x 方向進行 FFT 變換，再對 y 方向進行 FFT 變換。FFT 變換的順序， $x$ 優先或 $y$ 優先並不重要。一些軟體庫提供二維 FFT 變換。使用多維 FFT 還是多個一維 FFT 更有效取決於所求解的方程。通常，使用多維 FFT 編寫程式更容易，但在某些情況下效率不高，因為可以立即重用剛剛進行傅立葉變換的資料。

隱式-顯式方法

在這種方法中，公式 6 中給出的非線性項被顯式計算，而其餘項將被隱式寫成如下形式：

${\frac {{\hat {u}}_{k}^{n+1}-{\hat {u}}_{k}^{n}}{h}}=\epsilon \left((ik_{x})^{2}{\hat {u}}_{k}^{n+1}+(ik_{y})^{2}{\hat {u}}_{k}^{n+1}\right)+{\widehat {f(u^{n})}}_{k}$

${\hat {u}}_{k}^{n+1}\left(-\epsilon (ik_{x})^{2}-\epsilon (ik_{y})^{2}+{\frac {1}{h}}\right)={\frac {{\hat {u}}_{k}^{n}}{h}}+{\widehat {f(u^{n})}}_{k}$

${\hat {u}}_{k}^{n+1}={\frac {{\frac {{\hat {u}}_{k}^{n}}{h}}+{\widehat {f(u^{n})}}_{k}}{\left(-\epsilon (ik_{x})^{2}-\epsilon (ik_{y})^{2}+{\frac {1}{h}}\right)}}$

然後我們可以用 $f(u)$ 代替。

{\hat {u}}_{k}^{n+1}={\frac {{\hat {u}}_{k}^{n}\left({\frac {1}{h}}+1\right)-{\widehat {(u^{n})^{3}}}_{k}}{\left(-\epsilon (ik_{x})^{2}-\epsilon (ik_{y})^{2}+{\frac {1}{h}}\right)}}.

( 7)

用於生成圖 iii 的 Matlab 程式碼在列表 D 中。

( D)

使用隱式顯式時間步進法求解 2D Allen-Cahn 方程的 Matlab 程式程式碼下載

%Solving 2D Allen-Cahn Eq using pseudo-spectral with Implicit/Explicit
%u_t= epsilon(u_{xx}+u_{yy}) + u - u^3
%where u-u^3 is treated explicitly and epsilon(u_{xx} + u_{yy}) is treated implicitly
%BC = Periodic
%IC=v=sin(2*pi*x)+0.001*cos(16*pi*x;
clear all; clc;

%Grid
N = 256; h = 1/N; x = h*(1:N);
dt = .01;

%x and y meshgrid
y=x';
[xx,yy]=meshgrid(x,y);

%initial conditions
v=sin(2*pi*xx)+0.001*cos(16*pi*xx);
epsilon=.01;

%(ik) and (ik)^2 vectors in x and y direction
kx=(2*pi*1i*[0:N/2-1 0 -N/2+1:-1]);
ky=(2*pi*1i*[0:N/2-1 0 -N/2+1:-1]');
k2x=kx.^2;
k2y=ky.^2;

[kxx,kyy]=meshgrid(k2x,k2y);
        
for n = 1:500
    v_nl=v.^3; %calculates nonlinear term in real space
    %FFT for linear and nonlinear term
    v_nl = fft2(v_nl);
    v_hat=fft2(v);
    vnew=(v_hat*(1+1/dt)-v_nl)./ ...
       (-(kxx+kyy)*epsilon+1/dt); %Implicit/Explicit timestepping
    %converts to real space in x-direction
    v=ifft2(vnew);
    %Plots each timestep
    mesh(v); title(['Time ',num2str(n)]); axis([0 N 0 N -1 1]);
    xlabel x; ylabel y; zlabel u;
    view(43,22); drawnow;
end

( Dp)

使用隱式顯式時間步進法求解 2D Allen-Cahn 方程的 Python 程式。程式碼下載

#!/usr/bin/env python
"""
Solving 2D Allen-Cahn Eq using pseudo-spectral with Implicit/Explicit
u_t= epsilon(u_{xx}+u_{yy}) + u - u^3
where u-u^3 is treated explicitly and u_{xx} and u_{yy} is treated implicitly
BC = Periodic
IC=v=sin(2*pi*x)+0.5*cos(4*pi*y)
"""

import math
import numpy
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
from matplotlib.ticker import LinearLocator
import time
 
plt.ion()

# Setup the grid
N = 64; h = 1.0/N;
x = [h*i for i in xrange(1,N+1)]
y = [h*i for i in xrange(1,N+1)]
dt = 0.05
xx,yy = (numpy.mat(A) for A in (numpy.meshgrid(x,y)))

# Initial Conditions
u = numpy.array(numpy.sin(2*math.pi*xx) + 0.5*numpy.cos(4*math.pi*yy), dtype=float)

epsilon = 0.01

# (ik) and (ik)^2 vectors in x and y direction
I = complex(0,1)
k_x = 2*numpy.pi*numpy.array([I*n for n in range(0,N/2) + [0] + range(-N/2+1,0)])
k_y = k_x

kxx = numpy.zeros((N,N), dtype=complex)
kyy = numpy.zeros((N,N), dtype=complex)
for i in xrange(N):
    for j in xrange(N):
        kxx[i,j] = k_x[i]**2
        kyy[i,j] = k_y[j]**2
    
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
surf = ax.plot_surface(xx, yy, u,rstride=1, cstride=1, cmap=cm.jet,
        linewidth=0, antialiased=False)
fig.colorbar(surf, shrink=0.5, aspect=5)
plt.xlabel('x')
plt.ylabel('y')
plt.show()

v_hat = numpy.zeros((N,N), dtype=complex)
v_hat = numpy.fft.fft2(u)

for n in xrange(100):
    # calculate nonlinear term in real space
    v_nl = numpy.array(u**3, dtype=complex)
# FFT for nonlinear and linear terms
    v_nl = numpy.fft.fft2(v_nl)
    v_hat = (v_hat*(1+1/dt) - v_nl)
    v_hat=v_hat/(1/dt - (kxx+kyy)*epsilon) # Implicit/Explicit timestepping
    u = numpy.real(numpy.fft.ifft2(v_hat))
    # Remove old plot before drawing
    ax.collections.remove(surf)
    surf = ax.plot_surface(xx, yy, u,rstride=1, cstride=1, cmap=cm.jet,
         linewidth=0, antialiased=False)
    plt.draw()
plt.show()

練習

許多練習來自 Uecker^[4]。有關這些方程的另一個入門資料來源是 Trefethen 和 Embree^[5]。

Burgers 方程表示為： ${\frac {\partial u}{\partial t}}=\nu {\frac {\partial ^{2}u}{\partial x^{2}}}-u{\frac {\partial u}{\partial x}}$ 其中 $\nu \in \mathbb {R} ^{+}$ 且 $u$ 具有周期邊界條件。使用隱式顯式方法求解該方程。如果將 $\nu$ 取為較小的值，請確保使用足夠的網格點以獲得正確的數值解。檢查這一點的一種簡單方法是不斷增加網格點的數量並檢查解是否沒有變化。另一種檢查方法是計算傅立葉係數並檢查最高的係數是否衰減到機器精度。
Kuramoto-Sivashinsky 方程表示為： ${\frac {\partial u}{\partial t}}=-{\frac {\partial ^{2}u}{\partial x^{2}}}-{\frac {\partial ^{4}u}{\partial x^{4}}}-u{\frac {\partial u}{\partial x}}$ ${\frac {\partial u}{\partial t}}=-{\frac {\partial ^{2}u}{\partial x^{2}}}-{\frac {\partial ^{4}u}{\partial x^{4}}}-u{\frac {\partial u}{\partial x}}$ 其中 $u$ $u$ 具有周期邊界條件。
- 該方程式模擬了什麼，你預計它的解會表現出什麼樣的行為？
- 使用隱式-顯式方法求解該方程式的數值解。
一維 Gray-Scott 方程組如下： ${\frac {\partial u}{\partial t}}=d_{1}{\frac {\partial ^{2}u}{\partial x^{2}}}-uv^{2}+f(1-u),$ ${\frac {\partial u}{\partial t}}=d_{1}{\frac {\partial ^{2}u}{\partial x^{2}}}-uv^{2}+f(1-u),$ ${\frac {\partial v}{\partial t}}=d_{2}{\frac {\partial ^{2}v}{\partial x^{2}}}+uv^{2}-(f+k)v$ ${\frac {\partial v}{\partial t}}=d_{2}{\frac {\partial ^{2}v}{\partial x^{2}}}+uv^{2}-(f+k)v$ 其中 $d_{1}$ $d_{1}$ ， $d_{2}$ $d_{2}$ ， $f$ $f$ 和 $k$ $k$ 為常數。
- 該方程式模擬了什麼，你預計它的解會表現出什麼樣的行為？
- 使用隱式-顯式方法求解該方程式的數值解。嘗試幾個不同的 $d_{1}$ ， $d_{2}$ ， $f$ 和 $k$ 的值，並將所得模式與文獻中的發現進行比較。
二維 Swift-Hohenberg 方程如下： ${\frac {\partial u}{\partial t}}=-\Delta ^{2}u+2\Delta u+(\alpha -1)u-u^{3},$ ${\frac {\partial u}{\partial t}}=-\Delta ^{2}u+2\Delta u+(\alpha -1)u-u^{3},$
- 該方程式模擬了什麼，你預計它的解會表現出什麼樣的行為？
- 使用隱式-顯式方法求解該方程式的數值解，並嘗試幾個不同的 $\alpha$ 值。
二維 Gray-Scott 方程組如下： ${\frac {\partial u}{\partial t}}=d_{1}\Delta u-uv^{2}+f(1-u)$ ${\frac {\partial u}{\partial t}}=d_{1}\Delta u-uv^{2}+f(1-u)$ ${\frac {\partial v}{\partial t}}=d_{2}\Delta v+uv^{2}-(f+k)v$ ${\frac {\partial v}{\partial t}}=d_{2}\Delta v+uv^{2}-(f+k)v$ 其中 $d_{1}$ $d_{1}$ ， $d_{2}$ $d_{2}$ ， $f$ $f$ 和 $k$ $k$ 為常數。
- 該方程式模擬了什麼，你預計它的解會表現出什麼樣的行為？
- 使用隱式-顯式方法求解該方程式的數值解。
二維復 Ginzburg-Landau 方程表示如下： ${\frac {\partial A}{\partial t}}=A+(1+i\alpha )\Delta A-(1+i\beta )|A|^{2}A.$ ${\frac {\partial A}{\partial t}}=A+(1+i\alpha )\Delta A-(1+i\beta )|A|^{2}A.$ 您可以透過以下連結檢視該方程的入門教程：http://codeinthehole.com/static/tutorial/index.html
- 該方程式模擬了什麼，你預計它的解會表現出什麼樣的行為？
- 使用隱式-顯式方法，針對 $\alpha$ 和 $\beta$ 的多個值，找到該方程的數值解。

註釋

↑ Boyce 和 DiPrima (2010)
↑ Allen 和 Cahn (1979)
↑ 請注意，在程式設計時，您需要在每次計算下一個時間步長 $n+1$ 時更新非線性項 ( $u^{3}$ )。值得一提的是，對於每個時間步長，您需要從實空間轉換到傅立葉空間，再轉換回實空間，然後重複此過程。而對於熱方程，您可以在傅立葉空間中執行任意數量的時間步長，只有在記錄資料時才轉換回來。
↑ Uecker (2009)
↑ Trefethen 和 Embree

參考文獻

Allen, S.M.; Cahn, J.W. (1979). "A microscopic theory for antiphase boundary motion and its applications to antiphase domain coarsening". Acta Metallurgica. 27: 1085–1095.

DiPrima, R.C. (2010). Elementary Differential Equations and Boundary Value Problems. Wiley. {{cite book}}: Cite has empty unknown parameter: |lnguage= (help)

Trefethen, L.N., ed.,; Embree, K., ed., The (Unfninished) PDE coffee table book {{citation}}: Cite has empty unknown parameter: |coauthors= (help)CS1 maint: extra punctuation (link)

Uecker, H. A short ad hoc introduction to spectral methods for parabolic PDE and the Navier-Stokes equations. {{cite book}}: Cite has empty unknown parameter: |coauthors= (help)

[1] Boyce 和 DiPrima (2010)

[2] Allen 和 Cahn (1979)

[3] 請注意，在程式設計時，您需要在每次計算下一個時間步長 $n+1$ 時更新非線性項 ( $u^{3}$ )。值得一提的是，對於每個時間步長，您需要從實空間轉換到傅立葉空間，再轉換回實空間，然後重複此過程。而對於熱方程，您可以在傅立葉空間中執行任意數量的時間步長，只有在記錄資料時才轉換回來。

[4] Uecker (2009)

[5] Trefethen 和 Embree

[1]

[2]

[3]

[4]

[5]