数学实验与数学建模实践教程

1、基础实验一 数列极限与函数极限一、实验目的从刘徽的割圆术、裴波那奇数列研究数列的收敛性并抽象出极限的定义；理解数列收敛的准则；理解函数极限与数列极限的关系。二、实验材料1.1割圆术中国古代数学家刘徽在九章算术注方田章圆田术中创造了割圆术计算圆周率。刘徽先注意到圆内接正多边形的面积小于圆面积；其次，当将边数屡次加倍时，正多边形的面积增大，边数愈大则正多边形面积愈近于圆的面积。“割之弥细，所失弥少。割之又割以至不可割，则与圆合体而无所失矣。”这几句话明确地表明了刘徽的极限思想。以表示单位圆的圆内接正多边形面积，则其极限为圆周率。用下列Mathematica程序可以从量和形两个角度考察数列的收敛情况

2、： m=2;n=15;k=10; Fori=2,i<=n,i+, li_:=N2*SinPi/(3*2i),k;(圆内接正多边形边长) si_:=N3*2(i-1)*li*Sqrt1-(li)2/4,k;(圆内接正多边形面积) ri_:=Pi-si; di_:=si-si-1; Printi," ",ri," ",li," ",si," ",di t=Tablei,si,i,m,n(数组)ListPlott(散点图)1.2裴波那奇数列和黄金分割 由有著名的裴波那奇数列。如果令，由递推公式可得出，；。用下列Ma

3、thematica程序可以从量和形两个角度考察数列的收敛情况：n=14,k=10;Fori=3,i<=n,i+,t1=(Sqrt5+1)/2;t2=(1-Sqrt5)/2;fi_:=N(t1(i+1)-t2(i+1)/Sqrt5,k; (定义裴波那奇数列通项) rn=(5(1/2)-1)/2-fi-1/fi;Rn=fi-1/fi;dn=fi-1/fi-fi-2/fi-1; Printi," ",rn," ",Rn," ",dn;t=Tablei,fi-1/fi,i,3,nListPlott1.3收敛与发散的数列数列当时收敛，时发

4、散；数列发散。1.4函数极限与数列极限的关系用Mathematica程序m=0;r=10m;x0=0;fx_=x*Sin1/xPlotfx,x,-r,rLimitfx,x->x0观察的图象可以发现，函数在点处不连续，且函数值不存在，但在点处有极限。令，作函数的取值表，画散点图看其子列的趋向情况k=10;p=25;an_=1/n;tf=Tablen,Nfan,k,n,1,pListPlottfLimitfan,nInfinity,Direction1分别取不同的数列(要求)，重做上述过程，并将各次所得图形的分析结果比较，可知各子列的极限值均为上述函数的极限值。对于，类似地考察在点处的极限。

5、三、实验准备认真阅读实验目的与实验材料后要正确地解读实验，在此基础上制定实验计划（修改、补充或编写程序，提出实验思路，明确实验步骤），为上机实验做好准备。四、实验思路提示改变或增大，观察更多的项（量、形），例如，分别取50，100，200，；扩展有效数字，观察随增大数列的变化趋势，例如，分别取20，30，50；或固定50；或随增大而适当增加。对实验要思考，例如，定义中的指标与柯西准则中的指标间的差异；数列收敛方式；又例如，如何估计极限近似值的误差。改变函数及极限类型，例如，考虑六种函数极限，既选取极限存在也选取极限不存在的例子；改变数列，改变参数观察更多的量，考察形的变化趋势；扩展有效数字，提

6、高计算精度。要对实验思考，归纳数列敛散与函数敛散的关系。基础实验二 定积分数值计算一、实验目的学习定积分的数值计算方法，理解定积分的定义，掌握牛顿-莱布尼兹公式。二、实验材料计算定积分的近似值，可将积分区间等分而得矩形公式或也可用梯形公式近似计算如果要准确些，可用辛普森公式 对于，矩形公式、梯形公式、辛普森公式的Mathematica程序为a=0;b=1;k=10;fx_:=Sinx;d=NIntegratefx,x,a,b,k;（计算精确值）s1m_:=NSumfa+i*(b-a)/m*(b-a)/m,i,0,m-1,k;（取小区间左端点的矩形公式）s2m_:=NSumfa+(i+1/2)*

7、(b-a)/m*(b-a)/m,i,0,m-1,k; （取小区间中点的矩形公式）s3m_:=NSumfa+i*(b-a)/m*(b-a)/m,i,1,m,k; （取小区间右端点的矩形公式）s4m_:=NSum(fa+i*(b-a)/m+fa+(i+1)*(b-a)/m)/2*(b-a)/m,i,0,m-1,k; （梯形公式）s5m_:=N(b-a)/m/6*(fa+fb)+2*Sumfa+i*(b-a)/m,i,1,m-1+4*Sumfa+(i-1/2)*(b-a)/m,i,1,m),k;（辛普森公式）r1m_:=d-s1m;r2m_:=d-s2m;r3m_:=d-s3m;r4m_:=d-s4

