Matlab 中的并行谱数值方法/示例

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 等文献中找到）。其中 $n=64$ 是使用列表 B 中的 Matlab 程序获得的。

( i)

热方程的数值解，公式 1 使用后向欧拉法计算。

( 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()

练习

1. 编写程序使用 Crank-Nicolson 方法求解热方程。

2. 求解对流方程 $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 提出的，其表达式为

${\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 方程，我们发现

${\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

( ii)

使用隐式显式方法计算的一维 Allen-Cahn 方程（方程 2 ）的数值解，其中 $\epsilon =0.001$ 且 $u(x,t=0)=0.25\sin(x)$。

二维 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)

2D Allen-Cahn 方程的数值解，即方程 4 在时间 $t=500$ 时，且 $\epsilon =0.1$ 以及 $u(x,y,t=0)=\sin(2\pi x)+0.001\cos(16\pi x)$，使用隐式显式方法计算得出。

用于生成图 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。有关这些方程的另一个入门资料来源是 Trefethen 和 Embree。

1. 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$ 取为较小的值，请确保使用足够的网格点以获得正确的数值解。检查这一点的一种简单方法是不断增加网格点的数量并检查解是否没有变化。另一种检查方法是计算傅里叶系数并检查最高的系数是否衰减到机器精度。
2. 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}}$ 其中 $u$ 具有周期边界条件。
   • 该方程式模拟了什么，你预计它的解会表现出什么样的行为？
   • 使用隐式-显式方法求解该方程式的数值解。
3. 一维 Gray-Scott 方程组如下：${\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$ 其中 $d_{1}$，$d_{2}$，$f$ 和 $k$ 为常数。
   • 该方程式模拟了什么，你预计它的解会表现出什么样的行为？
   • 使用隐式-显式方法求解该方程式的数值解。尝试几个不同的 $d_{1}$，$d_{2}$，$f$ 和 $k$ 的值，并将所得模式与文献中的发现进行比较。
4. 二维 Swift-Hohenberg 方程如下：${\frac {\partial u}{\partial t}}=-\Delta ^{2}u+2\Delta u+(\alpha -1)u-u^{3},$
   • 该方程式模拟了什么，你预计它的解会表现出什么样的行为？
   • 使用隐式-显式方法求解该方程式的数值解，并尝试几个不同的 $\alpha$ 值。
5. 二维 Gray-Scott 方程组如下：${\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$ 其中 $d_{1}$，$d_{2}$，$f$ 和 $k$ 为常数。
   • 该方程式模拟了什么，你预计它的解会表现出什么样的行为？
   • 使用隐式-显式方法求解该方程式的数值解。
6. 二维复 Ginzburg-Landau 方程表示如下：${\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$ 的多个值，找到该方程的数值解。

参考文献

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)