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]