哈工大数值分析上机实验报告Word文档格式.docx

1、 end R=b-a;%求出误差k=k+1;endx=c%给出解Newton法及改良的Newton法源程序：% 输入函数f=input（请输入需要求解函数,s）%求解f（x）的导数df=diff（f）;%改良常数或重根数miu=2;%初始值x0x0=input（input initial value x0）;%迭代次数max=100;%最大迭代次数R=eval（subs（f,x0x）;%求解f（x0），以确定初值x0时否就是解while （abs（R）1e-8） x1=x0-miu*eval（subs（f,）/eval（subs（df, R=x1-x0; x0=x1; k=k+1;if （ev

2、al（subs（f,）max;%如果迭代次数大于给定值，认为迭代不收敛，重新输入初值 ss=input（maybe result is error,choose a new x0,y/n? if strcmp（ss,y x0=input（ k=0;k;%给出迭代次数x=x0；%给出解结果分析和讨论：1.用二分法计算方程在1，2的根。（,下同）计算结果为x= 1.523；f（x）= -3.4311e-007；k=18；由f（x）知结果满足要求，但迭代次数比拟多，方法收敛速度比拟慢。2.用二分法计算方程在1，1.5的根。x= 1.180；f（x）= 2.4815e-006；k=17；由f（x）知结

3、果满足要求，但迭代次数还是比拟多。3.用Newton法求解以下方程a） x0=0.5；x= 0.978；f（x）= 2.0313e-016；k=4；由f（x）知结果满足要求，而且又迭代次数只有4次看出收敛速度很快。b） x0=1；c） x0=0.45, x0=0.65； 当x0=0.45时，计算结果为x= 0.983；f（x）= -8.4584e-014；由f（x）知结果满足要求，而且又迭代次数只有4次看出收敛速度很快，实际上该方程确实有真解x=0.5。当x0=0.65时，计算结果为x= 0.000；f（x）=0；k=9；由f（x）知结果满足要求，实际上该方程确实有真解x=0.5，但迭代次数增

4、多，实际上当取x00.68时，x1，就变成了方程的另一个解，这说明Newton法收敛与初值很有关系，有的时候甚至可能不收敛。4.用改良的Newton法求解，有2重根，取 x0=0.55；并与3.中的c）比拟结果。当x0=0.55时，程序死循环，无法计算，也就是说不收敛。改时，结果收敛为x=0.286；f（x）=4.1127e-007；k=16；显然这个结果不是很好，而且也不是收敛至方程的2重根上。当x0=0.85时，结果收敛为x= 1.489；f（x）= 2.8737e-023；这次到达了预期的结果，这说明初值的选取很重要，直接关系到方法的收敛性，实际上直接用Newton法，在给定同样的条件和

5、精度要求下，可得其迭代次数k=15，这说明改良后的Newton法法速度确实比拟快。结论： 对于二分法，只要能够保证在给定的区间有根，使能够收敛的，当时收敛的速度和给定的区间有关，二且总体上来说速度比拟慢。Newton法，收敛速度要比二分法快，但是最终其收敛的结果与初值的选取有关，初值不同，收敛的结果也可能不一样，也就是结果可能不时预期需要得结果。改良的Newton法求解重根问题时，如果初值不当，可能会不收敛，这一点非常重要，当然初值适宜，一样情况下其速度要比Newton法快得多。实验报告二 Gauss列主元消去法求解线性方程组的方法很多，主要分为直接法和间接法。本实验运用直接法的Guass消去

6、法,并采用选主元的方法对方程组进展求解。目的和意义1.学习Gauss消去法的原理。2.了解列主元的意义。3.确定什么时候系数阵要选主元由于一般线性方程在使用Gauss消去法求解时，从求解的过程中可以看到，假设=0，那么必须进展行交换，才能使消去过程进展下去。有的时候即使0，但是其绝对值非常小，由于机器舍入误差的影响，消去过程也会出现不稳定得现象，导致结果不正确。因此有必要进展列主元技术，以最大可能的消除这种现象。这一技术要寻找行r，使得并将第r行和第k行的元素进展交换，以使得当前的的数值比0要大的多。这种列主元的消去法的主要步骤如下：1.消元过程对k=1,2,n-1,进展如下步骤。1）选主元，

7、记假设很小，这说明方程的系数矩阵严重病态，给出警告，提示结果可能不对。2）交换增广阵A的r，k两行的元素。 （j=k,n+1）3）计算消元 （i=k+1,n; j=k+1,n+1）2.回代过程对k= n, n-1,1,进展如下计算至此，完成了整个方程组的求解。Gauss消去法源程序：a=input（输入系数阵：nb=input（输入列阵b：n=length（b）;A=a bx=zeros（n,1）;%函数主体for k=1:n-1;%是否进展主元选取if abs（A（k,k）abs（t） p=r; else p=k; %交换元素 if p=k; for q=k:n+1; s=A（k,q）; A

