《数值分析》上机实验报告

数值分析上机实验报告1《数值分析》上机实验报告1.用Newton法求方程X7-X4+14=0在（0.1，1.9）中的近似根（初始近似值取为区间端点，迭代6次或误差小于0.00001）。1.1理论依据：设函数在有限区间[a，b]上二阶导数存在，且满足条件上的惟一解在区间平方收敛于方程所生的迭代序列迭代过程由则对任意初始近似值达到的一个中使是其中上不变号在区间],[0)(3,2,1,0,)(')()(],,[x|))(),((|,|,)(||)(|.4;0)(.3],[)(.20)()(.110......baxfxkxfxfxxxNewtonbabfafmirbacxfabcfxfbaxfbfxfkkkkkk令0)9.1()9.1(0)8(4233642)(0)16(71127)(0)9.1(,0)1.0(,1428)(3225333647ffxxxxxfxxxxxfffxxxf故以1.9为起点9.1)()(01xxfxfxxkkkk如此一次一次的迭代，逼近x的真实根。当前后两个的差=ε时，就认为求出了近似的根。本程序用Newton法求代数方程(最高次数不大于10)在（a,b）区间的根。21.2C语言程序原代码：#includestdio.h#includemath.hmain(){doublex2,f,f1;doublex1=1.9;//取初值为1.9do{x2=x1;f=pow(x2,7)-28*pow(x2,4)+14;f1=7*pow(x2,6)-4*28*pow(x2,3);x1=x2-f/f1;}while(fabs(x1-x2)=0.00001||x10.1);//限制循环次数printf(计算结果：x=%f\n,x1);}1.3运行结果：1.4MATLAB上机程序functiony=Newton(f,df,x0,eps,M)d=0;fork=1:Miffeval(df,x0)==0d=2;breakelsex1=x0-feval(f,x0)/feval(df,x0);ende=abs(x1-x0);x0=x1;ife=eps&&abs(feval(f,x1))=epsd=1;breakendend3ifd==1y=x1;elseifd==0y='迭代M次失败';elsey='奇异'endfunctiony=df(x)y=7*x^6-28*4*x^3;Endfunctiony=f(x)y=x^7-28*x^4+14;Endx0=1.9;eps=0.00001;M=100;x=Newton('f','df',x0,eps,M);vpa(x,7)1.5问题讨论：1.使用此方法求方解，用误差来控制循环迭代次数，可以在误差允许的范围内得到比较理想的计算结果。此程序的不足之处是，所要求解的方程必须满足上述定理的四个条件，但是第二和第四个条件在计算机上比较难以实现。2.Newton迭代法是一个二阶收敛迭代式，他的几何意义Xi+1是Xi的切线与x轴的交点,故也称为切线法。它是平方收敛的，但它是局部收敛的，即要求初始值与方程的根充分接近，所以在计算过程中需要先确定初始值。3.本题在理论依据部分，讨论了区间(0.1,1.9)两端点是否能作为Newton迭代的初值，结果发现0.1不满足条件，而1.9满足，能作为初值。另外，该程序简单，只有一个循环，且为顺序结构，故采用do-while循环。当然也可以选择for和while循环。42.已知函数值如下表：x12345f(x)00.693147181.09861231.38629441.6094378x678910f(x)1.79175951.94591012.0794452.19722462.3025851f’(x)f’(1)=1f’(10)=0.1试用三次样条插值求f(4.563)及f’(4.563)的近似值。2.1理论依据332211111111111()()()()()()()6666jjjjjjjjjjjjjjjjxxxxhxxhxxSxMMyMyMhhhh这里11jjjhxx，所以只要求出jM，就能得出插值函数S（x）。求jM的方法为：00111122112122212NNNNMdMdMd这里10000011111111116()6()(1,2,,1)61[()]1jjjjjjjjjNNNNNNjjjjjjjjjyydyhhyyyydjNhhhhdyyyhhhhhhhh最终归结为求解一个三对角阵的解。5用追赶法解三对角阵的方法如下：11112221222111111111nnnnnnnnnbcabclALUlbcalab1,,,nLdLUxdLdUx即若记则由得112111nnndlld，111111nnnnnxx综上可得求解方程Ax=d的算法：1111111111111,,,,1,2,3,,1,,1,,2,1iiiiiiiiiiiniiininibdlblcdlincxxxin2.2C语言程序代码：#includestdio.h#includemath.hvoidmain(){inti,j,m,n,k,p;doubleq10,p10,s4,g4,x0,x1,g0=1,g9=0.1;;doubles[10][10];doublea[10],b[10],c[10],d[10],e[10],x[10],h[9],u[9],r[9];doublef[10]={0,0.69314718,1.0986123,1.3862944,1.6094378,1.7917595,1.9459101,2.079445,2.1972246,2.3025851};printf(请依次输入xi:\n);for(i=0;i=9;i++)scanf(%lf,&e[i]);//求h矩阵for(n=0;n=8;n++)h[n]=e[n+1]-e[n];d[0]=6*((f[1]-f[0])/h[0]-g0)/h[0];6d[9]=6*(g9-(f[9]-f[8])/h[8])/h[8];for(j=0;j=7;j++)d[j+1]=6*((f[j+2]-f[j+1])/h[j+1]-(f[j+1]-f[j])/h[j])/(h[j]+h[j+1]);for(m=1;m=8;m++)u[m]=h[m-1]/(h[m-1]+h[m]);for(k=1;k=8;k++)r[k]=h[k]/(h[k-1]+h[k]);for(i=0;i=9;i++)//求u矩阵for(p=0;p=9;p++){s[i][p]=0;if(i==p)s[i][p]=2;}s[0][1]=1;s[9][8]=1;for(i=1;i=8;i++){s[i][i-1]=u[i];s[i][i+1]=r[i];}printf(三对角矩阵为:\n);for(i=0;i=9;i++)for(p=0;p=9;p++)//求r矩阵{printf(%5.2lf,s[i][p]);if(p==9){printf(\n);}}printf(根据追赶法解三对角矩阵得：\n);a[0]=s[0][0];b[0]=d[0];for(i=1;i9;i++){c[i]=s[i][i-1]/a[i-1];//求d矩阵a[i]=s[i][i]-s[i-1][i]*c[i];b[i]=d[i]-c[i]*b[i-1];if(i==8){p10=b[i];q10=a[i];}}x[9]=p10/q10;printf(M[10]=%lf\n,x[9]);for(i=9;i=1;i--){x[i-1]=(b[i-1]-s[i-1][i]*x[i])/a[i-1];printf(M[%d]=%lf\n,i,x[i-1]);}printf(可得s(x)在区间[4,5]上的表达式;\n);printf(将x=4.563代入得：\n);x0=5-4.563;x1=4.563-4;s4=x[3]*pow(x0,3)/6+x[4]*pow(x1,3)/6+(f[3]-x[3]/6)*(5-4.563)+(f[4]-x[4]/6)*(4.563-4);7g4=-x[3]*pow(x0,2)/2+x[4]*pow(x1,2)/2-(f[3]-x[3]/6)+(f[4]-x[4]/6);printf(计算结果：f(4.563)的函数值是：%lf\nf(4.563)的导数值是：%lf\n,s4,g4);}2.3运行结果：2.4问题讨论1.三次样条插值效果比Lagrange插值好，没有Runge现象，光滑性较好。2.本题的对任意划分的三弯矩插值法可以解决非等距节点的一般性问题。3.编程过程中由于定义的数组比较多，需要仔细弄清楚各数组所代表的参数，要注意各下标代表的含义，特别是在用追赶法计算的过程中。83.用Romberg算法求)00001.0(sin)75(32314.1允许误差dxxxxx.3.1理论依据：Romberg算法的计算步骤如下：（1）先求出按梯形公式所得的积分值(0)1[()()]2baTfafb（2）把区间2等分，求出两个小梯形面积之和，记为(1)1T，即(1)1[()()2()]42babaTfafbf这样由外推法可得(0)2T,(1)2(0)(1)(0)11(0)11221()421411()2TTTTT。（3）把区间再等分（即22等分），得复化梯形公式(2)1T，由(1)1T与(2)1T外推可得(2)(1)(1)112441TTT，2(1)(0)(0)2232441TTT，如此，若已算出2k等分的复化梯形公式()1kT，则由Richardson外推法，构造新序列()(1)(1)1441mkkkmmmmTTT，m=1,2,…,l,k=1,2,…,l-m+1,最后求得(0)1lT。（4）(0)(0)1llTT或(0)(0)1||llTT就停止计算，否则回到（3），计算(1)1lT，一般可用如下算法：1(0)12()(1)1111()(1)(1)1[()()]21{[(21)]}2224,1,2,,,1,2,,141lllllimkkkmmmmbaTfafbbabaTTfaiTTTmlklm9其具体流程如下，并全部存入第一列(0)1(1)T(0)2(3)T(0)3(6)T(1)1(2)T(1)2(5)T(2)1(4)T通常计算时，用固定l=N来计算，一般l=4或5即能达到要求。3.2C语言程序代码：#includemath.h#includestdio.hdoublef(doublex)//计算f(x)的值{doublez;z=pow(3,x)*pow(x,1.4)*(5*x+7)*sin(x*x);return(z);}main(){doublet[20][20],s,e=0.00001,a=1,b=3;inti,j,l,k;t[0][1]=(b-a)*(f(