数学实验报告 数学建模

数学实验报告

班级：应用物理091

学号：0902052016

姓名：李海瑞

2011.11.26

探索实验二非线性迭代

一，实验指导书解读

迭代是数学研究中的一个非常重要的工具，通过函数或向量函数由初始结点生成迭代结

点列，也可通过函数或向量函数由初值（向量）生成迭代数列或向量列。本实验主要分四步：首先利用蛛网图和迭代数列研究不动点的类型；

其次通过蛛网图和迭代数列研究Logistic映射，探索周期点的性质、认识混沌现象；

第三通过迭代数列或向量列求解方程（组）而寻求有效的求解方法；

最后，利用结点迭代探索分形的性质。

二,试验计划

2.1.迭代序列与不动点

对函数的迭代过程，我们可以用几何图象来直观地显示它——“蜘蛛网”。运行下列

Mathematica程序：

Clear[f]

f[x_] := (25*x - 85)/(x + 3); （实验时需改变函数）

Solve[f[x]==x , x] （求出函数的不动点）

g1=Plot[f[x], {x, -10, 20}, PlotStyle -> RGBColor[1, 0, 0],

DisplayFunction -> Identity];

g2=Plot[x, {x, -10, 10}, PlotStyle -> RGBColor[0, 1, 0],

DisplayFunction -> Identity];

x0=5.5; r = {};

r0=Graphics[{RGBColor[0, 0, 1],

Line[{{x0, 0}, {x0, x0}}]}];

For[i = 1, i <= 100, i++,

r=Append[r, 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}, （PlotRange控制图形上下范围）

DisplayFunction -> $DisplayFunction]

x[0]=x0;

x[i_]:=f[x[i-1]]; （定义序列）

t=Table[x[i],{i,1,10}]//N

ListPlot[t] （散点图）

如果只需迭代n次产生相应的序列，用下列Mathematica程序：

Iterate[f_,x0_,n_Integer]:=

Module[{ t={},temp= x0},AppendTo[t,temp];

For[i=1,i <= n, i++,temp= f[temp];

AppendTo[t,temp]];

t

]

f[x_]:= (x+ 2/x)/2;

Iterate[f,0.7,10]

2.2.Logistic映射与混沌

Mathematica程序：

IterGeo[a_, x0_] :=

Module[

{p1, p2, i, pointlist = {}, v= x0, fv= a*x0*(1 - x0)},

p1=Plot[ {a*x*(1 - x), x}, {x, 0, 1}, DisplayFunction -> Identity];

AppendTo[pointlist, {x0, 0}];

For[i = 1, i < 20, i++, AppendTo[pointlist, {v, fv}];

AppendTo[pointlist, {fv, fv}];

v= fv; fv= 4*v*(1 - v)];

p2=ListPlot[pointlist, PlotJoined -> True,

DisplayFunction -> Identity];

Show[{p1, p2}, DisplayFunction -> $DisplayFunction]

]

IterGeo[2.6, 0.3]

将区间（0，4］以某个步长a

?离散化，对每个离散的a值做迭代（2.2.2），忽略前50个迭代值，而把点()51,x

a，()52

,x

a显示在坐标平面上，最后形成的图形称为

a，…，()100

,x

Feigenbaum图。Mathematica程序：

Clear[f, a, x]; f[a_, x_] := a*x*(1 - x);

x0 = 0.5; r = {};

Do[

For[i = 1, i <= 300, i++,

x0 = f[a, x0];

If[i > 100, r = Append[r, {a, x0}]]

],

{a, 3.0, 4.0, 0.01}];

ListPlot[r]

在Logistic映射中，取4

a，任取两个初值使得它们之间的差的绝对值不超过0.l，运

=

行下列程序。

Sensitivity[n_Integer, x01_, x02_] :=

Module[

{pilist = {}, i, temp1=x01, temp2=x02},

For[i=1, i <= n, i++, temp1=4*temp1*(1-temp1);

temp2=4*temp2*(1-temp2);

AppendTo[pilist, {i, temp2-temp1}];

];

ListPlot[pilist, PlotJoined -> True]

]