8、（k,q）=A（p,q）; A（p,q）=s; end %判断系数矩阵是否奇异或病态非常严重 yipusilongdisp（矩阵奇异，解可能不正确） %计算消元,得三角阵 m=A（r,k）/A（k,k）; A（r,q）=A（r,q）-A（k,q）*m; %求解x x（n）=A（n,n+1）/A（n,n）; for k=n-1:-1:1; s=0; s=s+A（k,r）*x（r）; t=（A（k,n+1）-s） x（k）=（A（k,n+1）-s）/A（k,k）例：求解方程。其中为一小数，当时，分别采用列主元和不列主元的Gauss消去法求解，并比拟结果。记Emax为求出的解代入方程后的最大误差，按

9、要求，计算结果如下：当时，不选主元和选主元的计算结果如下，其中前一列为不选主元结果，后一列为选主元结果，下同。 0.391 0.651 2.972 2.163 2.451 2.721Emax= 9.2624e-010，0此时，由于不是很小，机器误差就不是很大，由Emax可以看出不选主元的计算结果精度还可以，因此此时可以考虑不选主元以减少计算量。时，不选主元和选主元的计算结果如下 1.877 0.348 1.807 2.174 3.731 2.609Emax= 2.4668e-005，0此时由Emax可以看出不选主元的计算精度就不好了，误差开场增大。 1.020 1.000 1.666 2.00

10、0 3.111 0000Emax= 0.503，0此时由Emax可以看出，不选主元的结果应该可以说是不正确了，这是由机器误差引起的。NaN 1NaN 2NaN 3Emax=NaN, 0不选主元时，程序报错：Warning: Divide by zero.。这是因为机器计算的最小精度为10-15，所以此时的就认为是0，故出现了错误现象。而选主元时那么没有这种现象，而且由Emax可以看出选主元时的结果应该是准确解。采用Gauss消去法时，如果在消元时对角线上的元素始终较大假设大于10-5，那么本方法不需要进展列主元计算，计算结果一般就可以到达要求，否那么必须进展列主元这一步，以减少机器误差带来的影

11、响，使方法得出的结果正确。实验报告三 Rung现象产生和克制由于高次多项式插值不收敛，会产生Runge现象，本实验在给出具体的实例后，采用分段线性插值和三次样条插值的方法有效的克制了这一现象，而且还取的很好的插值效果。1.深刻认识多项式插值的缺点。2.明确插值的不收敛性怎样克制。3.明确精度与节点和插值方法的关系。在给定n+1个节点和相应的函数值以后构造n次的Lagrange插值多项式，实验结果说明见后面的图这种多项式并不是随着次数的升高对函数的逼近越来越好，这种现象就是Rung现象。解决Rung现象的方法通常有分段线性插值、三次样条插值等方法。分段线性插值：设在区间a, b上，给定n+1个插

12、值节点a=x0x1xn=b和相应的函数值y0，y1，yn，求作一个插值函数，具有如下性质：1），j=0，1，n。2）在每个区间xi, xj上是线性连续函数。那么插值函数称为区间a, b上对应n个数据点的分段线性插值函数。三次样条插值：给定区间a, b一个分划：xN=b假设函数S（x）满足以下条件：1）S（x）在每个区间xi, xj上是不高于3次的多项式。2）S（x）及其2阶导数在a, b上连续。那么称S（x）使关于分划的三次样条函数。其中待插值的方程写成function的方式，如下y=1/1+25*x*x;写成如上形式即可，下面给出主程序 Lagrange插值源程序：n=input（将区间分为

13、的等份数输入：s=-1+2/n*0:n;%给定的定点，Rf为给定的函数x=-1:0.01:f=0;for q=1: l=1;%求插值基函数 for k=1: if k=q; l=l.*（x-s（k）./（s（q）-s（k）; l=l; f=f+Rf（s（q）*l;%求插值函数plot（x,f,r）%作出插值函数曲线grid on hold on分段线性插值源程序m=0;hh=0.001;for x=-1:hh: ff=0;%求插值基函数 switch k case 1 if xs（n）; l=（x-s（n）./（s（n+1）-s（n）; otherwise =s（k-1）&x=s（k）&=s（

14、k+1）; l=（x-s（k+1）./（s（k）-s（k+1）; else ff=ff+Rf（s（k）*l;%求插值函数值 m=m+1; f（m）=ff;end %作出曲线grid onhold on 三次样条插值源程序：采用第一边界条件%插值区间a=-1;b=1;%画图的步长s=a+（b-a）/n*0:%第一边界条件Rf（-1）,Rf（1）v=5000*1/（1+25*a*a）3-50/（1+25*a*a）4;%取出节点间距 h（k）=s（k+1）-s（k）;%求出系数向量lamuda,miula（k）=h（k+1）/（h（k+1）+h（k）; miu（k）=1-la（k）;%赋值系数矩阵A