8、m;r5m_:=d-s5m;（误差）t=Tables1m,r1m,s2m,r2m,s3m,r3m,s4m,r4m,s5m,r5m,m,100,1000,100 利用以上程序计算、，并对几个公式比较。如果函数在区间上连续，则在区间上可积。反之不然。设函数在上连续，而且是的一个原函数，则有牛顿-莱布尼兹公式。函数在不连续、不存在原函数，但在上可积；函数在不连续，但在上可积、存在原函数。此外函数处处不连续、不存在原函数，在任意区间(长度大于0)上不可积。求原函数并验证牛顿-莱布尼兹公式的Mathematica程序fx_:=Sinx;Integratef(x),x（求不定积分） Fx_:=%（定义原函

9、数）d=NIntegratef(x),x,a,b（求定积分） df=Fb-Fa （计算原函数的增量） r=d-df三、实验准备 认真阅读实验目的与实验材料后要正确地解读实验，在此基础上制定实验计划（修改、补充或编写程序，提出实验思路，明确实验步骤），为上机实验做好准备。四、实验思路提示先对一个函数，例如在区间0,1，在程序中改变（例如、）并适当扩展有效数字（例如、），运行程序计算定积分的近似值，分析误差。再考虑其它函数。最后对几个公式比较。 先对一个函数，例如在区间0,1, 运行程序计算。再考虑其它函数，例如指数函数、分段连续函数、。分析可积条件及牛顿-莱布尼兹公式成立的条件。基础实验三 盈亏

10、转折与投入产出一、实验目的理解微观经济学的基本理论与方法，利用线性代数的有关理论和方法解决经济管理中的盈亏转折及投入产出问题。二、实验材料 已知某企业的产品数量与成本的若干数据如下：产品数量（百件）61020成本数量（千元）104160370设每件产品的出厂价=20（千元百件），判断企业盈亏转折时的产品数量的变化范围及企业获取最大利润额时的产品数量。设成本函数，其中为待定系数；产值函数 ，于是利润函数成本函数的导函数称为边际成本，利润函数的导函数称为边际利润。设与为方程的两个实根与称为盈亏转折点。当并且时，即使企业获得利润的产品数量的范围为；而在这个范围之外企业不能获得利润。利润函数的唯一极大

11、值点即的最大值点，就是使企业获得最大利润的产品数量。 Mathematica程序data=6，104，10，160，20，370；（原始数据） InterpolatingPolynomialdata，x（求内插多项式，即成本函数） Lx_：=2 0*x-（0.5*x2 6*x+50）（利润函数） RootsLx=0，x（利润函数的零点） FindMaximumLx，x，0（求利润函数的极大值） MATLAB程序>>x=6，10，20； %产品数量 y=104，160，370； %成本数量 c=polyfit（x，y，2） %拟合二次多项式>>c=0，20，0-c； rr

12、oots(c) %求盈亏转折点 >>r1=roots（polyder（c）%求微分多项式的根 >>L=polyval（c，rl）%求最大利润2.2投人产出分析 某地区有三个重要产业，一个煤矿，一个发电厂和一条地方铁路。；。在某一周内，煤矿接到外地金额为50000元的定货，发电厂接到外地金额为 25000元的定货，外界对地方铁路没有需求，问三个企业在这一周内总产值多少才能满足自身及外界的需求？ 设为煤矿本周内的总产值，为电厂本周内的总产值，为铁路本周内的总产值，则 （3.1）记矩阵称为直接消耗矩阵，称为产出向量，称为最后需求向量，则方程组（3.1）表示为或 （3.2）其中

13、矩阵为单位矩阵，称为列昂杰夫矩阵，列昂杰夫矩阵为非奇异矩阵。设，矩阵称为完全消耗矩阵，它与矩阵一起在各个部门之间的投入产出中起平衡作用。矩阵可以称为投入产出矩阵，它的元素表示煤矿，电厂，铁路之间的投入产出关系。向量称为总投入向量，它的元素是矩阵的对应列元素之和，分别表示煤矿，电厂，铁路得到的总投入。由矩阵，向量，和，可得投入产出分析表如下：表3.1 投入产出分析表（单位：元）煤矿电厂铁路外界需求总产出煤矿错误！链接无效。电厂铁路总投入Mathematica程序A2=1,-,-,-,1（列昂杰夫矩阵）A2.x1，x2，x3=50000，25000，0（方程组）Solve,xl,x2,x3（解方程

14、组，求出总产出）X2=,0,0,0,0,0,0,C2=（IdentityMatrix3-A2）.X2（求投入产出矩阵）D2=1,1,l。C2（求总投入） MATLAB程序>>A= 0，0.65，0.55；0.25，0.05，0.10；0.25，0.05，0； %输入系数矩阵Y=50000，25000，0； %输入外界需求X（eye（size（A）-A）Y %求总产出 >>C=A*diag（X）%求投入产出矩阵 >>D=ones（size（X）*C %求总投入量 >>Y=80000，8000，800；%输入最终需求 >>X=（eye（s

15、ize（A）-A）；Y %求未来产出向量练习题 1、解线性方程组 （1） （2） 2、某机床厂的产品数量与成本的数据如下： 产品数量（台）51015成本（千元）100155220求其成本函数，盈亏转折点及最大利润值。 3、在某经济年度内，各经济部门的投入产出表如下： 单位：（10亿元）工业农业第三产业最后需求（）总产值（）工业6211625农业15第三产业31520假设经济年度工业，农业及第三产业的最后需求均为17（10亿元），预测经济年度工业，农业及第三产业的产出。基础实验四 矩阵与几何变换一、实验目的理解矩阵的一些性质和特征、矩阵与几何变换的关系，了解代数与几何的联系。二、实验材料几何变换

16、是计算机辅助图形设计中的基本技术，有着广泛的应用。为了在形式上将几何变换统一表示，引入了齐次坐标和齐次变换矩阵。二维直角坐标系中，直角坐标为的点的齐次坐标表示为，其中是不为零的常数，一个点的齐次坐标不是唯一的。齐次坐标和代表同一个点。二维直角坐标系中的齐次变换矩阵形式为，通常将其分块为，其中，。几何变换和矩阵的对应关系一个几何变换对应一个形为的齐次变换矩阵；反之，一个这种形状的可逆矩阵对应一个几何变换。矩阵的逆与它对应的几何变换的逆变换相对应，如旋转角度为的变换对应的齐次变换矩阵为其逆变换对应的齐次变换矩阵恰为此变换矩阵的逆矩阵平移量为的平移变换、比例系数与的比例变换、比例系数的整体比例变换、

17、错切系数为的向的错切变换、错切系数为的向的错切变换、关于轴对称的变换、关于直线对称的变换的齐次变换矩阵分别为、。2.3可逆矩阵的初等分解和几何变换的初等分解 一个形为的可逆矩阵可以表示为一系列初等矩阵的乘积；一个几何变换可以表示为一系列简单几何变换的乘积。 常见的简单几何变换有：平移量为的平移变换，平移量为的平移变换，向比例系数为的比例变换，向比例系数为的比例变换，错切系数为的向的错切变换，错切系数为的向的错切变换，关于轴对称的变换，关于直线对称的变换。练习1 几何变换能改变图形的形状、大小和方位，象比例放缩、平移、对称（亦称镜象）等，是计算机辅助图形设开的基本方法和工具。在二维直角坐标系中，

18、试用一个矩阵实现关于直线对称的齐次变换。该变换等于以下五个变换的乘积：平移量为的平移变换，旋转角度为的旋转变换，关于轴对称的变换，的逆变换，的逆变换。从而所求的矩阵为。Mathematica程序如下：a=- ArcTan3；T1=1，0，0，0，1，0，0，- 5，l；T2=Cosa，Sina，0，- Sina，Cosa，0，0，0，l；T3=1，0，0，0，- 1，0，0，0，1；T=NT1。T2。T3。InverseT2。InverseT1练习2 几何变换也是生成计算机动画的基本手段和方法。在二维直角坐标系中，中心坐标为( 0 , 0 ) ,半径为1的圆绕点旋转一周，要求每隔画一个圆，并给

19、出动画演示。绕点旋转角的齐次变换矩阵为Mathematica程序如下：x0=0；y0= 0；r=1；n= 30；xp=2；yp=2；a=2Pin；g3= ；g1= ParametricPlot xp+2*Sqrt2Cost，yp+2Sqrt2Sint,t，0，2Pi,AspectRatio->1，PlotRange->-3，6，-3，6，Plotstyle->RGBColor1，0，0，DisplayFunctity> Identity；T11，0，0，0，1，0，-xp，-yp，l；Fori=0，in,i+， T2=Cosa*i，Sina*i，0，-Sina*i，Co

20、sa*i，0，0，0，1；T； xc= Partx0，y0，1.TT，1； yc=Partx0，y0，l。TT，2； xt_=xc+r*Cost; yt_= yc+ r*Sint； g2=ParametricPlotxt，y t，t，0，2*Pi，AspectRatio->1， PlotRange->-3，6，-3，6，DisplayFunction-> Identity；g3=Appendg3，GraphicsDiskxc，yc，0.1； Showg1，g 2，g 3，DisplayFunction->DisplayFunction用Mathematica的动画播放功

21、能，播放这一系列图形，以理解旋转变换在此的作用，并注意它的用法。练习 3 证明矩阵能代表一个几何变换，并将其进行初等分解。能否代表一个几何变换，取决于这一矩阵是否可逆。求矩阵的行列式值的Mathematica程序T=1，2，0，3，4，0,5，6，1；DetT结果为。可见，矩阵是可逆的，能代表一个几何变换。又利用矩阵的初等分解方法可知，可以进行如下把变换为单位矩阵的过程，Mathematica程序为：T1=1，0，0，-3，1，0，0，0，1；（矩阵的第一行乘以-3加到第二行）T2=l，0，0，0，1，0，-5，0，1； T3=l，0，0，0，-12，0，0，0，1； T4=1，- 2，0，0

22、，l，0，0，0，1； T5=1，0，0，0，1，0，0，4，l；T5* T4* T3* T2* T1*T （矩阵变换为单位矩阵）初等行变换矩阵、分别左乘矩阵得单位矩阵。所以有 矩阵代表的几何变换可分解出五个简单变换的乘积。练习4 几何变换还是产生新图形的重要方法。令正方形围绕其中心旋转且放大，得到一组正方形。使得每个后继的正方形恰包含前一个，这样画出过程中的每个正方形，便形成了正方形的螺旋线。设初始正方形的四个顶点坐标分别为,每次旋转的角度为。很明显，此问题包含旋转和比例两个变换，每次旋转的角度为，比例系数经过简单推导可得为，又相应于这两个变换的齐次变换矩阵为,而且后一个正方形的各端点等于前

23、一个正方形相应的端点依次经过上述的旋转和比例变换后的结果。根据以上分析，可编出求解Mathematica程序如下： x10=-1；y10=-1；x20=1；y20 = -1；x30=l；y30=1；x40=-1；y40=1；a=Pi30；s Sina+Cosa；T1=Cosa，Sina，0，-Sina，Cosa，0，0，0，l；T2=s，0，0，0，s，0，0，0，1；TT=NT1。T 2G1=；Fori=1，i < 2 0,i+， v1=Nx10，y1 0，1。TT；v2=N x20，y20，1。TT；v3=Nx30，y30，l。TT；v4=N Fx40，y40，1。TT； xl= P

24、artv1，l；y1=Partvl，2；x2=Partv2，1；y2=Partv2，2； x3=Partv3，1；y3=Partv3，2；x4=Partv4，1；y4=Partv4，2； g1= Appendg1，GraphicsLinexl，y1，x2，y2，x3，y3，x4，y4，x1，y1，AspectRatio->1，PlotRange->-10，10，-10，10；x10=x1；y10=y1；x20=x2；y20=y2；x30=x3；y30=y3；x40=x4；y40=y4；Show g1输入并运行此程序，以理解这里的实现方法。思考题1、求二维直角坐标系中，绕点旋转的齐次

25、变换矩阵，其中点的齐次坐标为。2、求二维直角坐标系的中，关于直线对称的齐次变换矩阵。3、在二维直角坐标系中，中心坐标为，长短半轴长分别为2，1的椭圆，绕点旋转一周，要求每隔画一个椭圆，并给出动画演示。4、证明矩阵能代表一个几何变换，并将其进行初等分解。5、做出正六边形的螺旋线。提示：比例系数为，其中为每次旋转的角度，为正六边形的内角。基础实验五 数据拟合与曲线拟合一、实验目的 对于某个变化过程中的相互依赖的变量，可建立适当的数学模型，用于分析、预报、决策或控制该过程。对于两个变量可通过用一个一元函数去模拟这两个变量的取值，但用不同的方法可得到不同的模拟函数。 使用最小二乘法来进行数据拟合，用基

