假设我们在做一个化学反应实验，存在U和V两种化学物质，其中每单位U物质在两单位V物质的催化下，生成新的V物质。那么V物质的浓度将会发生怎样的变化？
设想一下，随着时间的推移，反应皿中 V物质的浓度应该是持续升高的，由于现实情况，比如气压、重力、气温等因素，V物质在容器中并不会均匀的分散在溶液中，而是“抱团分布”。
下面的视频就是0到5000秒，化学物质V的分布变化图，这种图一般会展现出斑点、条纹、迷宫等样式，所以叫做斑图。斑图的模拟研究起源于图灵，所以又称其图灵斑图。
回到V物质的斑图上，视频中的黄色区域即V物质浓度高的区域，起始时，黄色区域为一个“环状”，随着时间的推移，黄色区域逐渐扩散，最终成为了一个“迷宫”。
那么其他色彩区域是否指的是没有V物质的存在呢，并不是这样的，其他区域的颜色变化其实是指V物质的浓度不高，本视频中，“黄色——浅黄色——浅蓝色——蓝色”区域的V物质的浓度是逐渐降低的。
下面将详细解释图灵斑图的实现。
01
图灵斑图是什么？
图灵于1952年在他的这篇论文中，创造性地使用反应——扩散系统的数学模型来描述自然界中的斑图，比如虎纹，豹斑。这类可以用反应-扩散方程描述的斑图，被后人称为图灵斑图（Turing Pattern）。
与此同时，斑图不仅在动物的皮肤上展现，同样常见于半干旱地区的植被分布模拟，以及生物种群（捕食者和被捕食者）的分布上等等。
由于斑图是一种研究区域变化随时间变化的一种方法，其可以用于数学、统计学、生物学、地理学、物理学、化学、医学等各个专业。
老虎斑纹的出现究其根本是化学物质浓度的区域性差异所导致的，所以本文以Gray-Scott模型为例，讲述反应扩散方程的斑图的模拟与绘制。
02
反应扩散方程
下列方程组为反应扩散方程的一般形式：
其中，和为方程组中的扩散项，DU和DV分别为U和V的扩散系数，而f（U,V）与g（U,V）分别为U和V的生成率，图灵认为这两项为二次多项式。
上述内容比较抽象，下面直接以Gray-Scott模型为例进行讲述。
03
Gray-Scott模型
化学反应方程：
U和V为两种化学物质，通过反应产生新的V。因为V出现在第一反应的两边，所以它作为催化剂，对自身的生产起着催化剂的作用。而P为一种惰性产物。
该体系的总体行为由下面的公式描述，两个方程描述了两种化学物质的三种增减源。
上述方程组中：
F——进料率，代表补给率，常数；
k——去除率，常数；
04
斑图形成的条件
想要实现斑图就需要考虑其时空性质，即时间和空间的转变导致斑图的形成。
时空离散的捕食者系统斑图形成的条件具体来说为以下三点：
（1）稳定空间均匀定态的条件：
离散系统存在一个非平凡空间均匀定态，此定态对空间均匀扰动是稳定的；
（2） Neimark- Sacker失稳的条件：
稳定空间均匀定态在 Neimark- Sacker分岔作用下而变得不稳定；
（3）图灵失稳的条件：
稳定空间均匀定态对于至少一类空间非均匀扰动变得不稳定。
受篇幅所限，此处不讲述该模型的数值分析内容，该模型是符合这三点要求的。我们直接进行建模绘图。
05
斑图的形成
思考：
在绘制斑图时，我们可以将一张斑图划分为N*N个正方形格子，填充每个格子颜色，就可以绘制出一个斑图。而填充颜色是可以考虑通过建立一个N*N的矩阵，通过比较矩阵的每个元素与其他元素的值的大小可以确定颜色填充的差异。
因此，我们需要计算每个格子（矩阵元素）的值，此时最大的问题是如何确定边界上格子的值，比如矩阵第一列的元素的值，此时我们可以通过将这个矩阵视作无缝连接的矩阵，即矩阵第一行与最后一行相邻，第一列与最后一列相邻。
在matlab中可以通过circshift（）函数完成我们的构想。
Step 1：构造网格影响方程函数：
function out = my_laplacian(in)
out = -in ...
+ .20*(circshift(in,[ 1, 0]) + circshift(in,[-1, 0]) ...
+ circshift(in,[ 0, 1]) + circshift(in,[ 0,-1])) ...
+ .05*(circshift(in,[ 1, 1]) + circshift(in,[-1, 1]) ...
+ circshift(in,[-1,-1]) + circshift(in,[ 1,-1]));
Step 2：设置系数值与初始矩阵
f=.055;%进料率
k=.062;%去除率
da = 1;%U的扩散率
db = .5;%V的扩散率
% 网格的大小
width = 128;
% 5,000个模拟秒，每模拟秒4步
dt = .25;
stoptime = 5000;
设置初始矩阵：
function [t, A, B] = initial_conditions(n)
t = 0;
A = ones(n);
B = zeros(n);
B(51:60 ,51:70) = 1;
B(61:80,71:80) = 1;
Step 3：运行主程序
[t, A, B] = initial_conditions(width);
%主程序
tic
nframes = 1;
while t
anew = A + (da*my_laplacian(A) - A.*B.^2 + f*(1-A))*dt;
bnew = B + (db*my_laplacian(B) + A.*B.^2 - (k+f)*B)*dt;
A = anew;
B = bnew;
t = t+dt;
nframes = nframes+1;
end
%画图
axes('Position',[0 0 1 1])
axis off
hi = image(B);
hi.CDataMapping = 'scaled';
delta = toc;
disp([num2str(nframes) ' frames in ' num2str(delta) ' seconds']);
下面第一幅图为B（对应前文的V）的区域分布斑图
Step 4：拓展
如果我们想观察一段时间内B的分布区域的变化情况，那么我们可以通过绘制动图来实现。
[t, A, B] = initial_conditions(width);
targetframerate = 24;
frametime = 1/(24*60*60*targetframerate);
nextframe = now + frametime;
tic
nframes = 1;
while t
anew = A + (da*my_laplacian(A) - A.*B.^2 + f*(1-A))*dt;
bnew = B + (db*my_laplacian(B) + A.*B.^2 - (k+f)*B)*dt;
A = anew;
B = bnew;
hi.CData = B;
t = t+dt;
ht.String = ['Time = ' num2str(t)];
if now > nextframe
drawnow
nextframe = now + frametime;
end
nframes = nframes+1;
end
delta = toc;
disp([num2str(nframes) ' frames in ' num2str(delta) ' seconds']);
在原来的基础上，通过增加绘制动图的指令即可实现。
