南京邮电大学 数值计算实践报告

1、数值计算实践I、方程求根1、 实验目的熟悉和掌握Newton法,割线法,抛物线法的方法思路，并能够在matlab上编程实现2、 问题描述(1).给定一个三次方程,分别用Newton法,割线法,抛物线法求解.方程的构造方法:(a)根:方程的根为学号的后三位乘以倒数第二位加1再除以1000.假设你的学号为B06060141,则根为141*(4+1)/1000=0.564(b)方程:以你的学号的后三位数分别作为方程的三次项,二次项,一次项的系数,根据所给的根以及三个系数确定常数项.例如:你的学号是B06060141,则你的方程是x3+4x2+x+a0=0的形式.方程的根为0.564,因此有0.564

2、3+4*0.5642+0.564+a0=0,于是a0=-2.015790144你的方程为x3+4x2+x-2.015790144=0.(2)假设方程是sinx+4x2+x+a0=0的形式（三个系数分别是学号中的数字），重新解决类似的问题(3)构造一个五次方程完成上面的工作.四次方程的构造:将三次多项式再乘以(x-p*)2得到对应的五次多项式(p*为已经确定的方程的根,显然,得到的五次方程有重根).(4)将（2）中的方程同样乘以(x-p*)得到一个新的方程来求解注:(1)Newton法取0.5为初值,割线法以 0,1为初值,抛物线法以0,0.5,1为初值,(2)计算精度尽量地取高.终止准则：根据

3、来终止 (3)可供研究的问题：（一）的取值不同对收敛速度有多大的影响（二）将注（1）中的初值该为其它的初值，对收敛性以及收敛速度有无影响（三）)能否求出方程的所有的根（4）实验报告的撰写 实验报告包含的内容：（一）实验目的（二）问题描述（三）算法介绍（包括基本原理）（四）程序（五）计算结果（六）结果分析（七）心得体会3、 算法介绍在本问题中，我们用到了newton法，割线法，抛物线法。1.Newton法迭代格式为： 当初值与真解足够靠近，newton迭代法收敛，对于单根，newton收敛速度很快，对于重根，收敛较慢。2.割线法：为了回避导数值的计算，使用上的差商代替，得到割线法迭代公式：割线法

