注册 登录  
 加关注
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

一路

To find the final symmetry & beauty

 
 
 

日志

 
 
 
 

分形与混沌计算小计  

2016-09-17 12:26:28|  分类: 计算之美 |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |
Logistic 方程
Logistic 方程是从研究动物种群繁衍得出的简化模型,表达示如下:分形与混沌计算小计 - saturnman - 一路
   不用关心模型本身的意义如何,Logistic方程本身的性质就极其丰富有趣。对于x处在区间(0,1)上,r处在(0,4)上时,方程反复迭代过程中不会发散。可以从几个例子研究方程性质。
1.当r在不同区间时,方程迭代最终收敛特性差别极大,会表现有趣的分形性质
%matplotlib inline
import sys
import os
import numpy as nm
import matplotlib.pyplot as plt
def calc_stable_x(start,r):
x = start
for i in range(0,997):
x = r*x*(1-x)
return x
def docalc():
r = nm.arange(2.001,3.8,0.0001)
s = nm.random.random_sample(r.shape)
x = calc_stable_x(s,r)
fig = plt.gcf()
fig.set_size_inches(18.5, 10.5)
plt.plot(r,x,'.',color="r",markersize=2.0,label="Logistic Eq")
plt.xlabel("r")
plt.ylabel("x")
plt.legend()
plt.show()
docalc()

分形与混沌计算小计 - saturnman - 一路
 
当r<=3时,迭代方程最终会收敛到单值。随着r的增大,方程终值最终变为双周期震荡, 当r大概处在>3.6的范围时,对于十分微小差别的初值,随着迭代的进行,x获得的值也会分离得非常远,因此系统通过观察后期迭代结果完全无法推断系统初始状态。
%matplotlib inline
import sys
import os
import numpy as nm
import matplotlib.pyplot as plt
def calculate_trajectory(seed,r,count):
result = []
x = seed
f = lambda r,x:r*x*(1-x)
for i in range(0,count):
x=f(r,x)
result.append(x)
return result
def main():
r = 4.0
length = 60
x = range(0,length)
y1 = calculate_trajectory(0.2,r,length)
y2 = calculate_trajectory(0.20000001,r,length)
fig = plt.gcf()
fig.set_size_inches(18.5, 10.5)
plot1 = plt.plot(x,y1,color="r",label="seed=0.2,r=4.0,steps=60")
plot2 = plt.plot(x,y2,color="g",label="seed=0.20000001,r=4.0,steps=60")
plt.xlabel("steps")
plt.ylabel("iter value")
plt.legend()
plt.show()
main()

分形与混沌计算小计 - saturnman - 一路
 
上图展示的是初值分别为2和2.0000001时的迭代过程,最开始迭代值几乎相等,曲线几乎处在重合状态,但是大概经过20次迭代之后,曲线开始分离,当达到60次之后,系统完全忘记了最初的处置,进入一种混沌的状态。注意到如果进行解析计算,系统总是有一个稳定点,就是1-1/r,如果当r>3时,还有稳定点分形与混沌计算小计 - saturnman - 一路分形与混沌计算小计 - saturnman - 一路而对于实际系统和数值系统是没有办法达到任意精度的,任何非常小的差别都会造成迭代方程最终分离到系统的稳定周期终值点或是达到混沌状态。
%matplotlib inline
import sys
import os
import numpy as nm
import matplotlib.pyplot as plt
def calculate_trajectory(seed,r,count):
result = []
x = seed
f = lambda r,x:r*x*(1-x)
for i in range(0,count):
x=f(r,x)
result.append(x)
return result
def main():
r = 4.0
length = 100
x = range(0,length)
y1 = calculate_trajectory(1-1/r+0.000000001,r,length)
fig = plt.gcf()
fig.set_size_inches(18.5, 10.5)
plot1 = plt.plot(x,y1,color="r",label="seed=1-1/r+0.000000001,r=4.0,steps=60")
plt.xlabel("steps")
plt.ylabel("iter value")
plt.legend()
plt.show()
main()

分形与混沌计算小计 - saturnman - 一路
 
上图为初值为与1-1/r差距极小的位置,前20次迭代差距几乎没有,但是20次之后差距突然拉大,最终达到混沌的状态。分形计算中存在一个重要常数,我们如果可以计算每次周期分叉时的r的值,那么分形与混沌计算小计 - saturnman - 一路
 将趋于一个定值,让人吃惊的是此参数在非常多不同背景下出现,是一个通用常数,称为费根鲍姆常数,其值约 4.66920。此值人们普遍认为应该是一个像e,π一样的超越数,不过无论是解析计算此参数或是数值计算此参数,都是非常困难的,目前没有发现高效的算法。
洛伦兹吸引子
洛伦兹在研究天气预报的一组简化了的微分方程组时发现,在数值计算过程中,如果初值设置有微小的变化,方程演化将出现非常不同的结果。而无论实际计算过程还是实际的物理系统,都会在微小尺度上存在非常多的不确定性。因此就结果来说,系统的演化长期预测无论从理论上还是实际中,都将变得不可能。洛伦兹的发现影响深远,这也是蝴蝶效应一词的由来。在洛伦兹具体计算演化方程时,发现某些参数情况下,系统的参量演化会围着几个相空间点周围附近环绕,时而跳到另一相空间点,时而跳跃回去,根本无法准确预测其行为。被围绕的相空间点称为洛伦兹吸引子。
#转自:http://blog.sciencenet.cn/blog-265205-535333.html
# -*- coding: utf-8 -*-
%matplotlib inline
from scipy.integrate import odeint
import numpy as np
def lorenz(w, t, p, r, b):
# 给出位置矢量w,和三个参数p, r, b计算出
# dx/dt, dy/dt, dz/dt的值
x, y, z = w
# 直接与lorenz的计算公式对应
return np.array([p*(y-x), x*(r-z)-y, x*y-b*z])
t = np.arange(0, 30, 0.01) # 创建时间点
# 调用ode对lorenz进行求解, 用两个不同的初始值
track1 = odeint(lorenz, (0.0, 1.00, 0.0), t, args=(10.0, 28.0, 3.0))
track2 = odeint(lorenz, (0.0, 1.01, 0.0), t, args=(10.0, 28.0, 3.0))
# 绘图
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
fig = plt.gcf()
ax = Axes3D(fig)
ax.plot(track1[:,0], track1[:,1], track1[:,2])
ax.plot(track2[:,0], track2[:,1], track2[:,2])
plt.show()

