傅里叶谱方法-傅里叶谱方法求解二维 Schnakenberg 模型及其Matlab程序实现

3.3.4 二维 Schnakenberg 模型

Schnakenberg模型是一种描述化学反应动力学的数学模型，旨在研究化学反应中自组织现象和空间模式的形成。

该模型最早由德国化学家Heinrich Otto Wieland提出，后由德国数学家Theodor Schnakenberg在其博士论文中进行了推导和解析，并得到了该模型的名字。

Schnakenberg模型的意义在于，它可以模拟一系列反应中出现的自组织现象，例如化学波、化学斑图、化学分叉等。这些现象是由于反应物浓度分布在空间上的不均匀性导致的，通过Schnakenberg模型可以更好地理解这些现象，并为实验提供一定的指导和验证。

此外，Schnakenberg模型还可以应用于其他领域，例如生物学、物理学、生态学等。在这些领域中，该模型可以用于研究许多自组织现象，例如细胞分裂、病毒传播、生态系统演化等。

这个模型描述的是两种化学物质之间的反应过程。具体来说，化学物质 $u$ 会被消耗掉，同时生成化学物质$v$。反应速率的大小取决于$u$ 和$v$的浓度，而扩散系数则描述了物质在空间中的扩散速度。

Schnakenberg模型有许多有趣的特性和行为。例如，在某些参数范围内，这个模型的解可以显示出空间上的分布和时间上的演化，形成有序的斑图（patterns）。这些斑图包括静态和动态的形态，它们的形状和数量取决于模型参数的取值。这种有序结构的出现和演化，可以通过Schnakenberg模型来解释和理解实验室中观察到的许多化学反应和生物学现象。具体来说，当模型的参数取值适当时，Schnakenberg模型的解会形成一个包含不同浓度级别的周期性结构。这种周期性结构被称为Turing结构，是由Alan Turing在1952年提出的一种理论模型。Turing结构是许多自组织系统中的一种普遍现象，例如斑马条纹、蝴蝶翅膀上的斑点等。
在Schnakenberg模型中，Turing结构的出现是由反应和扩散之间的相互作用所导致的。具体来说，当扩散系数足够大时，化学物质可以在空间上扩散，使得浓度分布均匀。但当反应速率较慢时，这些化学物质在局部区域内聚集起来，并且在一定的空间尺度上发生反应。这个过程会导致局部区域内的浓度发生变化，从而形成周期性的结构。
需要注意的是，Turing结构的出现是非线性动力学的一个重要特征。这种结构的产生和演化需要考虑多个因素之间的相互作用，因此通常需要使用数值计算和数学分析的方法来研究其性质和行为。

需要注意的是，Schnakenberg模型是一种简化的模型，它假设了化学物质的浓度在空间和时间上均匀分布，没有考虑一些具体的生物或化学体系的细节和复杂性。因此，这个模型只能作为一种基本的参考模型，而不是具体体系的精确描述。

近年来，Schnakenberg模型被广泛应用于不同领域的研究中，包括流体力学、生物学、材料科学等。这个模型的理论和数值分析方法也在不断地发展和完善，为我们深入理解和探究复杂动态系统的性质和行为提供了有力的工具和平台。