26、本函数曲线及其变化模拟给定的曲线，理解拟合方法。二、实验材料2.1曲线拟合 （1）初等函数包括基本初等函数与它们经过加减乘除复合等运算后所得到的函数的图形及其变换。拟合函数为多项式情形理论上已经解决，称为拉格朗日插值多项式。（2）光滑曲线的有关内容，包括分段函数的连续性、一阶可导性与高阶可导性。（3）方程或方程组的求解，包括超越方程或方程组的近似解法，线性方程组的精确解。2.2最小二乘法给定平面上一组点(,)()作曲线拟合有多种方法，其中最小二乘法是常用的一种。最小二乘法的原理是：求,使达到最小。拟合时，选取一定的拟合函数形式，设拟合函数的基底函数为拟合函数为 确定使方差达到极小，此时得到的即

27、为所求。为使取到极值，将的表达式代入，对求的偏导数，令其等于零，得到方程组成的方程组，从中求解。当=1时，取拟合函数，此做法称为线性拟合，统计学上叫做线性回归。此时，临界方程组为从中解出与，有，其中，, 。Mathematica提供了最基本的数据拟合函数Fit，这个函数使用最小二乘法产生基函数的线性组合以构造出拟合函数。函数的参数表中包括三项：第一个参数是被拟合的数据；第二个参数是一个表，用于说明拟合用的基函数；第三个参数是拟合变量。2.3 线性拟合练习1 为研究某一化学反应过程中温度对产品得率(%)的影响,测得数据如下：100110120130140150160170180190455154

28、61667074788589试求其线性拟合曲线。 Mathematica程序:b1=100,45,110,51,120,54,130,61,140,66,150,70,160,74,170,78, 180,85,190,89 （将数据以表的形式输入）ft1=Fitb1,1,x,x （用Fit拟合，这里是线性拟合）gp=Plotft1,x,100,190,PlotStyle->RGBColor1,0,0 （作拟合曲线的图形）fp=ListPlotb1,PlotStyle->PointSize0.05,RGBColor0,0,1 （作散点图）Showfp,gp （显示点组与拟合曲线，作

29、图。下面为计算残差的程序）a= ;b= ; （a，b的值由上面的结果确定）fx_=a*x+b; （拟合函数）darata=Sum(b1i,2-fb1i,1)2,i,1,10(计算残差)改变b1去掉b1中的一组数据100，45，求改变后的线性拟合曲线去掉b1中的一组数据110，51，求改变后的线性拟合曲线去掉b1中的一组数据120，54，求改变后的线性拟合曲线去掉b1中的一组数据130，61，求改变后的线性拟合曲线去掉b1中的一组数据140，66，求改变后的线性拟合曲线去掉b1中的一组数据150，70，求改变后的线性拟合曲线去掉b1中的一组数据160，74，求改变后的线性拟合曲线去掉b1中的一组

30、数据170，78，求改变后的线性拟合曲线去掉b1中的一组数据180，85，求改变后的线性拟合曲线去掉b1中的一组数据190，89，求改变后的线性拟合曲线去掉b1中的两组组数据100，45，110，51求改变后的线性拟合曲线去掉b1中的两组组数据100，45，120，54求改变后的线性拟合曲线去掉b1中的两组组数据100，45，130，61求改变后的线性拟合曲线去掉b1中的两组组数据100，45，140，66求改变后的线性拟合曲线去掉b1中的两组组数据100，45，150，70求改变后的线性拟合曲线去掉b1中的两组组数据100，45，160，74求改变后的线性拟合曲线去掉b1中的三组组数据100

31、，45，110，51，190，89求改变后的线性拟合曲线去掉b1中的三组组数据100，45，160，74，180，85求改变后的线性拟合曲线去掉b1中的四组组数据110，51，160，74，170，78求改变后的线性拟合曲线去掉b1中的四组组数据120，54，150，70，190，89求改变后的线性拟合曲线增加b1中的一组数据125，58，求改变后的线性拟合曲线2.4 非线性拟合练习2 在某一化学反应里，由实验得到生物的浓度与时间（分）的关系如下123456789101112131415164.08.010.0求浓度与时间关系的拟合曲线。提示：先用ListPlot语句描点，观察点的分布情况，以