分形与混沌计算小计 - saturnman - 一路
 
关于分形图形计算

分形类的计算过程中有很多精美的图形,对计算机图形也有些影响,比如Mandelbrot集的计算,Mandelbrot集是复数迭代方程:分形与混沌计算小计 - saturnman - 一路从z=0开始迭代而让迭代结果不发散的复数c组成的集合。而不处在Mandelbrot集中的复数以它让迭代方程发散速度快慢着色,这样可以得到非常精美的图形。
1.Mandelbrot集图形计算
%matplotlib inline
import numpy
import matplotlib.pyplot

fig = matplotlib.pyplot.gcf()
fig.set_size_inches(150,75)

ITERATIONS = 500
#SIZE = lena.shape[0]
#注意:此参数会极大影响计算过程的内存消耗,对于小内存机器可以会让计算进程因内存不够而出错,请自行调节计算
SIZE = 3000
MAX_COLOR = 255.
x_min, x_max = -2.5, 1
y_min, y_max = -1.2, 1.2
x,y = numpy.meshgrid(numpy.linspace(x_min,x_max,2*SIZE),numpy.linspace(y_min,y_max,SIZE));
c = x + 1j*y
z = c.copy()
fractal = numpy.zeros(z.shape,dtype=numpy.uint8)+MAX_COLOR
for n in range(ITERATIONS):
mask = numpy.abs(z) <=10
z[mask] = z[mask]**2+c[mask]
fractal[(fractal==MAX_COLOR) & (~mask)] = (MAX_COLOR-1)*n/ITERATIONS
# Display the fractal
matplotlib.pyplot.imshow(fractal,cmap=matplotlib.pyplot.get_cmap('flag'))
matplotlib.pyplot.title('Mandelbrot')
matplotlib.pyplot.axis('off')

分形与混沌计算小计 - saturnman - 一路
 分形与混沌计算小计 - saturnman - 一路
 分形与混沌计算小计 - saturnman - 一路
 
2.分形曲线
#转自:http://blog.sciencenet.cn/blog-265205-535333.html
%matplotlib inline
# -*- coding: utf-8 -*-
#L-System(Lindenmayer system)是一种用字符串替代产生分形图形的算法
from math import sin, cos, pi
import matplotlib.pyplot as pl
from matplotlib import collections
class L_System(object):
def __init__(self, rule):
info = rule['S']
for i in range(rule['iter']):
ninfo = []
for c in info:
if c in rule:
ninfo.append(rule[c])
else:
ninfo.append(c)
info = "".join(ninfo)
self.rule = rule
self.info = info
def get_lines(self):
d = self.rule['direct']
a = self.rule['angle']
p = (0.0, 0.0)
l = 1.0
lines = []
stack = []
for c in self.info:
if c in "Ff":
r = d * pi / 180
t = p[0] + l*cos(r), p[1] + l*sin(r)
lines.append(((p[0], p[1]), (t[0], t[1])))
p = t
elif c == "+":
d += a
elif c == "-":
d -= a
elif c == "[":
stack.append((p,d))
elif c == "]":
p, d = stack[-1]
del stack[-1]
return lines
rules = [
{
"F":"F+F--F+F", "S":"F",
"direct":180,
"angle":60,
"iter":5,
"title":"Koch"
},
{
"X":"X+YF+", "Y":"-FX-Y", "S":"FX",
"direct":0,
"angle":90,
"iter":13,
"title":"Dragon"
},
{
"f":"F-f-F", "F":"f+F+f", "S":"f",
"direct":0,
"angle":60,
"iter":7,
"title":"Triangle"
},
{
"X":"F-[[X]+X]+F[+FX]-X", "F":"FF", "S":"X",
"direct":-45,
"angle":25,
"iter":6,
"title":"Plant"
},
{
"S":"X", "X":"-YF+XFX+FY-", "Y":"+XF-YFY-FX+",
"direct":0,
"angle":90,
"iter":6,
"title":"Hilbert"
},
{
"S":"L--F--L--F", "L":"+R-F-R+", "R":"-L+F+L-",
"direct":0,
"angle":45,
"iter":10,
"title":"Sierpinski"
},
]
def draw(ax, rule, iter=None):
if iter!=None:
rule["iter"] = iter
lines = L_System(rule).get_lines()
linecollections = collections.LineCollection(lines)
ax.add_collection(linecollections, autolim=True)
ax.axis("equal")
ax.set_axis_off()
ax.set_xlim(ax.dataLim.xmin, ax.dataLim.xmax)
ax.invert_yaxis()
fig = pl.gcf()
fig.patch.set_facecolor("w")
for i in xrange(6):
ax = fig.add_subplot(231+i)
draw(ax, rules[i])
fig.subplots_adjust(left=0,right=1,bottom=0,top=1,wspace=0,hspace=0)
pl.show()

分形与混沌计算小计 - saturnman - 一路
 
  评论这张
 
阅读(100)| 评论(0)
推荐 转载

历史上的今天

在LOFTER的更多文章

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2017