矩阵、随机化与分形图形

    Stetson大学的一个非常可爱的MM(以后本Blog将简称她为Stetson MM)和我分享了一个很神奇的东西。她们正在做一个线性代数的课题研究,题目的大致意思是“用矩阵来构造分形图形”。Stetson MM叫我试着做下面这个实验:对于一个坐标点(x,y),定义下面4个矩阵变换:
    
    然后,初始时令(x,y)等于(0,0),按照 T1 – 85%, T2 – 6%, T3 – 8%, T4 – 1% 的概率,随机选择一个变换对该点进行操作,生成的点就是新的(x,y);把它画在图上后,再重复刚才的操作,并一直这样做下去。我心里觉得奇怪,这为什么会得到分形图形呢?于是我写了一个简单的Mathematica程序:
list = {{0, 0}};
last = {{0}, {0}};
For[i = 0, i < 50000, i++, r = Random[];    If[r < 0.85, last = {{0.83, 0.03}, {-0.03, 0.86}}.last + {{0}, {1.5}},      If[r < 0.91, last = {{0.2, -0.25}, {0.21, 0.23}}.last + {{0}, {1.5}},        If[r < 0.99, last = {{-0.15, 0.27}, {0.25, 0.26}}.last + {{0}, {0.45}},          last = {{0, 0}, {0, 0.17}}.last + {{0}, {0}}        ]      ]    ];    list = Append[list, First[Transpose[last]]]; ] ListPlot[list, PlotStyle -> PointSize[0.002]]

    程序运行的结果真的是令我大吃一惊:竟然真的是一个分形图形!!我不禁再次对数学产生了一种崇敬和畏惧感!!

   


    静下心来仔细研究了一下,其实这个东西的道理很简单:这些变换的实质就是把图形转移为另一个位置下的另一个尺度,经过变换以后的每个点都仍然在这个图形内。不同的变换有着不同的作用。T1的作用很明显:把该点移动到下一片小树叶上相同的位置。我们用红色的线条标注了对某个点连续三次T1变换的路径。T2, T3的作用是,把这个点在整个大叶子上的位置“投射”到最底部的叶片上对应的位置,其中T2负责投射到左边,T3负责投影到右边。我们分别用蓝色箭头和绿色箭头来演示T2和T3的轨迹。比如,对大叶片的左边第三叶的中间某个点进行T2变换,得到左边第一叶的左边第三个更小的叶片的中间;如果再进行一次变换,则就变到了左边第一叶的左边第一叶。T3的作用也基本上类似,我就不再多说了。最后, T4的作用是把某个点初始化到(0, 0.17y),以后再经过一系列T1变换后就可以画出树叶中间的线条(的一部分)了,并等待某次T2或T3把该点变到最底层的叶片的中间线条上。有人会说,那这样变下去的话,岂不是所有点都只能在树叶中间的那根线条上?对!事实上,以后产生的每一个点都在某个“n级叶片”的中间那根线条上。

   

    其实,我们之前是见过类似的东西的。例如,Sierpinski三角形也有相似的结论:给出三角形的三个顶点,然后从其中一个顶点出发,每次随机向任意一个顶点移动1/2的距离(走到与那个顶点的连线的中点上),并在该位置作一个标记;无限次操作后所有的标记就组成了Sierpinski三角形。

