其次用牛顿切线法求根，取不同的初值x0（例如-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集的程序。对不同的参数?p,q?：（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、给定仿射变换
?1?Z??sZ?1 ?2?Z??sZ?1 以及相应的概率p1?p2?1，其中s,Z为复数取s?0.5?0.5i，绘制相应的IFS迭代2的吸引子的图形。取不同的s，观察图形的变化。
4、请你概括一下分形的一般特征。
5、二维迭代分形。考虑函数f(x,y)?y?gnxbx?c与g(x,y)?a?x构成的二维迭代分形称为Martin迭代。现观察其当a?45,b?2,c?300时，取初值所得到的二维迭代散点图。程序如下：
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）取参数为其它的值时会得到什么图形？
第三章 数学建模实验
数学建模实验一 行星运动的加速度
一、实验目的
1、了解开普勒定律和牛顿定律，研究行星运动的规律。 2、理解数值微分的运用。
二、实验材料
2.1问题
我们希望知道行星在任一位置所受力的方向和大小，并研究这个力与行星所在位置的关系。根据牛顿第二定律，只须根据开普勒定律研究行星运动时的加速度的变化情况，也就知道了行星运动时的受力情况。
图3.1.1 星运动的轨道
2.2预备知识
开普勒根据他同时代的天文学家第谷.布拉赫对行星运动长期大量观测的数据，以及他自己对行星运动的观测和巨大的复杂计算，得出了行星运动的三大定律（图3.1.1）：
第一定律 行星运动的轨道是椭圆，太阳位于这个椭圆的一个焦点。 第二定律 从太阳到行星的所引的矢径在相同时间内扫过相同的面积。
第三定律 行星绕太阳的周期的平方与其轨道长半轴的立方成正比。
牛顿在前人（开普勒、胡克、雷恩、哈雷）研究的基础上，凭借他超凡的数学能力证明，1687年在专著《自然哲学的数学原理》上发表了万有引力定律。
牛顿万有引力定律 任意两个质点通过连心线方向上的力相互吸引。该引力的大小与它们的质量乘积成正比，与它们距离的平方成反比，与两物体的化学本质或物理状态以及中介物质无关。
牛顿由开普勒定律导出了万有引力定律，反过来，由牛顿万有引力定律和牛顿第二定律可以推导出在太阳引力下的天体的运动定律。开普勒定律和牛顿定律之间
的相互推导可以由解微分方程来实现。
2.3建立模型
将行星看作一个质点，以太阳的位置O为原点，在行星运动的椭圆平面内建立直角坐标系，写出椭圆方程。行星在任一时刻t的位置P用坐标(x,y)表示。P的位置随t的变化而变化。它的两个坐标x?x(t),y?y(t)都是t的函数。行星的运动可以分解为x轴方向上的分运动与y轴方向上的分运动的合运动。这两个分运动分别用函数x?x(t),y?y(t)来描述。
行星由点P1(x(t1),y(t1))运动到P2(x(t2),y(t2))的时间?t?t2?t1与矢径
??OP设行星绕椭圆运行一周的1到OP2所扫过的椭圆内部分的面积S(P1,P2)成正比。
Tk?周期为T,椭圆面积为S,则?t?kS(P，其中。因此，P的位置(x,y)随,P)12St变化的函数关系是已知的，描述两个分运动的函数x?x(t),y?y(t)也是已知的。只要求出这两个分运动在任一时刻t的加速度ax(t),ay(t)，也就是求出了行星运动的加速度矢量a(t)在x轴方向和y轴方向上的分量，a(t)也就知道了。研究加速度矢量a(t)的方向和大小与行星的关系，就可以得出万有引力定律。
怎样求在任一时刻t?t0加速度的两个分量ax和ay？以求ax为例：在t0前后很短的一段时间里，可以近似地将x?x(t)所描述的运动看作匀加速运动，求出它的加速度作为t?t0时的加速度。按照高中物理知识，具体求法可为：设
t?1?t0?t1，t1?t?1很小，并且已经分别知道t?t?1,t0,t1时的x的值x?1,x0,x1，则从t?1到t0这段时间内的平均速度为
x?x?1 v?0.5?0
t0?t?1t?t0它就是在时刻t?0.5??1的瞬时速度。同样，从t0到t1的平均速度为
2x?x0 v0.5?1
t1?t0t?t1它就是在时刻t0.5?0的瞬时速度。于是加速度在x轴方向上的分量为
2v?v?0.5 ax(t0)?0.5
t0.5?t?0.5同样可求出加速度在y轴方向上的分量ay(t0)。
上面所求的加速度分量ax(t0)，ay(t0)就是函数x(t)和y(t)在t?t0对t的二阶导数。上面所说的方法实际上就是：在t?t0的附近的一个很短的区间内，将二阶导数认为是常数。将函数近似地看作不超过二次的多项式函数
x(t)?k0?k1t?k2t2。已知xi?x(ti)(i??1,0,1)就可以求出系数k0,k1,k2，从
而求出其二阶导数a?2k2作为t?t0时二阶导数的值。
本实验求函数x?x(t)在给定的点的二阶导数的方法称为数值微分法。一般地说，要求函数y?f(x)在给定点x?x0的导数，可以按导数的定义，取x的一个很小的增量?，计算差商
f(x0??)?f(x0)?
作为导数的近似值。这实际上是将函数在区间[x0,x0??]内近似地当做一次函数来求它的导数。如果希望更准确一些，取x?1?x0?x1使x1?x0和x0?x?1都很小而且比较接近，将函数在区间[x?1,x1]上近似地当做不超过二次的多项式函数
y?k2x2?k1x?k0，由函数在x?1,x0,x1三点的值求出待定的系数k0,k1,k2，然后就可求出x0点的一阶导数2k2x0?k1和二阶导数2k2。如果希望进一步提高精确度，可以在更多的点用更高次的多项式函数去拟合y?f(x)。 2.4练习题
设行星轨道椭圆半长轴和半短轴分别为a,b，半焦距为c?a2?b2，行星运
动周期为T，设太阳位置为定点O，行星位置为动点P。根据开普勒定律研究行星运动时的加速度的变化情况。
三、实验准备
认真阅读实验目的与实验材料后要正确地解读实验，在此基础上制定实验计划（修改、补充或编写程序，提出实验思路，明确实验步骤），为上机实验做好准备。
四、实验思路提示
4.1实验步骤
设行星轨道椭圆半长轴和半短轴分别为a,b，半焦距为c?a2?b2，行星运
动周期为T，设太阳位置为定点O，行星位置为动点P。
（1）以太阳位置O为原点建立直角坐标系，使椭圆的中心坐标为(?c,0)。写出椭圆的参数方程
?x?acos??c 0???2? ??y?bsin?Mathematica程序如下：
Clear[x,t]
c=Sqrt[a^2-b^2] x[t_]=a*Cos[t]-c y[t_]=b*Sin[t]
（2）假定行星从??0的点(a?c,0)开始沿逆时针方向（即?增加的方向）运动。按行星运动的先后顺序椭圆上依次取N个点（比如N?100）
Pi(xi,yi)(i?0,1,2,?,N?1),对应的参数0??0??1????N?1?2?，使所有的?i?1??i都很小，每两个相邻的点Pi、Pi?1非常接近（约定?N?2?,PN?P0）。
Mathematica程序如下：
td=Table[2*i*Pi/100,{i,0,100}];
（3）用三角形PiOPi到OPi?1所扫过的面积： i?1的面积?Si近似地代表矢径OP??1xiyi1?(xiyi?1?xi?1yi)
2xi?1yi?12 根据开普勒第二定律，行星从Pi到Pi?1所花的时间?ti与Si成正比。行星运行一圈的时间等于周期T，矢径所扫过的面积就是整个椭圆的面积S??ab，取常数??T/?ab，则?ti???Si。
?Si?Mathematica程序如下：
Clear[x1,x2,y1,y2,ds]
ds[x1_,y1_,x2_,y2_]=1/2*(x1*y2-x2*y1); S=Pi*a*b;
T=Sqrt[k*a^3]; u=T/S;
dt[x1_,y1_,x2_,y2_]=u*ds[x1,y1,x2,y2]
下面计算P1的加速度，先计算?t0,?t1，分别放在dt1和dt2中
dt1=dt[x[td[[1]]],y[td[[1]]],x[td[[2]]],y[td[[2]]]]
dt2=dt[x[td[[2]]],y[td[[2]]],x[td[[3]]],y[td[[3]]]]
（4）将行星从Pi?1到Pi?1的运动看作匀加速运动，根据Pi?1、Pi、Pi?1三点的坐标及?ti?1,?ti，计算出加速度在x轴方向上的分量和y轴方向上的分量，作为行星在Pi点的加速度分量，从而确定了行星在点Pi的加速度矢量。
例如，计算行星在Pi点的加速度在x轴方向上的分量如下：
已经分别知道t?ti?1,ti,ti?1时的x的值xi?1,xi,xi?1，则从ti?1到ti这段时间内的平均速度为
vi?0.5?它就是在时刻ti?0.5?xi?xi?1
ti?ti?1ti?1?ti的瞬时速度。同样，从ti到ti?1的平均速度为 2x?xi vi?0.5?i?1
ti?1?tit?ti?1它就是在时刻ti?0.5?i的瞬时速度。
2如果Pi?1、Pi、Pi?1三点分别是P0、P1、P2，Mathematica程序如下：
Clear[v,x1,x2,t1,t2,v1,v2]
v[x1_,x2_,t1_,t2_]=(x2-x1)/(t2-t1); t1=0;
t2=t1+dt1; t3=t2+dt2;