实验四函数的迭代、混沌与分形解读

实验四函数的迭代、混沌与分形
[实验目的]
1. 认识函数的迭代；
2. 了解混沌和分形.
迭代在数值计算中占有很重要的地位,了解和掌握它是很有必要的.本实验将讨论用Newton迭代求方程根的问题,以及迭代本身一些有趣的现象.
§1 基本理论
1.1 迭代的概念
给定某个初值,反复作用以同一个函数的过程称为迭代.函数f(x)的迭代过程如下:
x0,x1=f(x0),x2=f(x1),……..,x n=f(x n-1)…..,
它生成了一个序列{x n}迭代序列.
许多由递推关系给出的数列,都是递推序列.例如数列.
X0=1,x n=1+1/(1+x n-1) (n=1,2,…………..)
是由函数f(x)=1+1/(1+x)=(2+x)/(1+x)取初值为1所得的迭代序列.
1.2 迭代序列的收敛性
定理设函数f(x)满足:
(1)对任意x∈(a,b),f(x)∈(a,b);
(2)f(x)在(a,b)内可导,且存在常数L使得|f(x)'|=L<1,
则当初值x0∈(a,b)时,由f(x)生成的迭代序列收敛.
在迭代函数f(x)连续的条件下,如果迭代数列收敛,则它一定收敛于方程x=f(x)的根.该方程的根也称函数f(x)的不动点.
设x*为f(x)的不动点,f(x)'在x*的附近连续,若|f(x*)'|<1,则称不动点x*是稳定的;若f(x*)'=0,则称不动点x*是超稳定的.在超稳定点x*附近,迭代过程x n+1=f(x n)收敛到x*的速度是非常快的.
1.3 Newton迭代法
设函数g(x)具有一阶导数,且g(x)'≠0,则函数f(x)=x-g(x)/g(x)'的迭代称为Newton迭代,若函数f(x)存在不动点,则它一定是方程g(x)=0的根,故Newton迭代法可用来求方程g(x)=0的根.
§2 实验内容与练习
2.1 迭代的收敛
对于函数迭代,最重要的问题是迭代序列的收敛性.一般说,迭代序列是否收敛取决于迭代函数与初值.
作为一个例子，我们用来讨论用Newton迭代法求函数g（x）=（x-17）5/3（x-5）-2/3的
根，其Mathematica程序为：
Clear[g，x]；
g[x_]：=（x-17）^(5/3)*(x-5)^(-2/3);
f[x_]=Factor[x-g(x)]/D[g[x],x]];
x0=5.5;n=20;
For[i=1,i<=n,i++,
x0=N[f[x0]];
Print[i,”“,x0,”“,D[f[x],x]/.x->x0]
]
执行结果见表4.1.
表4.1的结果说明迭代序列收敛于g(x)的零点17.我们注意到程序中取的迭代处值为5.5,如果其它的数作为初值,所得的迭代序列是否收敛于17呢?我们可以取其它初值做实验,结果得到表4.2(表中第三列是迭代序列的前6位有效数字首次为17.0000的步数).
从表4.2中可看出,只要初值不取5,迭代序列都收敛于17,且收敛速度与初值的选取关系不大.
前面程序中使用的f(x)为g(x)的化简过的Newton迭代函数,用Mathematica命令可检查出它为(25x-85)/(x+3)(注意,这个式子扩充了原迭代函数在x=5,x=17处的定义),解方程f(x)=x.得到x=17,与x=5.即17和5是f(x)的两个不动点,有前面的讨论知这两个不动点是有区别的:对于17,不管初值取为多少(只要不为5),迭代序列总是收敛于它;而对于5,只要初值取为5时,迭代序列才以它为极限,这样一种现象在函数的迭代中普遍存在,为方便区分起见,我们给这样两种点各一个名称:像17这样的所有附近的点在迭代过程中都趋向于它的不动点,称为吸引点;而像5这样的所有附近的点在迭代过程中都远离它的不动点,称为排斥点.
上面的f(x)=(25x-85)/(x+3)是一个分式线性函数,对于一般的分式线性函数,迭代序列是否总是收敛呢?
练习1 编程判断函数f(x)=(x-1)/(x+1)的迭代序列是否收敛.
在上节我们已经指出,如果迭代序列收敛,一定收敛到函数的某个不动点,这就是说,迭代函数存在不动点是迭代序列收敛的必要条件.那么如果迭代函数存在不动点,迭代序列是否一定收敛呢?
练习2 先分别求出分式线性函数f1(x)=(x-1)/(x+3),f2(x)=(-x+15)/(x+1)的不动点,再编
程判断它们的迭代序列是否收敛.
运用上节的收敛定理可以证明:如果迭代函数在某不动点处具有连续的导数且导数值介于-1与1之间,那么取该不动点附近的点为初值所得到的迭代序列一定收敛到该不动点.
练习3 你能否说明为什么17是f(x)=(25x-85)/(x+3)的吸引点,而5是f(x)的排斥点?尽量多找些理由支持这个结论.
练习4 能否找到一个分式线性函数(ax+b)/(cx+d),使它产生的迭代序列收敛到给定的数?用这种办法计算2.
2.2迭代的”蜘蛛图”
对函数的迭代过程,我们可以用几何图象来直观地显示它.在xoy平面上,先作出函数y=f(x)与y=x的图象,对初值x0,在曲线y=f(x)上可确定一点p0,它以x0为横坐标,过p0引平行x轴的直线,设该直线与y=x交与点Ｑ1作平行于y轴的直线它与曲线y=f(x)的交点记为p1,重复上面的过程,就在曲线y=f(x)上得到点列p1,p2,……，如图４.１，不难知道，这些点的横坐标构成的序列x0,x1,x2,……，xn……就是迭代序列.若迭代序列收敛,则点列p1,p2,……趋向于y=f（x）与y=x的交点p*，因此迭代序列是否收敛，可以在图上观查出来，这种图因其形状像蜘蛛网而被称为“蜘蛛网”图。