30 条评论

  • Freeze

    4点……膜拜下
    问下mathematica语言matrix67是在哪里学的

  • Matrix78

    看来想做沙发一定要早啊

  • Matrix78

    MM不是很漂亮到是很有味道

  • Ai.Freedom

    这个MM的Blog好像见过

  • Matrix78

    Mathematica不是开源,我要赶紧删除.罚我我可负担不起

  • yangyanggoods

    只能说……我什么都没看懂……

  • pchu

    初三左右,我部文曲星就被人拿去做这类实验了~~

  • Ai.Freedom

    考虑买个Mathematica

  • dragon

    我用matlab验证了一下,并估计了高度的增长趋势。
    我觉得这就是生命–也许有一天,有机的生命形式会消失,但孕育于数学中的生命会一直存在下去。

  • cuiaoxiang

    你好,请问你的这段mathematica程序跑了多少时间?
    我用mathematica6.0 需要跑很长很长时间。。。
    不知道是不是哪里设置错了,我用C++写的代码很快就可以跑完了
    还是因为mathematica本身就这么慢?

  • collvey

    哦 我好像认识这个 这是著名的Barnsley’s fern吧

  • xiakexiake

    顾森我觉得你太扯了,去毛的北大中文系,我看不出你除了能泡到妞之外你在哪里还能学到什么。太扯太扯了,你整个大学四年貌似很充实,其实就是整个一个浪费整个一个空虚。我现在只想他妈的跟你强调老师的重要性。你现在整个不就是在做一个外力场下的随机行走么?走来走去能走成什么样?你有用过计算机模拟过吗?你不就是用数学来忽悠学物理的,用物理来忽悠学数学的,再用物理数学来忽悠所有其他人么?你在一片未知的大陆上遇见了不少风景,但你有本事去构建出全部美景的图像来吗?跟A说B,跟B说A,跟AB说CD,有意思么有意思吗?以前我觉得社会的蛀虫不过就那么几种,贪官奸商公务员,特么我现在发现你顾森居然又撑起了半边天。特么佩服佩服,特么我佩服个屁!

  • 燕仰

    楼上说什么北大中文系。

  • Pt-Cr

    15楼是不是太过激了……

  • spkobe

    回16楼:15楼是在说Matrix67呢,他是个文科生。
    不过15楼那娃确实过激了,没见过的就当补充下知识,当副业;好些人都身兼几专业的,学机电的,应用数学的搞编程去;等等。
    可能人觉得把喜欢的学科当主业学容易疲劳和失去激情吧;当副业反而能保持一直探索的热情

  • 莫莫

    回11楼:那是因为楼主的代码没有用上mathematica内置函数的缘故,由于每次Append时候都可能会分配额外空间,所以导致了慢。
    用mathematica内置函数NestList可以使得运行速度更快些:

    GetNextPoint[pt_] :=
    Module[{r, t1 = {{0.83, 0.03}, {-0.03, 0.86}}, p1 = {0, 1.5},
    t2 = {{0.2, -0.25}, {0.21, 0.23}}, p2 = {0, 1.5},
    t3 = {{-0.15, 0.27}, {0.25, 0.26}}, p3 = {0, 0.45},
    t4 = {{0, 0}, {0, 0.17}}, p4 = {0, 0}},
    r = Random[];
    If[r <= 0.85, t1.pt + p1,
    If[r <= 0.91, t2.pt + p2,
    If[r PointSize[0.002]]

  • 莫莫

    上述代码不知怎么的,上传过程中居然乱序了~,重来:

    GetNextPoint[pt_] :=
    Module[{r, t1 = {{0.83, 0.03}, {-0.03, 0.86}}, p1 = {0, 1.5},
    t2 = {{0.2, -0.25}, {0.21, 0.23}}, p2 = {0, 1.5},
    t3 = {{-0.15, 0.27}, {0.25, 0.26}}, p3 = {0, 0.45},
    t4 = {{0, 0}, {0, 0.17}}, p4 = {0, 0}},
    r = Random[];
    If[r <= 0.85, t1.pt + p1,
    If[r <= 0.91, t2.pt + p2,
    If[r PointSize[0.002]]

  • 莫莫

    这个网站,我服了!重新分两段发吧:
    GetNextPoint[pt_] :=
    Module[{r, t1 = {{0.83, 0.03}, {-0.03, 0.86}}, p1 = {0, 1.5},
    t2 = {{0.2, -0.25}, {0.21, 0.23}}, p2 = {0, 1.5},
    t3 = {{-0.15, 0.27}, {0.25, 0.26}}, p3 = {0, 0.45},
    t4 = {{0, 0}, {0, 0.17}}, p4 = {0, 0}},
    r = Random[];
    If[r <= 0.85, t1.pt + p1,
    If[r <= 0.91, t2.pt + p2,
    If[r <= 0.99, t3.pt + p3, t4.pt + p4]]]];

  • 莫莫

    这是调用上述函数50000次,得到一个列表,并画图:
    lst = NestList[GetNextPoint, {0, 0}, 100000];
    ListPlot[lst, PlotStyle -> PointSize[0.002]]
    经试验,速度快很多滴

  • esetlzn

    只要用append的循环肯定很卡••••••有空我也来试试!这个叶子我之前一直没画出来呵呵•••

  • chyanog

    这个应该比ls的Mathematica代码更快:

    v = {{{{0.85, 0.04}, {-0.04, 0.85}}, {0,
    1.6}}, {{{0.2, -0.26}, {0.23, 0.22}}, {0,
    1.6}}, {{{-0.15, 0.28}, {0.26, 0.24}}, {0,
    0.44}}, {{{0, 0}, {0, 0.16}}, {0, 0}}};
    g[p_] := #.p + #2 & @@ RandomChoice[{0.85, 0.07, 0.07, 0.01} -> v]
    Graphics[{Green, Point@NestList[g, {0, 0}, 10^4]}]

  • chyanog

    还可以更快:

    list = FoldList[#2.{#[[1]], #[[2]], 1} &, {0., 0.},
    RandomChoice[{85, 6, 8, 1}/100. ->
    N@{{{0.83, 0.03, 0}, {-0.03, 0.86, 1.5}}, {{0.2, -0.25,
    0}, {0.21, 0.23, 1.5}}, {{-0.15, 0.27, 0}, {0.25, 0.26,
    0.45}}, {{0, 0, 0}, {0, 0.17, 0}}}, 10^5]]; // AbsoluteTiming
    Graphics[{PointSize@Tiny, Point@list}]

  • Matrix89

    我是小号

  • myy

    我怎么运行程序的时候出现
    ??? for[i=0,i<50000,i++,r=random[];
    Error: Unbalanced or unexpected parenthesis or bracket.这一情况。怎么改

  • GZU_long

    考古啦,今天老师发的实验资料里面的附加题,打开网站直接开启时光机哈哈

回复给 莫莫 取消回复

6  +  4  =