Sensitivity[50, 0.1, 0.1001]

混沌不等于随机。实际上，在混沌区域之内，蕴涵着许多有序的规律。其Mathematica 程序：

distrib[n_Integer,m_Integer,x0_]:=

Module[{i,temp=x0,g1,f,k,c=Table[0,{i,m}]},For[i=1,i n,i++,temp=4*temp*(1-temp);

If[temp 1,c[[m]]++,c[[Floor[temp*m]+1]]++]];

f[k_]:=Graphics[{GrayLevel[0.5],Rectangle[{k-0.5,0},{k+0.5,c[[k]]}]}];

g1=Table[f[k],{k,1,m}];Show[g1,Axes True,PlotLabel->"x0=0.4"];]

n=100;m=20;x0=0.4;distrib[n,m,x0]

运行下列程序，听一听混沌的声音

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]

]

2.3. 方程求根

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]

线性方程组f x M x +=由f x M x n n

+=+1（ ,2,1,0=n ）给出迭代的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.分形

计算机绘出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]]

计算机绘出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]]

计算机绘出Sierpinski三角形的Mathematica程序：

redosierpinski[ptlist_List] :=

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]]

计算机绘出树木花草的Mathematica如下：

redotree[ptlist_List] :=

Block[{tmp = {}, i, ptnum = Length[ptlist]/2,

midpt1, midpt2, leftpt, rightpt},

For[i = 0, i < ptnum, i = i + 1,

midpt1 = (ptlist[[2i + 1]]*2 + ptlist[[2i + 2]])/3;

midpt2 = (ptlist[[2i + 1]] + ptlist[[2i + 2]]*2)/3;

leftpt = midpt1 + {{Cos[theta], -Sin[theta]},

{Sin[theta], Cos[theta]}}.{ptlist[[2i + 2]][[1]] - ptlist[[2i + 1]][[1]],

ptlist[[2i + 2]][[2]] - ptlist[[2i + 1]][[2]]}/3;

rightpt =

midpt2 + {{Cos[theta], Sin[theta]}, {-Sin[theta],

Cos[theta]}}.{ptlist[[2i + 2]][[1]] - ptlist[[2i + 1]][[1]],

ptlist[[2i + 2]][[2]] - ptlist[[2i + 1]][[2]]}/3;

tmp = Join[tmp, {ptlist[[2i + 1]], midpt1, midpt1, leftpt,

midpt1, midpt2, midpt2, rightpt,midpt2, ptlist[[2i + 2]]}]];

tmp]

showtree[ptlist_List] :=

Block[{tmp = {}, i, ptnum = Length[ptlist]/2},

For[i = 0, i < ptnum, i = i + 1,

AppendTo[tmp,

Line[{ptlist[[2i + 1]], ptlist[[2i + 2]]}]]];

Show[Graphics[tmp], AspectRatio -> 3/2/Sin[theta]]

]

theta = 30Degree;

showtree[Nest[redotree, {{0, 0}, {0, 1}}, 4]]

种颜色显示相应象素（或者用相应的灰度级显示）。

用第K

n mod

Mandelbrot集以及它的局部放大的Mathematica程序如下：

iter[x_,y_,lim_] := Block[{c,z,ct},c = x +I*y;z = c;ct = 0;

While[(Abs[z] < 2.0) && (ct < lim),++ct;z = z*z + c;];Return[ct];]

Mandelbrot1=DensityPlot[iter[x,y,50], {x,-2.0, 1.0}，{y,-1.5,1.5},

PlotPoints -> 120,Mesh -> False]

Mandelbrot2 = Show[Mandelbrot1, Graphics[Line[{{-0.9, -0.25},

{-0.7,-0.25}, {-0.7,-0.05},{-0.9,-0.05}, {-0.9，-0.25}}]]]

Mandelbrot3 = DensityPlot[iter[x,y,50],{x,-0.9,-0.7},{y,-0.25,-0.05},

PlotPoints -> 120,Mesh -> False]

