几种图灵斑（Turing Patterns）的简单matlab演示（BZ反应、Gray-Scott模型、LE模型）
0 引言
1 BZ震荡反应
2 Gray Scott模型
3 LE模型（CIMA反应）
惯例声明：本人没有相关的工程应用经验，只是纯粹对相关算法感兴趣才写此博客。所以如果有错误，欢迎在评论区指正，不胜感激。本文主要关注于算法的实现，对于实际应用等问题本人没有任何经验，所以也不再涉及。
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发现而得名。 模型可以简单被写为下面三个方程：
A + B → 2 A B + C → 2 B C + A → 2 C A+B\to2A \\ B+C\to2B \\ C+A\to2C \\ A+B→2AB+C→2BC+A→2C
每个量下一时刻的值，会随上一时刻的变量而改变，可以被写作： A t + 1 = A t + A t ( α B t − γ C t ) B t + 1 = B t + B t ( β C t − α A t ) C t + 1 = C t + C t ( α A t − β B t ) 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) \\ At+1=At+At(αBt−γCt)Bt+1=Bt+Bt(βCt−αAt)Ct+1=Ct+Ct(αAt−βBt)
这里取常数 α = 1.2 \alpha=1.2 α=1.2， β = 1 \beta=1 β=1， γ = 1 \gamma=1 γ=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模型是一种典型的反应扩散系统。其方程可以写作： U + 2 V → 3 V V → P U+2V \to 3V \\ V \to P U+2V→3VV→P 用来计算的微分方程可以写作： ∂ U ∂ t = D U Δ U − U V 2 + F ( 1 − U ) ∂ V ∂ t = D V Δ V + U V 2 − ( F + k )