斑图（pattern）是一类普遍存在于自然界、在时间或空间上具有某种规律的非均匀宏观结构。反应-扩散系统 (reaction-diffusion system) 是斑图理论中研究得最为广 泛的系统, 它起源于化学反应系统, 但又不局限于此, 还广泛应用于生物学、物理学、医学、金融学等。Schnakenberg 模型是反应-扩散系统中的一个有趣的模型, 数学形式如下:
{∂u∂t=(∂2∂x2+∂2∂y2)u+γ(a−u+u2v)∂v∂t=d(∂2∂x2+∂2∂y2)v+γ(b−u2v)(3-38)\\left\\{\\begin{array}{l} \\frac{\\partial u}{\\partial t}=\\left(\\frac{\\partial^2}{\\partial x^2}+\\frac{\\partial^2}{\\partial y^2}\\right) u+\\gamma\\left(a-u+u^2 v\\right) \\\\ \\frac{\\partial v}{\\partial t}=d\\left(\\frac{\\partial^2}{\\partial x^2}+\\frac{\\partial^2}{\\partial y^2}\\right) v+\\gamma\\left(b-u^2 v\\right) \\end{array}\\right.\\tag{3-38} ⎩⎨⎧​∂t∂u​=(∂x2∂2​+∂y2∂2​)u+γ(a−u+u2v)∂t∂v​=d(∂x2∂2​+∂y2∂2​)v+γ(b−u2v)​(3-38)
其中, $u$ 和 $v$ 可看做两种化学反应物质的浓度, $x、y$ 为空间坐标, $t$ 为时间, $a、b、d$ 、 $\gamma$ 为常数。对式 (3-38) 做二维傅里叶变换, 可将其转化为偏微分方程组:

{∂u^∂t=−(kx2+ky2)u^^+γ⋅F{a−F−1[u^]+F−1[u^]2F−1[v^]}∂v^∂t=−d(kx2+ky2)v^+γ⋅F{b−F−1[u^]2F−1[v^]}(3-39)\\left\\{\\begin{array}{l} \\frac{\\partial \\hat{u}}{\\partial t}=-\\left(k_x^2+k_y^2\\right) \\hat{\\hat{u}}+\\gamma \\cdot F\\left\\{a-F^{-1}[\\hat{u}]+F^{-1}[\\hat{u}]^2 F^{-1}[\\hat{v}]\\right\\} \\\\ \\frac{\\partial \\hat{v}}{\\partial t}=-d\\left(k_x^2+k_y^2\\right) \\hat{v}+\\gamma \\cdot F\\left\\{b-F^{-1}[\\hat{u}]^2 F^{-1}[\\hat{v}]\\right\\} \\end{array}\\right.\\tag{3-39} {∂t∂u^​=−(kx2​+ky2​)u^^+γ⋅F{a−F−1[u^]+F−1[u^]2F−1[v^]}∂t∂v^​=−d(kx2​+ky2​)v^+γ⋅F{b−F−1[u^]2F−1[v^]}​(3-39)

参数取值为: $a=0.1,b=0.8,d=26,\gamma =100$ 。为了得到靶型波, 将初始条件设置为: $u$ 在 $x-y$ 平面原点处为 1 , 在其他位置为 $0,v$ 在整个 $x-y$ 平面上均为 1 。利用傅里叶谱方法 求解该模型的代码如下。

主程序代码

clear all; close all;
L=16; N=64;
kx=2*pi/L*[0:N/2-1 -N/2:-1]; ky=kx;
[kX,kY]=meshgrid(kx,ky);
K2=kX.^2+kY.^2;
%初始条件
u=zeros(N); u(N/2,N/2)=1; v=ones(N);
ut=fft2(u); vt=fft2(v);
uvt=[ut(:); vt(:)];
%求解
a=0.1; b=0.8; d=26; gamma=100; t=[0:0.1:0.3];
[t,uvtsol]=ode45('schnakenberg',t,uvt,[],K2,N,gamma,a,b,d);
%画图
for n=1:4subplot(2,2,n)gca=pcolor(ifft2(reshape(uvtsol(n,1:N^2),N,N))); axis offset(gca,'LineStyle','none'), shading interptitle(['t=' num2str(t(n))]), axis('square'),
%     colormap('gray')
end

程序输出的结果如图所示，这与化学实验中观察到的靶型波一致。

Matlab 将计算结果制作成 gif 动画

引言

求解包含时间的偏微分方程 (组) 将得到随着时间变化的数值结果, 把这样的数据制作成 gif 动画并结合到幻灯片中, 在毕业答辩、小组讨论、课堂教学等场合有着广泛的应用。生动的彩色 gif 动画具有很强的表现力, 令人刮目相看, 大大提高了报告人所讲述理论结果的直观性、生动性、观赏性。

函数介绍

生成 gif 动画主要用到 $4$ 个函数: getframe、frame2im、rgb2ind、imwrite。

1. getframe 函数的一般调用形式为: F=getframe(h), 其作用是截取句柄为 h 的窗口内的一帧图像。

2. frame2im 函数的作用是把一帧截图转为图像数据。

3. rgb2ind 函数的作用是将 RGB 图像转换为索引图像, 一般调用形式为: [X, map]= rgb2ind(RGB,n) 。其中, X、map 分别为转换后的图像数据和颜色表数据, RGB 为转换前的图像数据, n 指定 map 中的颜色数。

4. imwrite 函数的作用是将图像数据写入图像文件, 一般调用形式为: imwrite(X,map,filename,fmt,Param1,Val1,Param2,Val2...)。其中, X、map 意义同上, filename 为文件名, fmt 为文件格式,Param1,Val1,Param2,Val2... 为若干可选参数及其取值。如：参数 LoopCount 为动画的循环播放次数, 这里设为 inf, 即无穷大。参数 DelayTime 为每帧间隔时间, 单位秒。参数 WriteMode 为写入文件的模式, 有覆盖 overwrite (默认) 和追加 append 两种选择。

Matlab源代码

生成 gif 动画的示例代码如下:

clear
clc
close
x=-1:0.02:1;y=x;
[X,Y]=meshgrid(x,y);
filename='test.gif';
for a=1:10u=a*exp(-10*(X.^2+Y.^2));mesh(x,y,u),axis([-1 1 -1 1 0 10]),drawnowim=frame2im(getframe(gcf));[A,map]=rgb2ind(im,256);if a==1%先以覆盖模式写入指定的gif文件imwrite(A,map,filename,'gif','LoopCount',Inf,'DelayTime',0.1);else%再以追加模式将每一帧写入gif文件imwrite(A,map,filename,'gif','WriteMode','append','DelayTime',0.1);end
end

运行代码之后在当前目录下生成 gif 文件, 该动画显示了一个三维高斯函数的峰值逐渐增大的过程。

为了更加直观的感受求解 Schnakenberg 模型得到的靶型波，此处结合上述操作将计算结果制作成如下 gif 动画。

修改后的程序源代码：

clear all; close all;
L=16; N=64;
kx=2*pi/L*[0:N/2-1 -N/2:-1]; ky=kx;
[kX,kY]=meshgrid(kx,ky);
K2=kX.^2+kY.^2;
%初始条件
u=zeros(N); u(N/2,N/2)=1; v=ones(N);
ut=fft2(u); vt=fft2(v);
uvt=[ut(:); vt(:)];
%求解
a=0.1; b=0.8; d=26; gamma=100; t=[0:0.01:1];
[t,uvtsol]=ode45('schnakenberg',t,uvt,[],K2,N,gamma,a,b,d);
%画图
filename = 'output.gif'; % 设置保存的gif文件名
for n = 1:length(t)gca=pcolor(ifft2(reshape(uvtsol(n,1:N^2),N,N))); axis offset(gca,'LineStyle','none'), shading interptitle(['t=' num2str(t(n))]), axis('square'),% colormap('gray')% 捕捉当前子图并转换为gif图像帧frame = getframe(gcf);im=frame2im(getframe(gcf));[A,map]=rgb2ind(im,256);% 写入gif图像文件if n == 1%先以覆盖模式写入指定的gif文件imwrite(A,map,filename,'gif','LoopCount',Inf,'DelayTime',0.01);else%再以追加模式将每一帧写入gif文件imwrite(A,map,filename,'gif','WriteMode','append','DelayTime',0.01);end
end