15、 for p=1: switch p case k A（k,p）=2; case k-1 A（k,p）=miu（p+1）; case k+1 A（k,p）=la（p-1）; otherwise A（k,p）=0;%求出d阵 d（k）=6*f2c（s（k） s（k+1） s（k+2）-miu（k）*v; case n-1 d（k）=6*f2c（s（k） s（k+1） s（k+2）-la（k）*v; d（k）=6*f2c（s（k） s（k+1） s（k+2）;%求解M阵M=Ad;M=v;M;v;%for x=a:b; if x=a; p=1; p=ceil（x-s（1）/（b-a）/n）; ff1

16、=0; ff2=0; ff3=0; ff4=0; ff1=1/h（p）*（s（p+1）-x）3*M（p）/6; ff2=1/h（p）*（x-s（p）3*M（p+1）/6; ff3=（Rf（s（p+1）-Rf（s（p）/h（p）-h（p）*（M（p+1）-M（p）/6）*（x-s（p）; ff4=Rf（s（p）-M（p）*h（p）*h（p）/6; f（m）=ff1+ff2+ff3+ff4 ;%作出插值图形x=a:k本实验采用函数进展数值插值，插值区间为-1，1，给定节点为xj=-1+jh,h=0.1，j=0,，n。下面分别给出Lagrange插值，三次样条插值，线性插值的函数曲线和数据表。图中只

17、标出Lagrange插值的十次多项式的曲线，其它曲线没有标出，从数据表中可以看出具体的误差。表中，L10（x）为Lagrange插值的10次多项式，S10（x），S40（x）分别代表n=10，40的三次样条插值函数，X10（x），X40（x）分别代表n=10，40的线性分段插值函数。xf（x）L10（x）S10（x） S40（x） X10（x） X40（x） -1.000 0.154 0.154 0.154 0.154 0.154 0.154 -0.000 0.239 1.920 0.040 0.239 0.910 0.239 -0.000 0.941 1.926 0.458 0.941 0.

18、665 0.941 -0.000 0.344 0.982 0.979 0.344 0.421 0.344 -0.000 0.176 0.176 0.176 0.176 0.176 0.176 -0.000 0.378 -0.674 0.744 0.378 0.882 0.378 -0.000 0.321 -0.250 0.866 0.321 0.588 0.321 -0.000 0.649 -0.418 0.849 0.649 0.294 0.649 -0.000 0.000 0.000 0.000 0.000 0.000 0.000 -0.000 0.788 0.257 0.713 0.78

19、8 0.000 0.788 -0.000 0.276 0.103 0.730 0.276 0.000 0.276 -0.000 0.825 0.267 0.883 0.825 0.000 0.825 -0.000 0.385 0.376 0.464 0.385 0.000 0.385 -0.000 0.231 0.080 0.860 0.231 0.000 0.231 -0.000 0.902 0.789 0.327 0.902 0.000 0.902 -0.000 0.000 0.340 0.431 0.000 0.000 0.000 -0.000 0.000 0.890 0.828 0.0

20、00 0.000 0.000 -0.000 0.824 0.073 0.810 0.824 0.000 0.824 0 1.000 1.000 1.000 1.000 1.000 1.000 0.000 0.824 0.073 0.810 0.824 0.000 0.824 0.000 0.000 0.890 0.828 0.000 0.000 0.000 0.000 0.000 0.340 0.431 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.902 0.789 0.327 0.902 0.000

21、0.902 0.000 0.231 0.080 0.860 0.231 0.000 0.231 0.000 0.385 0.376 0.464 0.385 0.000 0.385 0.000 0.825 0.267 0.883 0.825 0.000 0.825 0.000 0.276 0.103 0.730 0.276 0.000 0.276 0.000 0.788 0.257 0.713 0.788 0.000 0.788 0.000 0.649 -0.418 0.849 0.649 0.294 0.649 0.000 0.321 -0.250 0.866 0.321 0.588 0.32

22、1 0.000 0.378 -0.674 0.744 0.378 0.882 0.378 0.000 0.176 0.176 0.176 0.176 0.176 0.176 0.000 0.344 0.982 0.979 0.344 0.421 0.344 0.000 0.941 1.926 0.458 0.941 0.665 0.941 0.000 0.239 1.920 0.040 0.239 0.910 0.2391.000 0.154 0.154 0.154 0.154 0.154 0.154 从以上结果可以看到，用三次样条插值和线性分段插值，不会出现多项式插值是出现的Runge现象，

23、插值效果明显提高。进一步说，为了提高插值精度，用三次样条插值和线性分段插值是可以增加插值节点的方法来满足要求，而用多项式插值函数时，节点数的增加必然会使多项式的次数增加，这样会引起数值不稳定，所以说这两种插值要比多项式插值好的多。而且在给定节点数的条件下，三次样条插值的精度要优于线性分段插值，曲线的光滑性也要好一些。实验报告四 多项式最小二乘法对于具体实验时，通常不是先给出函数的解析式，再进展实验，而是通过实验的观察和测量给出离散的一些点，再来求出具体的函数解析式。又因为测量误差的存在，实际真实的解析式曲线并不一定通过测量给出的所有点。最小二乘法是求解这一问题的很好的方法，本实验运用这一方法实现对给定数据的拟合。1.