图4.2显示了分式线性函数f（x）=（25x-85）/（x+3）取初值为5.5的迭代过程，从图中可以看出该迭代是收敛的，且收敛到不动点17。

图4.2的“蜘蛛网”图可通过下面的程序获得
Clear[f]；图4.2 函数f（x）=（25x-85）/（x+3）的
f[x_]:=(25*x-85)/(x+3);
g1=Plot[f[x],{x,-10,20}, PlotStyle->RGBColor[1,0,0],
DisplayFunction->Identity];
g2=Plot[x,{x,-10,20}, PlotStyle->RGBColor[0,1,0],
DisplayFunction->Identity];
x0=5.5;r={};
r0=Graphics[{RGBColor[0,0,1],Line[{{x0,x0},{x0,f[x0]},{f{x0},f[x0]}}]
}]
];
x0=f[x0]
];
Show[g1,g2,r,r0,PlotRange->{-1,20},
DisplayFunction->$DisplayFunction]
练习五通过观察图4.2或通过改变初值重画f（x）=（25x-85）/（x+3）的蜘蛛网图，你是否能说明为什么该函数迭代的收敛速度与初值的选取关系不大？对于其他收敛的分式线形函数的迭代，是否有类似的结论？
f[x_]:=(25*x-85)/(x+3);
g1=Plot[f[x],{x,-10,20}, PlotStyle->RGBColor[1,0,0],
DisplayFunction->Identity];
g2=Plot[x,{x,-10,20}, PlotStyle->RGBColor[0,1,0],
DisplayFunction->Identity];
x0=5.5;r={};
r0=Graphics[{RGBColor[0,0,1],Line[{{x0,x0},{x0,f[x0]},{f{x0},f[x0]}}]
}]
];
x0=f[x0]
];
Show[g1,g2,r,r0,PlotRange->{-1,20},
DisplayFunction->$DisplayFunction]
练习五通过观察图4.2或通过改变初值重画f（x）=（25x-85）/（x+3）的蜘蛛网图，你是否能说明为什么该函数迭代的收敛速度与初值的选取关系不大？对于其他收敛的分式线形函数的迭代，是否有类似的结论？
2.3认识混沌
迭代序列若不收敛，它可能出现两种情况：
1.迭代次数充分大时，迭代序列出现周期性重复。

即存在自然数N，k>0, 使x N+k=x N，
这样迭代序列便成为：
x0，x1，。

x N，x N+1，。

，x N+k-1，x N，x N+1，。

，x N+k-1。

.
此时，x N，x N+1，。

，x N+k-1称为周期为k的循环；而初始点x0称为预周期点。

例如，对函数f(x)=-2+sin1.5x取初值x=0的迭代，可画出其蜘蛛网图如4.3所示，由该图可判断该迭代得到了一个周期为2的循环，0是该循环的一个预周期点。

2.序列没有规律`杂乱无章，称之为混沌。

例如，图4.4是函数f(x)= -2+sin1.5x取初值为
-0.7的迭代图，可看出该迭代产生了混沌。

图4.3 迭代出现循环图4.4迭代出现混沌
混沌具有两个特性：非随机性和对初始值的敏感性。

若初始值产生微小的误差，则该误差随迭代序列次数呈指数性增长，因此尽管迭代序列由初值和迭代函数完全决定，但随迭代次数的增加，它与随机序列并无多大差别，故混沌又称做确定性的随机运动。

练习6 通过观察图形进一步了解函数f(x)=a+sinbx的迭代（多取几组参数及初值）。

练习7 下列函数的迭代是否会产生混沌？
（1）2x , 0<=x<=1/2,
f(x)=
2(1-x), 1/2<x<=1;
(2)
f(x)=1/2(x+a/x).
练习8 函数f(x)=ax(1-x)(0<=x<=1)称为Logistic映射，试从“蜘蛛网”图观察它取初值为x0=0.5产生的迭代序列的收敛性，将观察记录填入表4.3，若出现循环，请指出它的周期。

2.4 人口增长的Logistic模型
Logistic映射来源于对人口变化规律的研究，下面我们简单加以介绍。

世界人口的增长，受到自然资源的约束，设地球允许承载的总人数为，以x（t）表示时刻为t时的人口总人数。

我们容易想到，在t时刻到时刻t+△t的人口平均增长速度，与t时刻的人口数成正比，也与容许增长的人口数成正比，故
[x（t+△t）-x（t）]/ △t=kx（t）（x m-x(t)）, (4.1)
其中k为比例系数，表示人口增长系数相关的系数。

设以一年为一时间单位，当△t=1时，
上面的式子变为
x（t+1）-x（t）= kx（t）（x m-x(t)），
即
x（t+1）=（1+kx m）x（t）（1-kx（t）/1+kx m
k*x（t+1）/（1+kx m）=（1+kx m）k/1+kx m x（t）（1-kx（t）/1+kx m
以n代替t，并令x n= k/1+kx m x（n），则有
x n+1=（1+kx m）x n（1-x n）
令α=1+kx m，可得
x n+1=αx n（1-x n）
此式正是由Logistic映射构成的迭代式。

因此分析Logistic映射的迭代行为，对掌握人口增长的变化规律有重大意义。

如果我们能很好的控制参数α，也就控制了人口增长方式，
是人口数朝着有利的方向发展。

2.5 Feignbaum图
对于Logistic映射f（x）=αx（1-x），我们来做一个实验：首先取α的值为3，在（0，1）中随机取一数作为初值
x0进行迭代，共迭代300次左右，丢弃起始的100次迭代的数据，在图上绘出所有的点（α，x n）（n>100）。

然后慢慢地增加α值，每增加一次，都重复前面的步骤，一直增加
到α=4为止，这样得到的图形，称为Feignbaum图
图4.5是由α取步长0.01所绘制的Feignbaum图，其Mathmatica程序如下：
Clear[f，u，x]；f[u-，x-]：=u*x*（1-x）；
X0=0.5；r={}；
Do[
For[i=1，i〈=300，i++，
x0=f[u，x0]；
If[I]100，r=Append[r，{u，x0}]]
]，
{u，3.0，4.0，0.01}]；
ListPlot[r]
为了获得更细致的图形，可将α的步长再缩少。

不过，由于循环次数较多，运行程序时，需要耐心地等待，一般要等数十分钟。

Feignbaum图对于分析函数f（x）=αx（1-x）的迭代行为非常有用。

从图4.5中可以看出：较左部分是一些清晰的曲线段，这说明对该范围内的任一α值而言，当迭代进行到100次以后，迭代所得的x n只取有限的几个值，表明迭代序列构成了一个循环，其周期等于竖直的直线与图形交点个数；从左到右，每一段曲线到一定位置同时一分为二，表明迭代序列的循环周期再这些位置处增长了一倍，因而曲线的这种分叉，被称为倍周期分支，有时也叫Pitch-Fork分支；随着α的增长，出现分支位置的间隔越来越少，大约在α=3.57左右，分支数已看不清楚，这是便出现了混沌，由此向右的区域被称为混沌区域，但是并非对所有大于3.57的a,函数迭代都出现混沌,在图中可看到,混沌区域中有一些空白带,这些
空白带由若干段曲线构成,说明对于相应的a,迭代出现周期循环,因此这些空白
带成为混沌区域中的周期窗口.例如当a=3.84时,迭代序列出现了周期为3的循环,
因此对应于a在3.84附近的区域就是一个周期为3的循环窗口.
练习9 在混沌区域中,还存在着其他的循环窗口,你能否找出?并通过计算
验证你的结论.
2.6对Logistic映射的进一步讨论
为讨论方便,我们称使迭代函数的一阶导数为0的点为临界点.对于Logistic
映射,通过简单的计算知,不管a是什么数,其临界点的皆为x=0.5,记
f n(x)=f(f(f…..(f(x)))) (其中右端有n 个f)
则有
fn(x1)=xn+1
与
fn(x1)=f(xn)f(xn-1)…f(x1).
若x1,x2,…xn+1是一周期为n的循环,且序列中出现0.5,则由fn(x1)=0
知该循环是超稳定的,既(0,1)中的任一点x0都是它的预周期点,且以x0为出值
进行迭代,收敛到该循环的速度超过指数增长的速度.
练习10 取初值为0.5,当a=3.5时,将得到一个周期为4的循环,该循环序列中含有0.5,实验证之.又问:以(0,1)内任一数为初值,是否都能得到这个循环?
其收敛速度与迭代次数间有何关系?
将n(0.5)视为a 的函数,在以a为横坐标,x为纵坐标的图上画出f1(0.5),
f2(0.5),f3(0.5),f4(0.5),f5(0.5)的图形见图4.6,
练习11 分析图4.6,你会发现一些什么样的结论?并说明理由,
图4.6 n=1,2,3,4,5是,fn(0.5)的图形
2.7二维迭代与分形
我们称由两个二元函数f(x,y)与g(x,y)取初值(x0,y0)构成的迭代
{xn+1=f(xn,yn),
yn+1=g(xn,yn)
为一个二维迭代.
二维迭代产生的序列也存在收敛性问题,但一般来说要比一元函数产生的迭代
序列的收敛性复杂的多.我们通常借助于图形进行观察,这种图形由二维点列(xn,yn)构成,
亦即二维散点图,
例1 由函数f(x,y)=y-sinx与g(x,y)=3.1-x取初值为(一。

2,0)构成的迭代,可通过
下面的程序获得其散点图:
a=3.1;
xn=1.2;yn=0;g={{xn,yn}};
For[n=1,n<=100,n++,
XN=xn,yN=yn;
Xn=yN-Sin[xN];yn=a-xN;
G=Append[g,{xn,yn}];];
ListPlot[g,AspectRatio->Automatic]
运行后得到图4.7,
虽然从图4.7中不能断定该序列是否收敛,但我们却发现这个图形本身非常
有特点,它由四个形状相似的小图形构成,这些小图形之间经过旋转,折叠后
图4.7 二维迭代散点图(之一)
能够相互重合,数学上称这类图形为分形图,
“分形”一词是在1975年由美国IBM公司数学家Benoit B,Mandelbrot首先提出来的.
它含有”碎化,分裂”之意,1982年Mandelbrot给分形下了一个通俗的定义:组成部分以某方式
与整体相似的形体叫分形,分形的两个最基本的性质是:比例性和置换不变性,比。

比例性指的是在一定表度范围内显微放大任何部分其不规则程度相同;而置换不变性是指每一部分经移位,旋转,缩放后与其他任意部位
相似,自然界许多物体,如植物,云团,雪花等都具有分形的性质,
科学家们猜测:自然界这些复杂的结构有可能是像二维迭代那样由简单的规律产生的,
bx 与g(x,y)=a-x构成的二维迭代称为Martin迭代.
例2 由函数f(x,y)=y-sgnx c
现观察其当a=45,b=2,c=-300时,取初值为(0,0)所得到的二维迭代散点图,
将例1的程序修改为:
a=45;b=2;c=-300;
xn=0;yn=0;g={{0,0}};
For{n=1,n<=5000,n++,
XN=xn;yN=yn;
Xn=yN-Sign[xN]*N[Sqrt[Abs[b*xN-c]]];
Yn=a-xN;
G=Append[g,{xn,yn}];];
ListPlot[g,PlotStyle->RGBColor[1,0,0],
AspectRatio->Automatic]
运行该程序,得到图4.8
图4.6 二维迭代散点图(之二)
该图形是著名的Martin图形的初期状态,随着迭代次数的提高,图形将会发生奇妙的变化.
练习12 对例2,试着提高迭代次数至26000,28000,100000,500000等观察图形有什么变化.
练习13 取参数a,b,c为其他的值回得到什么图形?参考表4.4
表4.4 Martin。

迭代参数表
1.D[f[x],x]/.x->x0
F(x)在x0处的导数
2.g1=Plot[f[x],{x,-10,20},PlotStyle->RGBColor[1,0,0],
DisplayFunction->Identity]
将g1定义为一个图形,该图形是f(x)在{-10,20}上的一段曲线弧(由
{x,-10,20}说明),语句PlotStyle->RGBColor{1,0,0}说明图形为红色,而语句
DisplayFunction->Identity使它不显示于屏幕.
3.r0=Graphics[{RGBColor[0,0,1],Line[……]}]
表示r0为蓝色的图形（折线），
4．r =Append[r，Graphics[…]
表明r是一图形集（由Graphics[…]定义），如此定义的r必须要有初值，
本实验中r的初值为空图，即r={}。

5．Show[g，DisplayFunction—>$DisplayFunction]
显示已定义的函数目标的图形g。