几种图灵斑（Turing Patterns）的简单matlab演示（BZ反应、Gray-Scott模型、LE模型）

惯例声明：本人没有相关的工程应用经验，只是纯粹对相关算法感兴趣才写此博客。所以如果有错误，欢迎在评论区指正，不胜感激。本文主要关注于算法的实现，对于实际应用等问题本人没有任何经验，所以也不再涉及。

0 引言

1952年艾伦·图灵在他的论文中the chemical basisof morphogenesis（形态发生的化学基础）中，给出图灵斑图的大概概念，从数学和化学的角度，揭示了生物体表面斑纹的产生机理。

这种图灵斑图产生机理是一种反应-扩散体系，之后在多种生物、化学、地理等学科被观测到，可以说是一种遍布于大自然的普适性的规律。

本文以三种经典模型为例，介绍图灵斑图的matlab数值仿真方法。由于本人不是相关领域的研究人员，所以主要注重于数值计算方面，而不涉及图灵斑图相关的化学反应机理、非线性分岔等。

相关参考资料如下：
[1] Morphology of Experimental and Simulated Turing Patterns (2009)
[2] Simulating the Belousov-Zhabotinsky reaction (2017) https://scipython.com/blog/simulating-the-belousov-zhabotinsky-reaction/
[3] 一个演示Gray–Scott model的网址： http://www.karlsims.com/rdtool.html
建议用Chrome浏览器打开
[4]matlab区域填充_图灵斑图与反应扩散方程——Matlab的实现
https://blog.csdn.net/weixin_30278943/article/details/112067158
[5] The Lengyel–Epstein Reaction Diffusion System (2020)
[6] Machine learning with Patterns based on Lengyel-Epstein model
https://github.com/standing-o/Machine-learning_with_Patterns_based_on_Lengyel-Epstein_model
[7] 混乱博物馆-图灵斑图：生命图案的奥秘 https://zhuanlan.zhihu.com/p/29118927
[8] Reaction-Diffusion Model as a Framework for Understanding Biological Pattern Formation (2010)
[9] Gray-Scott-Complex Patterns in a Simple System (1993)
[10] 一类离散反应扩散捕食系统的分岔和斑图自组织研究(杨洪举)

列的东西杂七杂八的，毕竟也不是啥正经论文，大概按照重要性顺序列了一下。

下面举的图灵斑仿真例子，实际上就是求解某个数学方程的解。所以下面代码需要有简单的数值分析基础，只要知道如何简单离散解偏微分方程就可以。数值方法后文中也会稍微有所涉及，但不会太详细。

1 BZ震荡反应

BZ反应是在1958年由Belousov和Zhabotinski发现而得名。
模型可以简单被写为下面三个方程：

$$\begin{gathered} A+B\to 2A \\ B+C\to 2B \\ C+A\to 2C \end{gathered}$$

每个量下一时刻的值，会随上一时刻的变量而改变，可以被写作：
$$\begin{gathered} A_{t+1}=A_t+A_t(\alpha B_t-\gamma C_t) \\ {B}_{t+1}=B_t+B_t(\beta C_t-\alpha A_t) \\ {C}_{t+1}=C_t+C_t(\alpha A_t-\beta B_t) \end{gathered}$$

这里取常数$\alpha =1.2$，$\beta =1$，$\gamma =1$。

且为了模拟反应扩散，下一步计算时，还需要对该变量取周围8个变量取平均，作为当前网格点的值。在matlab里，用imfilter函数来实现。

初始值定义为0~1之前的随机数。

代码如下：

clear
clc
close all
% Belousov-Zhabotinsky反应
%构建网格
dt=0.5;%时间步长
N=500;%网格总数量
dx=1;%网格大小
x=dx*(1:N);
[X,Y]=meshgrid(x,x);
%初始化，采用随机初始化
A0=rand(size(X));
B0=rand(size(X));
C0=rand(size(X));

%方程初常数
alpha=1.2;
beta=1;
gamma=1;

%微分方程求解
F=ones(3)/8;F(2,2)=0;
A_Old=A0;
B_Old=B0;
C_Old=C0;
for k_t=1:600
    %模拟扩散项
    A_Old=imfilter(A_Old,F,'circular');
    B_Old=imfilter(B_Old,F,'circular');
    C_Old=imfilter(C_Old,F,'circular');

    %1阶时间精度
    A_New=A_Old+A_Old.*(alpha*B_Old-gamma*C_Old)*dt;
    B_New=B_Old+B_Old.*(beta*C_Old -alpha*A_Old)*dt;
    C_New=C_Old+C_Old.*(gamma*A_Old- beta*B_Old)*dt;

    A_Old=A_New;
    B_Old=B_New;
    C_Old=C_New;
	%绘图
    if k_t>50 && mod(k_t,2)==1
        figure(1)
        clf
        pcolor(X,Y,A_New);shading interp
        caxis([0,1])
        pause(0.01)
    end
end

可以看到BZ震荡反应典型的两个图案：震荡和螺旋。

2 Gray Scott模型

Gray Scott模型是一种典型的反应扩散系统。其方程可以写作：
$$\begin{gathered} U+2V\to 3V \\ V\to P \end{gathered}$$
用来计算的微分方程可以写作：
∂ U ∂ t = D U Δ U − U V 2 + F ( 1 − U )   ∂ V ∂ t = D V Δ V + U V 2 − ( F + k )