首先放上实物图和模拟结果图：
一、蒙特卡罗方法
蒙特卡罗方法（Monte Carlo method），也称统计模拟方法，是二十世纪四十年代中期由于科学技术的发展和电子计算机的发明，而被提出的一种以概率统计理论为指导的一类非常重要的数值计算方法。是指使用随机数（或更常见的伪随机数）来解决很多计算问题的方法。
基本思想
当所求解问题是某种随机事件出现的概率，或者是某个随机变量的期望值时，通过某种"实验"的方法，以这种事件出现的频率估计这一随机事件的概率，或者得到这个随机变量的某些数字特征，并将其作为问题的解。
工作过程
蒙特卡罗方法的解题过程可以归结为三个主要步骤:构造或描述概率过程;实现从已知概率分布抽样;建立各种估计量。
蒙特卡罗方法解题过程的三个主要步骤:
(1)构造或描述概率过程
对于本身就具有随机性质的问题，如粒子输运问题，主要是正确描述和模拟这个概率过 程，对于本来不是随机性质的确定性问题，比如计算定积分，就必须事先构造一个人为的概率过程，它的某些参量正好是所要求问题的解。即要将不具有随机性质的问题转化为随机性质的问题。
(2)实现从已知概率分布抽样
构造了概率模型以后，由于各种概率模型都可以看作是由各种各样的概率分布构成的，因此产生已知概率分布的随机变量(或随机向量)，就成为实现蒙特卡罗方法模拟实验的基本手段，这也是蒙特卡罗方法被称为随机抽样的原因。随机数是我们实现蒙特卡罗模拟的基本工具。
(3)建立各种估计量
一般说来，构造了概率模型并能从中抽样后，即实现模拟实验后，我们就要确定一个随机变量，作为所要求的问题的解，我们称它为无偏估计。建立各种估计量，相当于对模拟实验的结果进行考察和登记，从中得到问题的解。
二、皮蛋花纹的形成
松花晶体为氢氧化镁水合晶体[Mg(OH)2].镁离子的来源除蛋内含有少量外，主要来自蛋壳和蛋膜。在浸泡松花蛋的过程中，壳和膜上的Mg离子受NaOH、NaCl及H2O等的作用，部分被溶解并随之进入蛋内。当蛋内Mg离子浓度达到足以同OH-离子化合形成大量Mg(OH)2时，由于是处于蛋白质凝胶体内的特殊环境下，它们就形成水合晶体。
由于布朗运动，晶体颗粒受到无规则的碰撞作用，可以看成一个随机过程，晶体颗粒在这种随机过程中靠近结合，最终形成了晶体花纹。使用蒙特卡罗方法对皮蛋花纹的模拟就是基于这个模型。为了简化，首先考虑核心的随机过程，而将物理的状态变化进行简化。
简化模型：给定一个晶核，晶核附近区域的晶体粒子随机运动，在很接近晶核时，结合到一起，晶核持续生长变成晶花。
过程如下图所示：
三、程序实现
设计主要分为三个部分：
（1）在晶核周围的区域产生随机分布的晶体粒子；
（2）这些随机分布的晶体粒子进行随机运动；
（3）粒子和晶核的相遇结合以及结合之后新晶核的更新。
下面给出最简单的模拟程序（结果已经在文章开头给出）：
'''----------------------------------------------------------------------------------------------------------
# Name: Diffusion-limited Aggregation(DLA) simulation
# Author: Wenchao Liu
# Date: 30/07/2017
# Description: To simulate the diffusion limited aggrehation(DLA) to the particles
# by using Mento Carlo method ( a random method).
# Debug in Python 3.5.3.
# This program was modified in 02/08/2017.
-----------------------------------------------------------------------------------------------------------'''
from math import sqrt
from random import randint
import matplotlib.pyplot as plt
# 定义随机行走函数
def Random_walk():
# 产生随机数 0,1,2,3，分别代表随机往东，西，南，北走出一步
walk = randint(0,3)
return walk
# 定义一个产生随机运动粒子的效应器
def Create_p(i, r,x_core, y_core):
xn = 0
yn = 0
# 随机产生粒子的初始坐标
x = randint(-r, r)
y = randint(-r, r)
while (sqrt(x ** 2 + y ** 2) < r):
walk = Random_walk()
if walk == 0:
x = x + 1
elif walk == 1:
x = x - 1
elif walk == 2:
y = y - 1
else:
y = y + 1
for j in range(i):
if (abs(sqrt((x - x_core[j]) ** 2 + (y - y_core[j]) ** 2) == 1)):
xn = x
yn = y
break
return xn, yn
# 粒子聚集的核
N = int(input('请输入模拟的次数： '))
x_core = [0] * N
y_core = [0] * N
for i in range(1,N):
'''模拟初期阶段，核比较小，粒子与核相遇需要的时间太长，
这里做了一个逐步的近似，加快这一过程'''
if i < 50 or i == 50:
r = 10 # 产生粒子的区域半径
elif 50 < i < 180 or i == 180:
r = 20
elif 180 < i < 650 or i == 650:
r = 30
elif 650 < i < 1000 or i == 1000:
r = 40
elif 1000 < i <2000 or i == 2000:
r = 50
elif 2000 < i < 3500 or i == 3500:
r = 60
elif 3500 < i < 5000 or i == 5000:
r = 70
elif 5000 < i < 10000 or i == 10000:
r = 80
else:
r = 150
x_core[i], y_core[i] = Create_p(i, r,x_core, y_core)
#break
#print(x_core, y_core)
plt.figure()
#plt.plot(xn, yn, 'b-', xn, yn,'g.', xn[0], yn[0], 'ro', xn[len(xn)-1], yn[len(xn)-1], 'mo')
plt.plot(x_core, y_core, 'm+',x_core[0], y_core[0], 'r*')
plt.axis([-r-20,r+20,-r-20,r+20])
plt.show()
对这个模型的实现自然有许多可以改进的地方。一个是模拟过程中的计算效率问题，主要是粒子与晶核相遇结合的这个过程的算法改进；另一个是考虑物理状态在这些过程中的变化。
