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}$

剩下的就是對所有時間步長完成後計算的解進行逆傅立葉變換，將其轉換回實空間。這是一個線性偏微分方程，因此只需要在最後進行一次逆傅立葉變換。我們將在後面看到，對於非線性偏微分方程，情況有所不同。該方法的 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$ for $x<\pi$ 和 $u(t=0,x)=1$ for $x\geq \pi$
直到時間 $T=1$ 。您可以透過使用變數分離或假設解的形式為 $u(x,t)=f(x+t)$ 並推斷出 $f$ 是什麼，以滿足初始條件來實現。在這兩種情況下，請使用前向尤拉、後向尤拉和 Crank-Nicolson 時間步長方案。在計算完每種情況下的精確解後，檢查這三種方法中最終時間的最大誤差如何取決於時間步長。

非線性方程

一維 Allen-Cahn 方程

到目前為止，我們只處理了線性方程。現在是時候進行非線性 PDE 了。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)$ 並繪製由列表中Matlab程式碼計算的空間時間演化 C 。計算結果見圖 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)
↑ 請注意，在程式設計時，您需要更新非線性項 ( $u^{3}$ ) 每次您想要計算下一個時間步長 $n+1$ 。值得一提的原因是，對於每個時間步長，您都需要從實空間轉換到傅立葉空間，然後再轉換回實空間，然後重複此過程。而對於熱方程，您可以在傅立葉空間中執行任意數量的時間步長，只需在記錄資料時再轉換回實空間即可。
↑ Uecker (2009)
↑ Trefethen and Embree

參考文獻

Allen, S.M.; Cahn, J.W. (1979). "反相界運動的微觀理論及其在反相疇粗化中的應用". Acta Metallurgica. 27: 1085–1095.

Boyce, W.E.; DiPrima, R.C. (2010). 常微分方程及其邊值問題. Wiley.

Trefethen, L.N., ed.,; Embree, K., ed., (未完) PDE 咖啡桌書{{citation}}: CS1 maint: extra punctuation (link)

Uecker, H. 關於拋物線型偏微分方程和 Navier-Stokes 方程的譜方法簡明 ad hoc 介紹.

[1] Boyce 和 DiPrima (2010)

[2] Allen 和 Cahn (1979)

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

[4] Uecker (2009)

[5] Trefethen and Embree

[1]

[2]

[3]

[4]

[5]