文档介绍：该【《数值分析》上机实验报告 】是由【小果冻】上传分享，文档一共【21】页，该文档可以免费在线阅读，需要了解更多关于【《数值分析》上机实验报告 】的内容，可以使用淘豆网的站内搜索功能，选择自己适合的文档，以下文字是截取该文章内的部分文字，如需要获得完整电子版，请下载此文档到您的设备，方便您编辑和打印。数值分析上机实验报告1?数值分析?-X4+14=0在〔,〕中的近似根〔初始近似值取为区间端点,〕。1理论依据:设函数在有限区间[a,b]上二阶导数存在,,逼近x的真实根。当前后两个的差<=ε时,就认为求出了近似的根。本程序用Newton法求代数方程(最高次数不大于10)在〔a,b〕区间的根。:#include<>#include<>main(){doublex2,f,f1;doublex1=;//{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)>=||x1<);//限制循环次数printf("计算结果:x=%f\n",x1);}:=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;End>>x0=;>>eps=;>>M=100;>>x=Newton('f','df',x0,eps,M);>>vpa(x,7):,用误差来控制循环迭代次数,可以在误差允许的范围内得到比较理想的计算结果。此程序的缺乏之处是,所要求解的方程必须满足上述定理的四个条件,但是第二和第四个条件在计算机上比较难以实现。,他的几何意义Xi+1是Xi的切线与x轴的交点,故也称为切线法。它是平方收敛的,但它是局部收敛的,即要求初始值与方程的根充分接近,所以在计算过程中需要先确定初始值。,讨论了区间(,)两端点是否能作为Newton迭代的初值,,,能作为初值。另外,该程序简单,只有一个循环,且为顺序结构,故采用do-while循环。当然也可以选择for和while循环。:x12345f(x)(x)’(x)f’(1)=1f’(10)=()及f’()的近似值。,所以只要求出,就能得出插值函数S〔x〕。求的方法为:这里最终归结为求解一个三对角阵的解。5用追赶法解三对角阵的方法如下:,综上可得求解方程Ax=d的算法::#include<>#include<>voidmain(){inti,j,m,n,k,p;doubleq10,p10,s4,g4,x0,x1,g0=1,g9=;;doubles[10][10];doublea[10],b[10],c[10],d[10],e[10],x[10],h[9],u[9],r[9];doublef[10]={0,,,,,,,,,};printf("请依次输入xi:\n");for(i=0;i<=9;i++)6scanf("%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("%",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=:\n");7x0=5-;x1=-4;s4=x[3]*pow(x0,3)/6+x[4]*pow(x1,3)/6+(f[3]-x[3]/6)*(5-)+(f[4]-x[4]/6)*(-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()的函数值是:%lf\nf()的导数值是:%lf\n",s4,g4);}:,没有Runge现象,光滑性较好。。,需要仔细弄清楚各数组所代表的参数,要注意各下标代表的含义,特别是在用追赶法计算的过程中。:Romberg算法的计算步骤如下:〔1〕先求出按梯形公式所得的积分值〔2〕把区间2等分,求出两个小梯形面积之和,记为,即这样由外推法可得,。〔3〕把区间再等分〔即等分〕,得复化梯形公式,由与外推可得,,如此,假设已算出等分的复化梯形公式,那么由Richardson外推法,构造新序列,m=1,2,…,l,k=1,2,…,l-m+1,最后求得。〔4〕或<就停止计算,否那么回到〔3〕,计算,一般可用如下算法:9其具体流程如下,并全部存入第一列通常计算时,用固定l=N来计算,一般l=4或5即能到达要求。:#include<>#include<>doublef(doublex)//计算f(x)的值{doublez;z=pow(3,x)*pow(x,)*(5*x+7)*sin(x*x);return(z);}main(){doublet[20][20],s,e=,a=1,b=3;inti,j,l,k;t[0][1]=(b-a)*(f(b)+f(a))/2;//下为romberg算法t[1][1]=(b-a)*(f(b)+2*f((b+a)/2)+f(a))/4;t[0][2]=(a*t[1][1]-t[0][1])/(4-1);j=3;for(l=2;fabs(t[0][j-1]-t[0][j-2])>=e;l++){for(k=1,s=0;k<=pow(2,l-1);k++)s+=f(a+(2*k-1)*(b-a)/pow(2,l));//判断前后两次所得的T(0)的差是否符合要求,如果符合精度要求那么停止循环t[l][1]=(t[l-1][1]+(b-a)*s/pow(2,l-1))/2;for(i=l-1,j=2;i>=0;i--,j++)t[i][j]=(pow(4,j-1)*t[i+1][j-1]-t[i][j-1])/(pow(4,j-1)-1);}if(t[0][1]<e)printf("t=%\n",t[0][1]);elseprintf("用Romberg算法计算函数所得近似结果为:\nf(x)=%\n",t[0][j-1]);}: