distrib[n_Integer,m_Integer,x0_]:=
-temp);
n=100;m=20;x0=0.4;distrib[n,m,x0]
->\
f[k_]:=Graphics[{GrayLevel[0.5],Rectangle[{k-0.5,0},{k+0.5,c[[k]]}]}];
另一个说明混沌不是随机的事实是,混沌区域有许多有序的窗口。将这些窗口放大可以看到令人振奋的自相似现象,同时还有许多周期轨道。在 Feigenbaum图的右部,你应当能看到一个由三条曲线穿过的空白带,它是一个“周期为 3的窗口”。你能找到其它窗口吗?它们的周期是什么?窗口里有什么图案?这些窗口跟上题的第二问中的k周期轨道有什么关系? 运行下列程序,听一听混沌的声音
PlayChaos[n_Integer, x0_] :=
Module[
{t = {}, i, temp = x0},
For[i = 1, i <= n, i++, temp = 4*temp*(1 - temp); AppendTo[t, Floor[temp*100]]];
ListPlay[t, PlayRange -> {0, 100}, SampleRate -> 5] ]
和函数f?x??ax?1?x?一样有着混沌行为的函数还很多。其中较简单的有“帐篷函数”和“锯齿函数”。“帐篷函数”T?x? 定义为
T?x???“锯齿函数”S?x?定义为 S?x???2x?2x?1?0?x?0.5
0.5?x?12x?2?2x?0?x?0.5
0.5?x?1容易验证,帐篷函数和锯齿函数有下列关系:
0?x?1 T?T?x???T?S?x?? 令h?x??sin2?x2,帐篷函数与f?x??4x?1?x?有下列关系:
fk?h?x???h?Tk?x??
2.3 方程求根
对于代数方程g(x)=0,其根可用下列程序求得
Solve[g(x)= = 0 , x] 也可用下列程序求得
g[x_]:=expr
Plot[g[x],{x,a,b}]
FindRoot[g(x)= = 0 , {x,x0}]
将方程g(x)?0改写为等价的方程x?f(x),然后选取一初值利用(2.2.1)作迭代,迭代数列{xn}收敛的极限就是方程g(x)?0的解。即求方程g(x)?0的根等价于求函
数f(x)的不动点。注意:由g(x)?0可得不同的等价的方程x?f(x)。例如取2x?1x3?1,2,也可取为32x?1。 g(x)?x?2x?1,而f(x)可取为2x
由于不动点分吸引型和排斥型,因此g(x)?0的根为f(x)的排斥不动点时,就
3不能通过迭代函数f(x)以及一个初值x0利用(2.2.1)迭代求解,因为此时得到的数列{xn}不收敛。这时需要新的方法,Newton切线法就是其中之一。其迭代数列Mathematica程序如下:
Iterate[f_,x0_,n_Integer]:= Module[{ t={},temp= x0}, AppendTo[t,temp];
For[i=1,i <= n, i++,temp=N[x0-f[x0]/h[x0]]; AppendTo[t,temp]]; t ]
f[x_]:=x^3-2*x+1; h[x_]=Dt[f[x],x]; Iterate[f,4,10]
而要通过几何直观观察,可由如下Mathematica程序实现:
Clear[f]
f[x_] := x^3-2*x+1;
g1 = Plot[f[x], {x,2, 5}, PlotStyle -> RGBColor[1, 0, 0], DisplayFunction -> Identity]; x0 = 4; r = {}; h[x_]=Dt[f[x],x];
For[i = 1, i <= 100, i++,
If[h[x0]≠0,x1=N[x0-f[x0]/h[x0],20]];
r = Append[r, Graphics[{RGBColor[0, 0, 1], Line[{{x0, 0},
{x0, f[x0]}, {x1, 0}}] }]]; x0 =x1 ];
Show[g1, r, PlotRange -> {-20, 20},
DisplayFunction -> $DisplayFunction]
用迭代法求解线性方程组的思想与方程求根的方法是类似的,用迭代方法也可以求更加复杂的非线性方程组的解。由于非线性方程组可能有许多解(甚至有无穷多个解),因此对它的求解比线性方程组的求解要面临更多的挑战。
?????? 线性方程组x?Mx?f由xn?1?Mxn?f(n?0,1,2,?)给出迭代的Mathematica程序:
LSIterate[m_,f_List,f0_List,n_Integer]:= Module[
{i,var=f0,t=Table[{},{i,n}]}, For[i=1,i<=n,i++,t[[i]]=var;var=m.var+f]; t ]
m={{0.2,0.3},{0.4,0.2}};f={1,1};f0={0,0}; LSIterate[m,f,f0,20]
2.4 分形
早在上世纪末及本世纪初,一些数学家就构造出一些边界形状极不光滑的图形。
由于这类图形长期以来被视为“不可名状的”或“病态的”,因而,只有当人们需要反例时才想到它们。这类图形的构造方式都有一个共同的特点,即最终图形F都是按照一定的规则R通过对初始图形F不断修改得到的。其中最具有代表性的图形是Koch曲线,Koch曲线的构造方式是:
给定一条直线段F0(图2.2.l上左),将该直线三等分,并将中间的一段用以该线段为边的等边三角形的另外两条边替代,得到图形F1,称F1为生成元(图2.2.l上右)。然后,再对图形F1中的每一小段都按上述方式修改(图2.2.l下左),以至无穷,则最后得到的极限曲线F?lim。 Fk即是所谓的Koch曲线(图2.2.l下右)k??
图2.2.l Koch曲线的生成
计算机绘出Koch曲线的Mathematica程序:
redokoch[ptlist_List] :=
Block[{tmp = {}, i, pnum = Length[ptlist]},
For[i = 1, i < pnum, i = i + 1, tmp = Join[tmp, {ptlist[[i]], ptlist[[i]]*2/3 + ptlist [[i + 1]]/3,
(ptlist[[i]] + ptlist [[i + 1]])/2 + {ptlist[[i]] [[2]] - ptlist[[i + 1]] [[2]], ptlist[[i + 1]][[1]] - ptlist[[i]][[1]]}*Sqrt[3]/6,
ptlist[[i]]/3 + ptlist[[i + 1]]*2/3,ptlist[[i + 1]]}]]; tmp] lnko01 = {{0, 0}, {1, 0 }};
Show[Graphics[Line[Nest[redokoch, lnko01, 5]], AspectRatio -> Sprt[3]/6]]
分形的基本特性完全由生成元决定,因此,给定一个生成元,就可以生成各种各样的分形图形。以下是几个经典的分形图形及其生成Mathematica程序:
图2.2.2 Minkowski“香肠”
计算机绘出Minkowski“香肠”的Mathematica程序:
redominkowski[ptlist_List] :=
Block[{tmp = {}, tmp1, i, pnum = Length[ptlist]}, For[i = 1, i < pnum, i = i + 1,
tmp1 = {ptlist[[i]][[2]] - ptlist[[i + 1]][[2]], ptlist[[i + 1]][[1]] - ptlist[[i]][[1]]}/4; tmp = Join[tmp, {ptlist[[i]],
ptlist[[i]]*3/4 + ptlist[[i + 1]]/4,
ptlist[[i]]*3/4 + ptlist[[i + 1]]/4 + tmp1, ptlist[[i]]/2 + ptlist[[i + 1]]/2 + tmp1, ptlist[[i]]/2 + ptlist[[i + 1]]/2,
ptlist[[i]]/2 + ptlist[[i + 1]]/2 - tmp1, ptlist[[i]]/4 + ptlist[[i + 1]]*3/4 - tmp1, ptlist[[i]]/4 + ptlist[[i + 1]]*3/4, ptlist[[i + 1]]}] ]; tmp ]
redomk1[ptlist_list] :=
Block[{tmp = {ptlist[[1]][[2]] - ptlist[[2]][[2]], ptlist[[2]][[1]] - ptlist[[1]][[1]]}/4} {ptlist[[1]],
ptlist[[1]]*3/4 + ptlist[[2]]/4,
ptlist[[1]]*3/4 + ptlist[[2]]/4 + tmp, ptlist[[1]]/2 + ptlist[[2]]/2 + tmp, ptlist[[1]]/2 + ptlist[[2]]/2,
ptlist[[1]]/2 + ptlist[[2]]/2 - tmp, ptlist[[1]]/4 + ptlist[[2]]*3/4 - mp, ptlist[[1]]/4 + ptlist[[2]]*3/4, ptlist[[2]] }]
redomk2[ptlist_list] :=
Block[{tmp = {}, i, pnum = Length[ptlist]}, For[i = 1, i < pnum, i = i + 1,
tmp = Join[tmp, redomk1[{ptlist[[i]], ptlist[[i + 1]]}]]]; tmp ]
In01 = {{0, 0}, {1, 0}};
Show[Graphics[Line[Nest[redominkowski, In01, 4]], AspectRatio -> 1/GoldenRatio]]
图2.2.3 Sierpinski三角形
计算机绘出Sierpinski三角形的Mathematica程序:
Block[{tmp = {}, i, pnum = Length[ptlist]/3}, For[i = 0, i < pnum, i = i + 1,
tmp = Join[tmp, {ptlist[[3i + 1]],
(ptlist[[3i + 1]] + ptlist[[3i + 2]])/ 2, (ptlist[[3i + 1]] + ptlist[[3i + 3]])/ 2, (ptlist[[3i + 1]] + ptlist[[3i + 2]])/2,
ptlist[[3i + 2]],
(ptlist[[3i + 2]] + ptlist[[3i + 3]])/ 2, (ptlist[[3i + 1]] + ptlist[[3i + 3]])/ 2, (ptlist[[3i + 2]] + ptlist[[3i + 3]])/2,
ptlist[[3i + 3]]}]]; tmp]
showsierpinski[ptlist_List] :=
Block[{tmp = {}, i, pnum = Length[ptlist]/3}, For[i = 0, i < pnum, i = i + 1
AppendTo[tmp,Polygon[{ptlist[[3*i + 1]], ptlist[[3*i + 2]], ptlist[[3*i + 3]]}]]];
Show[Graphics[tmp], AspectRatio -> 1/GoldenRatio]] po1 = {{-1, 0}, {1, 0}, {0, Sqrt[3]}};
showsierpinski[Nest[redosierpinski, po1, 4]]
redosierpinski[ptlist_List] :=
图2.2.4龙曲线
图2.2.5 Hilbert曲线