数学实验报告班级:应用物理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}]//NListPlot[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,xa,()52,xa显示在坐标平面上,最后形成的图形称为a,…,()100,xFeigenbaum图。
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映射中,取4a,任取两个初值使得它们之间的差的绝对值不超过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 + 1AppendTo[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]]种颜色显示相应象素(或者用相应的灰度级显示)。