32、确定拟合函数。(1)用多项式函数拟合的Mathematica程序:Cleargp,fp;b2=1,4,2,6.4,3,8.0,4,8.4,5,9.28,6,9.5,7,9.7,8,9.86,9,10.0,10,10.2,11,10.32,12,10.42,13,10.5,14,10.55,15,10.58,16,10.6gp=ListPlotb3,PlotStyle->RGBColor0,1,0,PointSize0.04 ft2=Fitb3,Tablexi,i,0,4,x （用四次曲线拟合） fp=Plotft2,x,0,17,PlotStyle->RGBColor1,0,0 S

33、howgp,fp fx_=expr; （用拟合的多项式函数来定义f(x)） darata=Sum(b2i,2-fb2i,1)2,i,1,16(计算残差) (2)用函数作拟合，求拟合曲线。作变换,，拟合函数变形为。Mathematica程序为:fxx_:=1/xfyy_:=Logy nb=Tablefxb2i,1,fyb2i,2,i,1,16ft3=Fitnb,1,x,x（拟合）f4=a*Expb/x （a，b的值由上面的结果确定） t1=Plotf4,x,1,18,PlotStyle->RGBColor1,0,0t2=ListPlotb2,PlotStyle->RGBColor0,

34、1,0,PointSize0.05Show%,%(3)用作拟合。Mathematica程序为:gy_:=1/ysb=Tableb2i,1,gb2i,2,i,1,16ft5=Fitsb,1,1/x,x f5=1/ft5t3=Plotf5,x,1,16,PlotStyle->RGBColor0,0,1 Showt1,t2,t3（在一张图上比较一下用两种方法得到的函数曲线） （4）用分段函数作拟合。2.4 思考题 1、在钢线碳含量（）对于电阻（）的效应的研究中，得到以下数据：1518192126求其线性拟合曲线.2、一种合金在某种添加剂的不同浓度（）之下作抗压强度（）实验，得数据如下：以模型作

35、曲线拟合。3、下表是某年美国轿车价格的调查资料, 以模型作曲线拟合:使用年数12345678910平均价格2615194314941087765538484290226204 4、待拟合数据如下234567891011121314151610 试作非线性拟合。第二章 探索数学实验探索实验一 素数一、实验背景与实验目的 德国数学家高斯说过，数学是科学的女王，而数论则是数学的女王。在数论这一充满了趣味而布满荆棘的领域中，有关素数的问题(如著名的Goldbach猜想)始终是最富有魅力最吸引人的研究问题。本实验将探索素数的规律及其相关的某些有趣问题。具体地，请同学们研究下列素数问题：（1）素数表的构造

36、；（2）素数的判别；（3）最大的素数；（4）构造生成素数的公式；（5）素数的分布。我们希望通过本节的学习激发同学们对数论的好奇心，让同学们被自然数的神奇规律而折服，同时让同学们认识到探索自然数规律的艰难性。二、实验材料2.1素数的判别与求解如果一个大于1的自然数只能被l及它本身整除，则该数称为素数，否则称为合数。从数学史的黎明时期开始，数学家们就一直在探索自然数的奥秘。远在古希腊时代，欧几里得就证明了每一个合数都可以分解为若干个素数的乘积，并且在不计较素数排列顺序时这种分解是唯一的，这就是所谓的算术基本定理。算术基本定理表明，素数是构造自然数的基石，正如物质的基本粒子一样。正是由于素数如此重要

37、的地位才使得一代又一代数学家努力地探询素数的规律。首先，一个最基本的问题是素数到底有多少个?会不会在某一充分大的自然数以后就没有素数呢?建议同学们先做以下实验：假设已知前个素数，它们按从小到大的顺序排列为。对，计算。问：1)是否都是素数? 2)如果不是素数，是否含有不同于，的素因子?Mathematica 程序：NumPn_Integer:= Modulei,Num, Num=ProductPrimei,i,1,n+1; PrintNum; PrintPrimeQ Num; PrintFactorIntegerNum DoNumPn,n,1,20根据实验结果，猜测素数是否有无穷多个?试给出证明

38、。关于素数的下一个基本问题是:如何求出小于某一给定整数的所有素数?古希腊另一位学者Eratosthenes给出了解决这一问题的方法，这一方法被后人称为Eratosthenes筛法。Eratosthenes筛法的基本思想是，将自然数列从2开始按顺序排列至某一整数。首先，从上述数列中划除所有2的倍数(不包括2)。在剩下的数中，除2外最小的是3。接着，从数列中划掉所有3的倍数(不包括3)。然后在剩下的数中，再划去5的倍数。这个过程一直进行下去，则最后剩下的数就是不超过的所有素数。借用Eratosthenes筛法，经过众多学者的艰辛努力，D.N.Iehmer于1914年编织出了10 000 000以内

