求知 文章 文库 Lib 视频 iPerson 课程 认证 咨询 工具 讲座 Model Center   Code  
会员   
要资料
 
 

利用python进行数据分析
搭建python环境
统计 1880-2010 年间 全美婴儿姓名的趋势
Ipython & Ipython Notebook
NumPy Basics: Arrays and Vectorized Computation
Pandas(Python Data Analysis Library)
数据加载、储存与文件格式
绘图和可视化(Matplotlib)
时间序列
经济,金融数据应用
补充例子
国王与囚徒
利用python进行科学计算
分形与混沌之-Mandelbrot集合
分形与混沌之-迭代函数系统(IFS)
分形与混沌之-蔡氏电路模拟
对于μ子实验数据进行快速处理与分析
37%法则 - "非诚勿扰" 中的博弈
关于时间/日期的转换
深入研究
一切从游戏开始:完整的一个 python to hack 实例!
习题:扑克牌发牌
 
 

分形与混沌之-迭代函数系统(IFS)
1897 次浏览
14次  

分形与混沌之-迭代函数系统(IFS)

迭代函数系统是一种用来创建分形图案的算法,它所创建的分形图永远是绝对自相似的。下面我们直 接通过绘制一种蕨类植物的叶子来说明迭代函数系统的算法:

有下面4个线性函数将二维平面上的坐标进行线性映射变换:

x(n+1)= 0; y(n+1) = 0.16 * y(n)

x(n+1) = 0.2 * x(n) ? 0.26 * y(n); y(n+1) = 0.23 * x(n) + 0.22 * y(n) + 1.6

x(n+1) = ?0.15 * x(n) + 0.28 * y(n); y(n+1) = 0.26 * x(n) + 0.24 * y(n) + 0.44

x(n+1) = 0.85 * x(n) + 0.04 * y(n); y(n+1) = ?0.04 * x(n) + 0.85 * y(n) + 1.6

所谓迭代函数是指将函数的输出再次当作输入进行迭代计算,因此上面公式都是通过坐标 x(n),y(n) 计 算变换后的坐标 x(n+1),y(n+1)。

现在的问题是有4个迭代函数,迭代时选择哪个函数进行计算呢?

我们为每个函数指定一个概率值,它们依次为1%, 7%, 7%和85%。选择迭代函数时使用通过每个函数的 概率随机选择一个函数进行迭代。

上面的例子中,第四个函数被选择迭代的概率最高。

最后我们从坐标原点(0,0)开始迭代,将每次迭代所得到的坐标绘制成图,就得到了叶子的分形图案。 下面的程序演示这一计算过程:

# -*- coding: utf-8 -*-
# by whyx 2016/6

import numpy as np
import pylab as pl
from math import log
import time
from matplotlib import cm

# 蕨类植物叶子的迭代函数和其概率值
eq1 = np.array([[0,0,0],[0,0.16,0]])
p1 = 0.01

eq2 = np.array([[0.2,-0.26,0],[0.23,0.22,1.6]])
p2 = 0.07

eq3 = np.array([[-0.15, 0.28, 0],[0.26,0.24,0.44]])
p3 = 0.07

eq4 = np.array([[0.85, 0.04, 0],[-0.04, 0.85, 1.6]])
p4 = 0.85


def ifs(p, eq, init, n):
# 进行函数迭代
# p: 每个函数的选择概率列表
# eq: 迭代函数列表
# init: 迭代初始点
# n: 迭代次数
#返回值: 每次迭代所得的X坐标数组, Y坐标数组, 计算所用的函数下标

# 迭代向量的初始化
pos = np.ones(3, dtype=np.float)
pos[:2] = init

# 通过函数概率,计算函数的选择序列
p = np.add.accumulate(p)
rands = np.random.rand(n)
select = np.ones(n, dtype=np.int)*(n-1)
for i, x in enumerate(p[::-1]):
select[rands<x] = len(p)-i-1

# 结果的初始化
result = np.zeros((n,2), dtype=np.float)
c = np.zeros(n, dtype=np.float)

for i in xrange(n):
eqidx = select[i] # 所选的函数下标
tmp = np.dot(eq[eqidx], pos) # 进行迭代
pos[:2] = tmp # 更新迭代向量

# 保存结果
result[i] = tmp
c[i] = eqidx

return result[:,0], result[:, 1], c

start = time.clock()

p1, p2, p3, p4 = 0.01, 0.07, 0.07, 0.85
#p1, p2, p3, p4 = 0.11, 0.17, 0.17, 0.55

x, y, c = ifs([p1,p2,p3,p4],[eq1,eq2,eq3,eq4], [0,0], 100000)
print time.clock() - start
pl.figure(figsize=(6,6))
pl.subplot(121)
pl.scatter(x, y, s=1, c="g", marker="s", linewidths=0)

pl.axis("equal")
pl.axis("off")
pl.subplot(122)
pl.scatter(x, y, s=1,c = c, marker="s", linewidths=0)
pl.axis("equal")
pl.axis("off")
pl.subplots_adjust(left=0,right=1,bottom=0,top=1,
wspace=0,hspace=0)
pl.gcf().patch.set_facecolor("white")
pl.show()

 


您可以捐助,支持我们的公益事业。

1元 10元 50元





认证码: 验证码,看不清楚?请点击刷新验证码 必填



1897 次浏览
14次
 捐助