Julia 集图形的Mathematica 程序：

julia[x_,y_,lim_,cx_,cy_] :=Block[{z, ct = 0}, z = x + I*y;

While[(Abs[z] <2.0) && (ct < lim), ++ct; z = z*z + (cx +I*cy);];Return[ct];] julia1=DensityPlot[julia[x, y,50, 0.27334, 0.00742],{x,-1.5,1.5}, {y,-1.5,1.5}, PlotPoints -> 120, Mesh -> False]

julia2 = Show[julia1,Graphics[Line[{{-0.7, -0.1}, {-0.3,-0.1},{-0.3,0.3},{-0.7,0.3},{-0.7,-0.1}}]]]

julia3 = DensityPlot[julia[x,y,50,0.27334,0.00742],{x,-0.7,-0.4},{y,-0.1,0.3}, PlotPoints -> 120,Mesh -> False]

给出IFS 迭代产生的分形的L 级灰度图象。

Mathematica 程序如下：

pl=0.3;aaa=1/2 +I1//N;f1[z_]:=(z+aaa)/2; p2=0.3;bbb=0//N;f2[z_]:=(z+bbb)/2; p3=0.4;ccc=1//N; f3[z_]:=(z+ccc)/2; f[z_]:=Block[{tmp},tmp=Random[];