4、的收敛阶虽然低于newton法，但迭代以此只需计算一次函数值，不需计算其导数，所以效率高，实际问题中经常应用。3. 抛物线法：可以通过三点做一条抛物线，产生迭代序列的方法称为抛物线法。其迭代公式为：其中 是一阶均差和二阶均差。收敛速度比割线法更接近于newton法。对于本问题的解决就以上述理论为依据。终止准则为：本题中所有精度取1e-8。4、 程序计算结果问题一根据所给的要求，可知待求的方程为：牛顿法建立newton_1.m的源程序，源程序代码为：function y=newton_1(a,n,x0,nn,eps1) x(1)=x0; b=1; i=1; while(abs(b)eps1*x(

5、i) i=i+1; x(i)=x(i-1)-n_f(a,n,x(i-1)/n_df(a,n,x(i-1); b=x(i)-x(i-1); if(inn)error return; end end y=x(i);建立n_f.m的源程序来求待求根的实数代数方程的函数，源程序代码为：function y=n_f(a,n,x)% 待求根的实数代数方程的函数 y=0.0; for i=1:(n+1) y=y+a(i)*x(n+1-i); end建立n_df.m的源程序来方程的一阶导数的函数，源程序代码为：function y=n_df(a,n,x)% 方程的一阶导数的函数 y=0.0; for i=1:

6、n y=y+a(i)*(n+1-i)*x(n-i); end在matlab软件中执行下列语句并得到的最终结果截图：割线法建立gexian.m的源程序，源程序代码为function x=gexian(f,x0,x1,e) if nargine i=i+1; z=x-(feval(f,x)*(x-y)/(feval(f,x)-feval(f,y); y=x; x=z; end i 在matlab软件中执行下列语句并得到的最终结果截图：抛物线建立paowuxian.m的源程序，源程序代码为：function x=paowuxian(f,x0,x1,x2,e) if nargine i=i+1; h1

7、=y-z; h2=x-y; c1=(feval(f,y)-feval(f,z)/h1; c2=(feval(f,x)-feval(f,y)/h2; d=(c1-c2)/(h1+h2); w=c2+h2*d; xi=x-(2*feval(f,x)/(w+(w/abs(w)*sqrt(w2-4*feval(f,x)*d); z=y; y=x; x=xi; endi在matlab软件中执行下列语句并得到的最终结果截图：研究一：只改变初值由上述结果可知，方程的解在0.2附近，所以将牛顿法为0.2；割线法的初值设为0，0.4；抛物线法的初值设为0，0.2，0.4；牛顿法根据问题1中牛顿法的程序，在mat

8、lab软件中执行下列语句并得到的最终结果截图：割线法根据问题1中割线程序，在matlab软件中执行下列语句并得到的最终结果截图：抛物线法根据问题1中抛物线法程序，在matlab软件中执行下列语句并得到的最终结果截图：研究二 只改变精度将精度由1e-8改为1e-50和1e-100观察迭代次数有何变化牛顿法：根据问题1中的牛顿法的程序，在matlab软件中执行下列语句并得到的最终结果截图：精度为1e-50时精度为1e-100时割线法根据问题1中的割线法的程序，在matlab软件中执行下列语句并得到的最终结果截图：精度为1e-50时精度为1e-100时抛物线法根据问题1中的抛物线法的程序，在matl

9、ab软件中执行下列语句并得到的最终结果截图：精度为1e-50时精度为1e-100时、研究结论在只改变初值时，当初值定得越靠近初值，迭代次数就越少。在只改变精度时，当精度越来越大时，迭代次数并几乎不变。综上所述，初值对迭代次数的影响比较大，精度对迭代次数影响不大。问题二问题描述根据所给的要求，可知待求的方程为：问题分析仍然利用（1）中方法求解这一问题，并利用图解法找到初值，通过观察图像，将newton法初值设为：0.1，割线法初值设为：0,0.2。抛物线法初值设为：0,0.1,0.2。图像见下图：Newton法问题一的牛顿法的求解只适用于线性方程，所以在问题二中用其他方法来求解方程。建立newt

10、on1.m源程序，源程序代码为：function x=newton1(fn,dfn,x0,e)if narginei=i+1；x0=x;x=x0-feval(fn,x0)/feval(dfn,x0);end在matlab软件中执行下列语句并得到最终结果截图割线法利用问题一中的割线法程序，在matlab软件中执行下列语句抛物线法利用问题一中的抛物线法程序，在matlab软件中执行下列语句问题三问题描述按照题目要求对五次方程进行构造为：问题分析仍然利用一中方法求解这一问题，并利用图解法找到初值，通过观察图像，将牛顿法的两组初值为0以及0.5，割线法初值设为：0,0.5以及0，0.2。抛物线法初值设

11、为：两组初值为0,0.5,1以及0,0.25,0.5。牛顿法利用问题二中的程序，在matlab软件中执行下列语句：初值为0时初值为0.5时割线法利用问题一中割线法的程序，在matlab软件中执行下列语句：初值为0，0.5时初值为0，0.2时抛物线法利用问题一中抛物线法的程序，在matlab软件中执行下列语句：初值为0，0.5，1时初值为0，0.25，0.5时问题四问题描述根据题目要求对方程进行构造为：问题分析仍然利用问题一中方法求解这一问题，并利用图解法找到初值，通过观察图像，newton法初值选取了两组初值为0以及0.5，割线法初值设为：0，0.5和0，0.3。抛物线法初值设为：两组初值为0

12、,0.2,0.4以及0，0.5，1。牛顿法利用问题二中的程序，在matlab软件中执行下列语句：初值为0时初值为0.5时割线法利用问题一中割线法的程序，在matlab软件中执行下列语句：初值为0，0.5时初值为0，0.3时抛物线法利用问题一中抛物线法的程序，在matlab软件中执行下列语句：初值为0，0.1，0.2时初值为0，0.5，1时5、 计算结果及分析Newton法割线法抛物线法问题一0.2240.2240.224问题二0.14660.14660.1466问题三0.2240.2240.224问题四0.224和0.1465772585571010.224和0.14660.224和0.146

13、6 问题一 迭代步数 问题二 迭代步数 问题三迭代步数 问题四迭代步数Newton法 6 5 37 8 9割线法 8 5 50 13 10抛物线法 9 5 63 14 13结果分析将Newton法，割线法，抛物线法进行比较可以看到在本文题中，三种方法计算得到的最终结果基本相同，但是迭代步数有较大差别，综合看来Newton法迭代步数最少，割线法次之，抛物线法最次。在各个问题的研究中，我通常都会采用不同的初值，发现不同初值会对应不同的迭代次数，另外针对问题一，我选用了不同的精度，发现迭代次数并没有很大的变化，因而一个初值的选定可以对迭代次数产生很大的影响，而精度的改变对迭代次数的影响很小。对每个算

14、法单独来看，显然选择初值不同对于迭代步数影响较大，对于找到根也会有影响。因此应该先通过画图确定根的大致位置，给出在其附近的初值。6、 心得体会在实现这三个算法的过程中，本身编程较易实现，最重要的是对算法本身的理解，只有真正理解算法的含义才能更快更好的实现程序。II、离散型最小二乘和连续型最小二乘问题一、实验目的掌握并能够利用离散型最小二乘和连续性最小二乘求解问题。二、问题描述1：以函数生成的6个数据点（0.25,23.1）,（1.0,1.68）, （1.5,1.0）,（2.0,0.84）,（2.4,0.826）,（5.0,1.2576）为例，运行程序得到不同阶对应的曲线拟合产生的多项式函数。2

15、. 例题：计算f(x)=exp(x)在-1，1上的二、三次最佳平方逼近多项式。并画图进行比较。三、问题分析问题1是离散最小二乘问题。最小离散最小二乘就是根据一批有误的数据（）i=1,2，n确定参数，并通过均方误差来比较曲线拟合的优劣，在本题中通过画图来比较不同阶方程拟合效果的优劣。选择两种方法实现离散最小二乘。方法一，建立normal equation(法方程组)，求解k次多项式系数。法方程组构造方法： 方法二：由于在matlab中存在ployfit函数，可以即为方便的用k次多项式拟合。问题2是一个连续型最小二乘法的应用实例。对于给定的函数，如果存在使得则称S*(x)是f (x)在集合中的最佳

16、平方逼近函数。显然，求最佳平方逼近函数的问题可归结为求它的系数，使多元函数取得极小值，也即点()是I (a0, ，an)的极点。由于I (a0, a1, ，an)是关于a0, a1, ，an的二次函数，利用多元函数取得极值的必要条件， (k = 0, 1, 2, , n)即得方程组如采用函数内积记号那么，方程组可以简写为.（1）这是一个包含n + 1个未知元a0, a1, , an的n + 1阶线性代数方程组，写成矩阵形式为 （2）此方程组叫做求aj (j = 0, 1, 2, , n)的法方程组。显然，其系数行列式就是克莱姆行列式Gn = Gn (j0, j1, , jn)。由于j0, j1

17、, , jn线性无关，故Gn 0，于是上述方程组存在唯一解。从而肯定了函数f (x)在中如果存在最佳平方逼近函数，则必是 .（3）四、实验程序问题一法一：离散最小二乘法function f=normalequation(n)syms tr=0;xx=0.25,1.0,1.5,2.0,2.4,5.0;yy=23.1,1.68,1.0,0.84,0.826,1.2576;x=zeros(n,n);y=zeros(n,1);for i=1:n for k=1:6 x(i,i)=xx(k).(2*i-2)+x(i,i); y(i,1)=yy(k).*xx(k).(i-1)+y(i,1); endend

18、for i=2:n for j=1:i-1 for k=1:6 x(i,j)=xx(k).(i+j-2)+x(i,j); x(j,i)=x(i,j); end endendz=xy;for i=1:n r=z(i,1)*t(i-1)+r;endvpa(r,5)法二:用matlab中的ployfit函数对k次多项式进行拟合。建立number2.m文件，带入如下：x=0.25,1.0,1.5,2.0,2.4,5.0;y=23.1,1.68,1.0,0.84,0.826,1.2576;plot(x,y,o)p1=polyfit(x,y,1) %利用线性拟合xi=-0.25:0.01:6.0;y1=p

19、olyval(p1,xi);plot(x,y,o,xi,y1,r:);hold on;p2=polyfit(x,y,2) %二次拟合y2=polyval(p2,xi);plot(x,y,o,xi,y2,m);hold on;p3=polyfit(x,y,3) %三次拟合y3=polyval(p3,xi);plot(x,y,o,xi,y3,b:);hold on;p4=polyfit(x,y,4) %四次拟合y4=polyval(p4,xi);plot(x,y,o,xi,y4,c);hold on;p5=polyfit(x,y,5) %五次拟合y5=polyval(p5,xi);plot(x,y

20、,o,xi,y5,g);hold off;legend(初始点,y=-2.3840*x+9.8499,初始点,y=2.1793*x2-15.8845*x+22.3999,初始点,y=-2.0143*x3+18.3499*x2-45.2167*x+32.7386,初始点,y=1.4828*x4-16.0435*x3+56.3154*x2-79.4994*x+39.6688,初始点,y=-0.9961*x5+12.1743*x4-55.2367*x3+116.6262*x2-116.7266*x+45.8090);问题二：最佳平方逼近算法（1）输入被逼近函数f(x)和对应的逼近区间a,b并选择逼近

21、函数系和权函数；（2）解方程组（1）或（2），其中方程组的系数矩阵和右端的项由式（3）得到；（3）由式（3）得到函数的最佳平方逼近。将上述算法编写成MATLAB程序共需三个程序：第一个程序（函数名Sapproach.m）计算最佳逼近函数的系数：function S=Sapproach(a,b,n) %定义逼近函数global i;global j; if nargin fun=exp(x); fplot(fun,-1,1) hold on xi=-1:0.1:1; yi=polyval(S,xi); plot(xi,yi,r:)得到以下结果：当三次逼近时，有以下结果绘制图形如下： fun=ex

22、p(x); fplot(fun,-1,1) hold on xi=-1:0.1:1; yi=polyval(S,xi); plot(xi,yi,r:)得到如下所得图形：六、结果分析显然离散最小二乘中次数较高拟合的效果较好，但由于本问题中初始点较少，并不能完全反应函数本身的特点，同时在本问题中，选择了两种方式，结果接近，也可以互相验证正确性。在连续最小二乘法的实验中较顺利的达到了预期的结果。从试验结果看出三次逼近没有二次逼近效果理想，验证了最佳平方逼近理论。七、心得体会 在离散最小二乘与连续最小二乘程序设计的过程中，要么通过均差来表示函数拟合的优劣，要么通过图像来评价，本问题中选择了图像，更加清

23、晰直观。III、数值积分1、 实验目的熟悉并掌握数值积分的方法，重要训练复化梯形公式，复化simpson公式以及romberg积分。2、 问题描述问题三数值积分椭圆周长的计算。考虑椭圆,为计算其周长,只要计算其第一象限的长度即可.用参数方程可以表示为,计算公式为为计算方便,我们可以令,即计算下面的积分 (可以归结为上面的形式)采用复化梯形公式，复化Simpson公式以及Romberg积分的方法计算积分给出通用程序，该通用程序可以计算任何一个函数在任意一个区间在给定的精度下的数值积分。程序输出为计算出的数值积分值以及计算函数值的次数。利用你的程序计算5个椭圆的周长。3、 算法介绍首先利用给出的各

24、迭代公式，设计程序。在matlab对话框中输入要计算的函数，给出区间和精度。复化梯形的迭代公式为：复化simpson迭代公式为：Romberg迭代公式为：4、 算法程序复化梯形公式和复化simpson公式我们放在jifen.m中。（%标记后的程序可用来把b看为变量时的算法实现）%复化梯形公式function y=jifen(f,n,a,b) （说明：f表示任一函数，n精度，a，b为区间）fi=f(a)+f(b);h=(b-a)/n;d=1;%function f=jifen(n,a,b,c)%syms t%y=sqrt(1+(c2-1)*cos(t)2);%ya=subs(y,t,a);%yb

25、=subs(y,t,b);%fi=ya+yb;for i=1:n-1 x=a+i*h; fi=fi+2*f(x); d=d+1; %yx=subs(y,t,x); %fi=fi+2*yx;endf4=h/2*fi,d%复化simposon公式f1=0;f2=0;dd=1;for i=1:n-1 dd=dd+1; if rem(i,2)=0; x1=a+i*h; f1=f1+f(x1); else rem(i,2)=0; x2=a+i*h; f2=f2+f(x2) ; end endf3=(h/3)*(f(a)+4*f1+2*f2+f(b),dd romberg积分建立romberg.m文件fu

26、nction y=romberg(f,n,a,b) （说明：f表示任一函数，n精度，a，b为区间）z=zeros(n,n);h=b-a;z(1,1)=(h/2)*(f(a)+f(b);f1=0;for i=2:n for k=1:2(i-2) f1=f1+f(a+(k-0.5)*h); end z(i,1)=0.5*z(i-1,1)+0.5*h*f1; h=h/2; f1=0; for j=2:i z(i,j)=z(i,j-1)+(z(i,j-1)-z(i-1,j-1)/(4(j-1)-1); endendz,n5、 实验结果（1）对于复化梯形公式和复化simpson公式，我们运行下列语句并得

27、到结果：题中将椭圆中的未知量a取为1，b取为0.5。f4为复化梯形公式得到的1/4个椭圆周长，f3为复化simpson公式得到的1/4个椭圆周长。从而得到椭圆的周长为4*1.2111=4.8444 。（2） 对于romberg，运行下列语句并最终得到结果为：题中将椭圆中的未知量a取为1，b取为0.5。用romberg算法对椭圆的周长进行计算的到椭圆的周长为4.8442 。6、 结果分析我们计算了当椭圆长轴为1，短轴为0.5时的周长。通过上述三种方法的计算可以看到，结果相差不大。根据椭圆周长的一个计算公式我们可以得到L=4.283。因此三种方法都较好的接近真值。7、 心得体会我们应该熟练掌握这三

28、种方法，才能在编程时正确快速的写出迭代公式。同时在一种思想的前提下可以寻找多种方法实现算法，互相验证，从而得到更准确的结果。IV、newton差值和hermite差值方法1、 实验目的掌握并能够利用newton差值和hermite差值方法解决问题。2、 问题描述上述函数的导数为采用三种方法中最好的方法计算这一积分(1)利用数值积分的方法给出在(可以直接计算精确值的，用精确值)，用Newton插值方法得到5个椭圆的周长(2)利用数值积分的方法给出在(可以直接计算精确值的，用精确值)，用Hermite插值方法得到5个椭圆的周长(3) 选做题：利用以及导数更多的值来进行插值，插值误差会有什么变化？(

29、4)选做题：采用其它的插值方法改进插值的效果。3、 算法介绍a确定，对于给定的b值都对应着一个椭圆，在本问题中用newton插值法和hermite得到的多项式代替椭圆周长公式中的进行积分，首先画出图像，选择初始点。图像的实现代码见picture1.m。newton插值法迭代公式：Hermite法迭代公式：作图时建立picture.m文件画出和其导数图像。（注：此图像为b=0.5时）根据图像，我们选取我们选取0,0.3,0.6,0.9,1.2,1.5为初始点。4、 算法程序（1）建立newtondedai1.m文件。function z=newtondedai1(f,n)syms xia=zer

30、os(n,n);x=0 0.3 0.6 0.9 1.2 1.5;y=feval(f,x);a(:,1)=y;for i=2:n for j=2:i a(i,j)=(a(i,j-1)-a(i-1,j-1)/(x(1,i)-x(1,i-j+1); endendt=xi-x(1,1);p=a(1,1);for i=2:n p=p+a(i,i)*t; t=t*(xi-x(1,i);endp=collect(vpa(p)（2）建立hermite3.m文件。function z=hermite3(fn,dfn,n)syms xia=zeros(2*n,2*n);x=0.3 0.6 0.9; xx=0.3

31、0.3 0.6 0.6 0.9 0.9;y=feval(fn,x);yy=feval(dfn,x); for i=1:n z(1,2*i)=x(1,i); z(1,2*i-1)=x(1,i); a(2*i,1)=y(1,i); a(2*i-1,1)=y(1,i); a(2*i,2)=yy(1,i); end for i=1:n if 2*i6 a(2*i+1,2)=(a(2*i+1,1)-a(2*i,1)/(z(1,2*i+1)-z(1,2*i); end end for i=3:(2*n) for j=3:i a(i,j)=(a(i,j-1)-a(i-1,j-1)/(z(1,i)-z(1,i

32、-j+1); end end p=a(1,1); t=xi-xx(1,1); for i=2:2*n p=a(i,i)*t+p; t=t*(xi-xx(1,i); end p=vpa(collect(p) %p=vpa(collect(p),8)5、 实验结果问题（1）newton插值法，首先得到当a=1，b=0.5时用newton插值得到的多项式逼近椭圆周长表达式。执行下列语句其中p所表示的多项式就是newton插值方法得到的多项式。若保留3位小数，既接下来将得到的p多项式，调用复化梯形公式和复化simpson公式。说明即对newton插值法在（0,2pi）上积分得到椭圆周长为4.8440和4.8436 。问题（2）hermite插值法首先得到当a=1，b=0.5时用hermite插值得到的多项式逼近椭圆周长表达式。其中p所表示的多项式就是hermite插值方法得到的多项式。接下来将得到的p多项式，调用复化梯形公式和复化simpson公式。即对hermite插值法在（0,2pi）上调用复化梯形公式和复化simpson公式积分得到椭圆周长为