数值分析上机实验报告

数值分析上机实验报告《数值分析》上机实验报告1.用Newton 法求方程 X 7-X 4+14=0在（0.1，1.9）中的近似根（初始近似值取为区间端点，迭代6次或误差小于0.00001）。

1.1 理论依据：设函数在有限区间[a ，b]上二阶导数存在，且满足条件{}αϕ上的惟一解在区间平方收敛于方程所生的迭代序列迭代过程由则对任意初始近似值达到的一个中使是其中上不变号在区间],[0)(3,2,1,0,)(')()(],,[x |))(),((|,|,)(||)(|.4;0)(.3],[)(.20)()(.110......b a x f x k x f x f x x x Newton b a b f a f mir b a c x f ab c f x f b a x f b f x f k k k k k k ==-==∈≤-≠>+令)9.1()9.1(0)8(4233642)(0)16(71127)(0)9.1(,0)1.0(,1428)(3225333647>⋅''<-=-=''<-=-='<>+-=f f x x x x x f x x x x x f f f x x x f故以1.9为起点⎪⎩⎪⎨⎧='-=+9.1)()(01x x f x f x x k k k k 如此一次一次的迭代，逼近x 的真实根。

当前后两个的差<=ε时，就认为求出了近似的根。

本程序用Newton 法求代数方程(最高次数不大于10)在（a,b ）区间的根。

1.2 C语言程序原代码：#include<stdio.h>#include<math.h>main(){double x2,f,f1;double x1=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||x1<0.1); //限制循环次数printf("计算结果：x=%f\n",x1);}1.3 运行结果：1.4 MATLAB上机程序function y=Newton(f,df,x0,eps,M)d=0;for k=1:Mif feval(df,x0)==0d=2;breakelsex1=x0-feval(f,x0)/feval(df,x0);ende=abs(x1-x0);x0=x1;if e<=eps&&abs(feval(f,x1))<=epsd=1;breakendendif d==1y=x1;elseif d==0y='迭代M次失败';elsey= '奇异'endfunction y=df(x)y=7*x^6-28*4*x^3;Endfunction y=f(x)y=x^7-28*x^4+14;End>> x0=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循环。

2.已知函数值如下表：试用三次样条插值求f(4.563)及f ’(4.563)的近似值。

2.1 理论依据332211111111111()()()()()()()6666j j j j j j j jj j j j j j j j x x x x h x x h x x S x M M y M y M h h h h ---------------=++-+-这里11j j j h x x --=- ，所以只要求出j M ，就能得出插值函数S （x ）。

求j M 的方法为：00111122112122212N N N N M d M d M d μλμλμλ--⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦这里10000011111111116()6()(1,2,,1)61[()]1j j j j jj j j j N N N N N N j jj j j j j j j y y d y h hy y y y d j N h h h h d y y y h h h h h h h h μλμ+----------⎧'=-⎪⎪⎪--=-=-⎪+⎪⎨⎪'=--⎪⎪⎪==-=⎪++⎩最终归结为求解一个三对角阵的解。

用追赶法解三对角阵的方法如下：11112221222111111111n n n n n nnn n b c a b c l A LU l b c a l a b γβγββγβ-----⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥===⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦1,,,n L d LUx d L d Ux δδδδδδ⎡⎤=⎧⎢⎥===⎨⎢⎥=⎩⎢⎥⎣⎦即若记则由得 112111nn n d l l d δδ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦ ， 111111n n n n n x x βγδβγβδ--⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦综上可得求解方程Ax=d 的算法：1111111111111,,,,1,2,3,,1,,1,,2,1i i i i i i ii i i in i i i n i n i b d l b l c d l i n c xx x i n αβδββδδδδββ+++++++++⎧====-⎪⎪⎪=-=-⎨⎪-⎪===-⎪⎩2.2 C 语言程序代码：#include<stdio.h>#include<math.h>void main() {int i,j,m,n,k,p;double q10,p10,s4,g4,x0,x1,g0=1,g9=0.1;; double s[10][10];double a[10],b[10],c[10],d[10],e[10],x[10],h[9],u[9],r[9];double f[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];d[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;i<9;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);g4=-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现象，光滑性较好。