Which[tmp

ShowIFS[z0_,MeshRange_List,divi_List,nmax_]:=

Block[{i,j,z=z0,a=divi[[1]],b=divi[[2]],temp1,temp2,mumax=0}, For[i=a,i>=1,i--,For[j=b,j>=1,j--,mu[i, j]= 0]]; For[i=nmax,i>=1,i--,

temp1=Floor[divi[[1]]*(Re[z]- MeshRange[[1]] [[1]])/( MeshRange[[2]] [[1]]- MeshRange[[1]][[1]])]+1; temp2=

Floor[divi[[2]]*(Im[z]- MeshRange[[1]] [[2]])/( MeshRange[[2]] [[2]]- MeshRange[[1]] [[2]])]+1; mu[temp1,temp2]++; z=f[z]; ];

For[i=a,i>=1,i--,

For[j=b,j>= 1,j--,mumax=Max[mumax, mu[i,j]]]]; mul=Table[GrayLevel[1-N[mu[j,i]]/mumax],{i,a},{j,b}]; Show[Graphics[RasterArray[mu1]]] ]

ShowIFS[0+I0,{{-0.1,-0.1},{1.1,1.1}},{150,150},100000]

三， 实验思路

3.1.迭代序列与不动点

首先对函数3

85

25+-=

x x y 研究不动点，需要 （1）Plot 中{x,-10,20}改为{x,-50,50}；对PlotRange 中{ -1,20}改为{-50,50} （2）（2）x0=5.5中5.5分别改为-30，-50： （3）对t=Table[x[i],{i,1,10}]//N 中10分别改为100； （4）对i<=100中100分别改为200，500。

运行程序后观察蛛网图与散点图！一看数列是否收敛？如收敛，极限是多少？收敛速度

是快是慢？二看蛛网图中的轨道是否趋于平衡点？与平衡点处曲线的斜率有没有关系？三看初值对结果有没有影响？

3.2. Logistic 映射与混沌

就Logistic 映射，对a =0.5，1，1.2，2，2.1，2.9，2.999，3，3.001，3.2，3.235，3.236，3.237，3.44等，分别取x 0= 0，0.2，0.5，0.8，1.0运行程序，观察结果。

观察结果就是看数列是否收敛，蛛网图中的轨道是否趋于平衡点，与a 的关系！对a 的定义范围[0，4]分成若干个区间，就初值（属于（0，1）时）看数列是否收敛，蛛网图中的轨道是否趋于平衡点？可用散点图认识。 对Logistic 映射讨论下列问题：

1）找出一个a 值，它对应的迭代具有2周期点。这种性质依赖于初值吗？你能找到多个a 值具有这种性质吗？

2）你能对任意的k 找到一个a 值，使得它对应的迭代具有k 周期点吗？哪些k 值能给出k 周期点？在每种情况下，结果是否依赖于初值的选取？

3）如果某个a 值能给出周期点，它是否一定是吸引的周期点？你能否找到排斥的周期点？

4）试着从理论上分析：()x f 的不动点是什么？对哪些a 值迭代收敛到每个不动点？哪些初值收敛到不动点？哪些初值导致发散？对周期点做类似的分析。

3.3. 方程求根

对于方程0123

=+-x x ，首先分别考虑函数21

3+x ，212x

x -，312-x ，取不同的初值x0

（例如-9,-5,-2,-1,-0.5,0,0.5,1, 2 ,8,10 等），生成相应的数列。若数列有极限，则极限即为所求根。

其次用牛顿切线法求根，取不同的初值x 0（例如-9,-5,-2,-1,-0.5,0,0.5,1, 2,8,10等），由此数列的极限即为所求根。

思考：为何数列收敛？分析原因！此外，初值x0决定数列的极限吗？ 最后考虑用其它方法求根。

而对于sinx -x +1=0及其它方程求根，类同上述步骤。

3.4. 分形

1、用Koch 曲线的生成元做迭代得到的极限图形称为Koch 雪花曲线。 （1）试用计算机画出Koch 雪花曲线；

（2）试计算雪花曲线的边长及面积，它们是否有限？你如何解释所得出的结论？ （3）雪花曲线是否光滑(即每一点是否有切线存在)? （4）其它的一些分形是否具有类似的性质？

2、编写绘制Julia 集的程序。对不同的参数()q p ,：（0，1），（-1，0），（0.11，0.66），（-0.102 81，0.957 23），（-1.25，-0。01）观察Julia 集的变化。取Julia 集的不同局部放大，你能看到某种自相似现象吗？

绘制Mandelbrot 集。然后，任意选取它的一个局部将其放大，然后再将放大图形的不同局部放大。由此观察Mandelbrot 集与Julia 集有何关系。进一步，取参数μ位于Mandelbrot 集的不同部位（如内部、外部、边界、某个苞芽的内部等），现察相应的Julia 集的形状的变化。

3、给定仿射变换

()11+=Z s Z ω ()12-=Z s Z ω

以及相应的概率2

1

21=

=p p ，其中Z s ,为复数取i s 5.05.0+=，绘制相应的IFS 迭代的吸引子的图形。取不同的s ，观察图形的变化。

4、请你概括一下分形的一般特征。

5、二维迭代分形。考虑函数c bx gnx y y x f --=),(与x a y x g -=),(构成的二维迭代分形称为Martin 迭代。现观察其当300,2,45===c b a 时，取初值所得到的二维迭代散点图。程序如下：

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]

问题：

（1）试着提高迭代次数至20000，50000，100000等观察图形有什么变化。 （2）取参数为其它的值时会得到什么图形？

四，实验过程与结果

4.1.迭代序列与不动点

Mathematica 程序：

Clear[f]

f[x_] := (25*x - 85)/(x + 3); （实验时需改变函数） Solve[f[x]==x , x] （求出函数的不动点）

g1=Plot[f[x], {x, -10, 20}, PlotStyle -> RGBColor[1, 0, 0], DisplayFunction -> Identity];

g2=Plot[x, {x, -10, 10}, PlotStyle -> RGBColor[0, 1, 0], DisplayFunction -> Identity]; x0=5.5; r = {};

r0=Graphics[{RGBColor[0, 0, 1], Line[{{x0, 0}, {x0, x0}}]}]; For[i = 1, i <= 100, i++,

r=Append[r, 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}, （PlotRange 控制图形上下范围） DisplayFunction -> $DisplayFunction] x[0]=x0;

x[i_]:=f[x[i-1]]; （定义序列） t=Table[x[i],{i,1,10}]//N ListPlot[t] （散点图） 运行结果：