39、的素数表。利用Eratosthenes筛法，请同学们手工编写100以内的素数表。据此，思考如何估计利用筛法编写10000000以内的素数表的工作量?利用Eratosthenes筛法，通过计算机编程求10000以内的所有素数。Mathematica 程序：Sieven_Integer:= Module t=,i,temp, Fori=2,i<=n,i+,AppendTot,i; Fori=1,Primei<=Sqrtn,i+,temp=Primei; t=Selectt,(#1=temp|Mod#1,temp!=0)& tSieve1000筛法是用乘法寻找素数。实际上，也可以

40、用除法判别一个数是否是素数。而且，用除法的效率可能会更高。假设我们已经找到了前个素数，为了寻找下一个素数我们从开始依次检验每一个整数，看是否能被某个，整除。如果能被前面的某个素数整除，则为合数。否则即为下一个素数，实际上，为了提高算法的效率，人们不需用前面的每一个素数去试除，而只需用不超过的素数去除就可以了。利用试除方法，通过计算机编程求10000以内的所有素数。Mathematica 程序：DivPrimen_Integer:= Module t=,i,j,temp,divided, Fori=2,i<=n,i+, j=1;divided=False; WhilePrimej<=

41、Sqrti&&(!divided), temp=Primej; divided=(Modi,temp=0);j=j+1; If!divided,AppendTot,i ; tDivPrime10000试将试除法与筛法进行比较，哪一个更有效?在以上Mathematica 程序后分别加上程序TimingSieve10000，TimingDivPrime10000。虽然从理论上来说，Eratosthenes筛法和试除方法可以求出所有的素数，但是通过上面的实验可能会发现，利用这些方法构造大的素数表是不切实际的。例如，要构造以内的素数表，利用当今最快的计算机也得需要近大约一亿年。实际上，

42、数学工作者一直在致力扩大素数表的范围。寻找更大的素数不仅是数学爱好者的乐趣，更是实际应用工作者(如密码编码者)努力追寻的事。数十年前，象l l1(23个1)以及的素性判别问题难倒了许多睿智的数学家。当今，数学家研究出了一套非常复杂而高深的技巧来对数的素性进行检验。寻找到的最大素数的高峰一股劲地向上无限攀升。截止到2008年，所找到的最大素数是，这是第 46个 梅森(Mersenne)素数，其十进制形式有1300万位！为了领略一下隐藏在这些高深技巧下面的数学思想，请同学们做以下实验。对，观察被整除所得的余数。从观察结果你能得出什么结论?再取其它的整数 (如3，4，5)，观察被整除的情况。特别注意

43、观察当为素数时的结果。能否因此确信你的结论?进一步，你所得出的结论的逆命题是否成立?由此，用你的结论能否给出判别一个数是否是素数的判别方法? 对互质的整数及，求使得除的余数为1的最小整数。当视为素数时，观察与之间的关系，你能得到什么结论?类似地，对做进一步的观察，你能否确信你的结论?你所得出的结论的逆命题是否成立?Mathematica 程序： Mn_Integer:=Moduley,k,m=2; k=m(n-1);x=Modk,n; Printn," ",PrimeQn," ",x," ",GCDm,n DoMn,n,2,200上述

44、实验表明，给出一个简明的素数判别准则并不容易，通常需要将更多高深的技巧与之结合才能给出判别素数的更有效方法。对于具有特殊结构的数的素性判别，有更加快捷的方法。这方面最引人注目的例子是Mersenne数，即形如的数。利用Mersenne数可以构造出非常大的素数，如前面指出的最大素数。对，判断哪些Mersenne数是素数?如果为合数，Mersenne数是素数还是合数?如果为素数，Mersenne数是否一定是素数?Mersenne数素性判别的Mathematica 程序：Mersennen_Integer:= Module M,i,u=4, If!Prime Qn, False, M=2n-1; F

45、ori=1,i<n-1,i+,u=Modu2-2,M;, Ifu=0,True,False 从实验中可以看到，Mersenne素数是极其稀少的。借助大型计算机的威力，截止2008年，数学家仅发现了46个Mersenne素数。可以看出，Mersenne素数的分布是极不规则的！如何通过判断是否是素数?数学家Lucas与Lehmer早在1930年就给出了具体可行的判别方法。假定为素数，定义数列，。如果，则从为素数；否则，从为合数。前面列出的Mersenne素数就是在此算法的基础上发现的。有人猜测，如果本身是Mersenne素数，则此数也一定是Mersenne素数。你是否相信这一猜测?还有人相信

46、Mersenne素数有无穷多个，你的意见呢?2.2生成素数的公式虽然我们在本段给出了生成素数表的算法，但要找到大的素数确是异常的艰辛。你也许会想，要是能找到一个正好生成全部素数的公式就好了。例如，是否存在单变量整系数的多项式，它只生成素数并且生成所有的素数?更一般地，是否存在一个生成素数的多变量函数公式?如果这样的公式不存在，能否找到一个虽不能给出全部但能给出无穷多个素数(且只给出素数)的公式?实际上，早在十七、十八世纪，大数学家Fermat，Euler等就研究过这类公式。1640年Fermat在给Mersenne的信中指出，对所有的整数，永远是素数。的确，都是素数。然而，好事到此为止。173

47、2年，大数学家Euler指出，不是素数，他并且找到了的因子分解。此后，人们分别证明了，与都是合数，并得到了它们的素因子分解。实际上有人猜测当时都是合数。由此可见，Fermat当初的猜测离真实结果相差有多远。Fermat数与正多边形做图有紧密的联系。这一惊人的发现是伟大的德国数学家Gauss得到的。Gauss在19岁那年证明了：一个正边形可用直尺与圆规作图的充要条件是，或者，其中为不同的Fermat数。特别地，正17边形可以用直尺与圆规做出。1772年Euler研究过用单变量的整系数多项式来生成素数，其最著名的例子是。对，计算。它们是否都给出素数?在l0 000以内的素数中，由公式给出的素数占多

48、少?类似地，对公式及做同样的判别。你自己能否给出一个类似的公式? Mathematica 程序： t=0; Eun_Integer := Moduley ,y= n2 + n +41; IfPrimeQyTrue, t=t+1; DoEun ,n,1,100 Printt公式有一个非常奇特的性质。为揭示这一特性，考察它的二次求根公式的判别式。163有什么特别的地方?有！请看三个互不相干的无理数产生了一个非常接近整数的数，确实让人吃惊。你也许会想这不是一个巧合，幕后一定隐藏着什么。你确是猜对了，它与Gauss关于类数问题的工作紧密相关。实际上，是使 (其中都是整数)这样的数系具有唯一因子分解的最

49、大的值。你认为是否存在单变量的整系数多项式，使得对每个都是素数?你能否证明你的结论?从上面的努力似乎可以看出，利用单变量函数难以得到素数的生成公式。那么，多变量函数是否可行呢?作为Hilbet第十问题的一个推论，马蒂雅舍维奇证明了：存在一个多元多项式，其正值构成的集合恰好是素数的全体。遗憾的是，他并没有给出怎样具体地构造这样的多项式。后经众多数学家的努力，终于在1977年构造出了一个具有26个变量25次的素数生成多项式！2.3 素数的分布如果将素数在数轴标出来，同学们会发现素数的分布很不规则。通过以下实验做进一步的观察。（1）用表示不超过的素数的个数，表示区间内素数的个数。试计算，以及，。从计

50、算结果看，随着整数范围的扩大，者数是越来越稀还是越来越密？选取一些更长的区间，再尝试以上同样的实验。（2）将素数从小到大顺序排列，用表示相邻素数间的间隔。计算(如1000，10000等)，然后将点标在平面坐标系中。你能从中观察到素数的间隔有什么规律呜?譬如，素数的间隔值有哪些?它们各重复多少次?哪些间隔值的重复次数多?最大间隔值是多少?随增大，最大间隔值是否也随之增大?根据上述实验，你对素数的分布能做什么样的猜测?譬如，间隔差为2的素数是否有无穷多个?更一般地，间隔差为某个偶数的素数对是否有无穷多个?是否存在相邻的素数，其间隔值可以任意大?你能否证明你的这些猜测?从上面的实验了解到，素数的分布是极不规则的。它虽然延数轴分布越来越稀疏，但有时素数之间的间隔又很小。不过，就总的趋势而言，固定区间长度内的素数的个数是越来越少。下面，我们进一步研究素数个数的变化规律。在二维坐标平面上标出点列， (取不同的，如1000，10 000等)，也可以用折线将点