{{x 5},{x 17}}

对Plot 中{x,-10,20}改为{x,-50,50}；对PlotRange 中{ -1,20}改为{-50,50} 运行结果：

10

55101520

5

10

15

20

51015

8

10

12

1416402020

40

40

20

20

40

5101520

8

10

12

14

16

5

10

15

20

8

10

12

14

16

20

15

10

5

402020

16

14

12

10

8

5101520 x0=5.5中5.5分别改为-30，-50，运行结果：

20

15

10

5

302010102030

对t=Table[x[i],{i,1,10}]//N 中10分别改为100；运行结果：

对i<=100中100分别改为200，500，运行结果：

5101520

8

10

12

14

16

1055101520

5

10

15

20

1055101520

5

10

15

20

20

406080100

17

17

17

17

实验观察：

1.一看数列是否收敛？如收敛，极限是多少？收敛速度是快是慢？

16

14

12

10

8

5101520 20

15

10

5

5101520

2.二看蛛网图中的轨道是否趋于平衡点？与平衡点处曲线的斜率有没有关系？

3.三看初值对结果有没有影响？

2.2Logistic映射与混沌

Mathematica程序：

IterGeo[a_, x0_] :=

Module[

{p1, p2, i, pointlist = {}, v= x0, fv= a*x0*(1 - x0)},

p1=Plot[ {a*x*(1 - x), x}, {x, 0, 1}, DisplayFunction -> Identity];

AppendTo[pointlist, {x0, 0}];

For[i = 1, i < 20, i++, AppendTo[pointlist, {v, fv}];

AppendTo[pointlist, {fv, fv}];

v= fv; fv= 4*v*(1 - v)];

p2=ListPlot[pointlist, PlotJoined -> True,

DisplayFunction -> Identity];

Show[{p1, p2}, DisplayFunction -> $DisplayFunction]

]

IterGeo[2.6, 0.3]

运行结果：

Mathematica程序：

Clear[f, a, x]; f[a_, x_] := a*x*(1 - x);

x0 = 0.5; r = {};

Do[

For[i = 1, i <= 300, i++,

x0 = f[a, x0];

If[i > 100, r = Append[r, {a, x0}]]

],

{a, 3.0, 4.0, 0.01}];

ListPlot[r]

运行结果：

1

0.8

0.6

0.4

0.2

0.20.40.60.8

其Mathematica 程序：

Sensitivity[n_Integer, x01_, x02_] :=

Module[

{pilist = {}, i, temp1=x01, temp2=x02},

For[i=1, i <= n, i++, temp1=4*temp1*(1-temp1); temp2=4*temp2*(1-temp2);

AppendTo[pilist, {i, temp2-temp1}]; ];

ListPlot[pilist, PlotJoined -> True] ]

Sensitivity[50, 0.1, 0.1001] 运行结果：

1020304050

1.0

0.5

0.5

Mathematica程序：

distrib[n_Integer,m_Integer,x0_]:=

Module[{i,temp=x0,g1,f,k,c=Table[0,{i,m}]},For[i=1,i n,i++,temp=4*temp*(1-temp);

If[temp 1,c[[m]]++,c[[Floor[temp*m]+1]]++]];

f[k_]:=Graphics[{GrayLevel[0.5],Rectangle[{k-0.5,0},{k+0.5,c[[k]]}]}];

g1=Table[f[k],{k,1,m}];Show[g1,Axes True,PlotLabel->"x0=0.4"];]

n=100;m=20;x0=0.4;distrib[n,m,x0]

运行结果：

实验观察：1.

5

10

15

2.5

57.5

1012.5

15x0 0.4

51015

20

406080100x0 0.5

5

1015

2

46810x0 0.8

2.3 方程求根

可由如下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] 运行结果： 线性方程组f x M x +=由f x M x n n

+=+1（ ,2,1,0=n ）给出迭代的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]

实验观察：

1345

20

10

10

20

2.4 分形

计算机绘出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]] 运行结果：

计算机绘出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,