数值分析-上机实验(精简版)

《数值分析》上机实习报告学院：专业：姓名：学号：指导教师：使用语言：C语言使用软件：VC++6.02013年12月20日第一题：列主元素法求方程组根1.题目用列主元消去法求解Ax=b。2.理论依据和应用条件:所谓列主元消去法是，对矩阵作恰当的调整，选取绝对值最大的元素作为主元素。然后把矩阵化为上三角阵，再进行回代，求出方程的解。此题采用宏定义将其一般化了，对于类似的方程都可以运用此法求解。3.计算程序：#include<stdio.h>#include<math.h>#defineN9/*宏定义N=9*/voidmain(){ intn,t,k,i,j; doublem,h,x[N],A[N][N+1]=/*A和b合并成A[N][N+1]*/ {{12.38412,2.115237,-1.061074,1.112336,-0.113584,0.718719,1.74238,3.067813,-2.031743,2.1874369}, {2.115237,19.141823,-3.125432,-1.012345,2.189736,1.563849,-0.784165,1.112348,3.123124,33.992318}, {-1.061074,-3.125432,15.567914,3.123848,2.031454,1.836742,-1.056781,0.336993,-1.010103,-25.173417}, {1.112336,-1.012345,3.123848,27.108437,4.101011,-3.741856,2.101023,-0.71828,-0.037585,0.84671695}, {-0.113584,2.189736,2.031454,4.101011,19.897918,0.431637,-3.111223,2.121314,1.784317,1.784317}, {0.71819,1.563849,1.836742,-3.741856,0.431637,9.789365,-0.103458,-1.10356,0.238417,-86.612343}, {1.742382,-0.784165,-1.056781,2.101023,-3.111223,-0.103458,14.7138465,3.123789,-2.213474,1.1101230}, {3.067813,1.112348,0.336993,-0.71828,2.121314,-1.103456,3.123789,30.719334,4.446782,4.719345}, {-2.031743,3.123124,-1.010103,-0.037585,1.784317,0.238417,-2.213474,4.446782,40.00001,-5.6784392}}; n=N;n=n-1;for(k=0;k<n;k++) {m=fabs(A[k][k]);t=k;for(i=k+1;i<=n;i++)if(fabs(A[i][k])>m)/*选取列主元m*/ { m=fabs(A[i][k]); t=i;/*记录行号t*/ }for(j=k;j<=n+1;j++)/*换行*/ { h=A[k][j]; A[k][j]=A[t][j]; A[t][j]=h; }for(i=k+1;i<=n;i++)/*消元并计算*/for(j=k+1;j<=n+1;j++)A[i][j]=A[i][j]-A[i][k]*A[k][j]/A[k][k]; } x[n]=A[n][n+1]/A[n][n];/*回带计算x*/for(k=n-1;k>=0;k--) { x[k]=A[k][n+1];for(j=k+1;j<=n;j++)x[k]=x[k]-A[k][j]*x[j];x[k]=x[k]/A[k][k]; }for(i=0;i<=n;i++)/*输出x*/printf("x[%d]=%f;\n",i+1,x[i]);}4.实验结果5.问题讨论〔1〕首先，我们将A和b的数值定义到一个数组中了，这方便的程序的调用；〔2〕其次，程序的实现，关键在于主元的选取和换行，需要仔细的思考；〔3〕再次，本人将程序一般化了，类似方程组，都可以通过此算法实现。第二题：三次样条插值求近似值题目函数值如下表：x12345f(x)00.6931471.09861231.38629441.6094378x678910f(x)1.79175951.94591012.0794452.19722462.3025851f'(x)f'(1)=1f'(10)=0.1试用三次样条插值求f(4.563)和f'(4.563)的近似值2.理论依据和应用条件:依据三弯矩插值法列出相应的方程，然后将x的值带入〔在求解三弯矩方程的参数时，应用追赶法求解方程组求出相应参数值〕。运用到课本上4.7.2公式以及S〔x〕的一阶导数的公式，4.7.4公式，4.7.5公式以及追赶法的相关公式。3.计算程序：#include<stdio.h>#defineN10/*计算10个节点*/intx[N]={1,2,3,4,5,6,7,8,9,10},i,h=1;/*定义全局变量x[N],i,h*/doublef[N]={0,0.69314718,1.0986123,1.3862944,1.6094378, 1.7917595,1.9459101,2.079445,2.1972246,2.3025851};voidmain(){doubleM(intk),s,s1;/*调用doubleM(intk)函数*/ floatx0,x1,y; intj;printf("X=");scanf("%f",&y);j=int(y-1);/*因为数组的第一个值是x[0],不是x[1]*/ x0=(float)x[j]; x1=(float)x[j+1];s=(x1-y)*(x1-y)*(x1-y)/(6*h)*M(j)+(y-x0)*(y-x0)*(y-x0)/(6*h)*M(j+1)/*求S(X)*/ +(f[j]-M(j)*h*h/6)*(x1-y)/h+(f[j+1]-M(j+1)*h*h/6)*(y-x0)/h;printf("S=%f\n",s);s1=-(x1-y)*(x1-y)*M(j)/(2*h)+(y-x0)*(y-x0)*M(j+1)/(2*h)/*求S'(X)即一阶导数*/ -(f[j]-M(j)*h*h/6)/h+(f[j+1]-M(j+1)*h*h/6)/h;printf("S1=%f\n",s1);}doubleM(intk)/*定义doubleM(intk)函数*/{ doubleu[N-1],v[N-1],b[N],l[N-1],r[N],z[N],n[N];/*v[N-1]即书上；b[N]即书上d；n[N]即书上x[N]*/doublefd0=1.0,fd9=0.1;floatd=2;for(i=0;i<=N-2;i++)/*追赶法求值*/if(i<N-2)/*求u[N-1]*/u[i]=(float)h/(h+h);elseu[i]=1;for(i=0;i<=N-2;i++)/*求v[N-1]*/if(i==0)v[i]=1;elsev[i]=1-u[i-1];for(i=0;i<=N-1;i++)/*求b[N]*/ if(i==0)b[i]=6*((f[i+1]-f[i])/h-fd0)/h;else {if(i==N-1)b[i]=6*(fd9-(f[i]-f[i-1])/h)/h;elseb[i]=6*((x[i]-x[i+1])*f[i-1]+(x[i+1]-x[i-1])*f[i]+(x[i-1] -x[i])*f[i+1])/((x[i-1]-x[i])*(x[i]-x[i+1])*(x[i-1]-x[i+1]));}for(i=0;i<=N-2;i++)/*求l[N-1]*/if(i==0)l[i]=u[i]/d;elsel[i]=u[i]/(d-l[i-1]*v[i-1]);for(i=0;i<N;i++)/*求r[N]*/if(i==0)r[i]=d;elser[i]=d-l[i-1]*v[i-1];for(i=0;i<N;i++)/*求z[N]*/if(i==0)z[i]=b[i];elsez[i]=b[i]-l[i-1]*z[i-1];for(i=N-1;0<=i;i--)/*求n[N]*/if(i==N-1)n[i]=z[i]/r[i];elsen[i]=(z[i]-v[i]*n[i+1])/r[i];return(n[k]);/*返回值即M(intk)*/}4.计算结果5.问题讨论在做数组相关的编程时，注意数组下标是从0开始的；编程的过程中，要慎用全局变量；算法的设计，尽量简洁；追赶法的使用过程中，注意对应项的转换和实现。第三题：牛顿迭代法求方程近似根题目用Newton法求方程在(0.1,1.9)中的近似根(初始近似值取为区间端点，迭6次或误差小于0.00001).2.理论依据和应用条件:根据Newton迭代法的2.3.1公式，本算法设计一般化了，类似方程，皆可以通过此算法实现。3.计算程序：#include<stdio.h>#include<math.h>voidmain(){floatx,x1,f,f1;x=1.9;do/*设置一个循环*/ {x1=x; f=(x1*x1*x1-28)*x1*x1*x1*x1+14; f1=(7*x1*x1*x1-112)*x1*x1*x1; x=x1-f/f1; }while(fabs(x-x1)>1e-6);/*判断是否满足精度*/ printf("Therootis：%.6f\n",x);}4.计算结果5.问题讨论〔1〕在运用Newton法时，确定初值的选取；〔2〕编程的过程要用到数学函数时，必须输入#include<math.h>。第四题：Romberg算法求积分1.题目用Romberg算法求.2.理论依据和应用条件:运用梯形公式5.1.2和Romberg算法的5.4.1公式、5.4.4公式、5.4.5公式、5.4.6公式求解。对于同类积分，皆可以通过这个算法求得。3.计算程序#include<stdio.h>#include<math.h>#defineN10/*先取最大值10，假设不行直接改N的值即可*/doubleb=3,a=1;voidmain(){ doublet[N],s[N-1],c[N-2],r[N-3];doublef(doublex1),sum;intk,i;for(i=0;i<=N-1;i++){if(i==0)/*开始求T[N]，从0开始考虑*/ t[i]=(f(a)+f(b))*(b-a)/2.0;else{for(k=0;k<=pow(2,i-1)-1;k++)if(k==0) sum=f(a+(b-a)/(pow(2,i-1))/2);elsesum=sum+f(a+(b-a)/(pow(2,i-1))/2+k*(b-a)/(pow(2,i-1)));/*前i项的求和*/t[i]=(t[i-1]+(b-a)/(pow(2,i-1))*sum)/2;}if(i>0) s[i-1]=(t[i]*4.0-t[i-1])/3.0;/*求s[N-1]*/if(i>1) c[i-2]=16*s[i-1]/15-s[i-2]/15;/*求c[N-2]*/if(i>2){r[i-3]=(64*c[i-2]-c[i-3])/63.0;/*求r[N-3]*/if((fabs(r[i-3]-r[i-4]))<(1e-5))break;/*误差小于1e-5时，结束*/}}printf("Theresultis%f\n",r[i-3]);/*输出结果*/}doublef(doublex1)/*定义函数f(x)*/{doublet2;t2=pow(3,x1)*pow(x1,1.4)*(5*x1+7)*(sin(x1*x1));return(t2);}4.计算结果5.问题讨论由于是顺序求值，不会产生冗余运算，先允许N为10，假设没结果；那么可以直接改变N的值即可；注意数组的下标是从0开始的，注意转换；不同数据类型的数据运算时，注意转换；算方法的设计过程中，注意算法的通用型，不要就题论题。第六题：定步长四阶Runge-Kutta法求方程组题目用定步长四阶Runge-Kutta求解h=0.0005，打印yi(0.025)，yi(0.045)，yi(0.085)，yi(0.1)，(i=1，2，3)2.理论依据和应用条件:运用四阶Runge-Kutta公式6.3.8，对于定步长的方程组，求解时，h是定值。由于是方程组，可以将看成一个向量。3.计算程序#include<stdio.h>#include<math.h>#defineN4/*宏定义N*/voidmain(){doubleh=0.0005,x=0,k[5][4],Y[4],y1[4],f[4];doubleF(doublex1,doubley[4],doublef[4]);/*调用函数F*/ doubleb[N]={0.025,0.045,0.085,0.1};inti,j,t,c;for(t=0;t<=N-1;t++)/*Runge-kutta开始*/ {c=int(b[t]/h);{for(j=0;j<=3;j++)Y[j]=0;/*初始条件为0*/for(i=1;i<=c;i++) {for(x=0;x<=b[t];x=x+h) for(j=1;j<=3;j++)/*求K1值*/ y1[j]=Y[j]; F(x,y1,f); for(j=1;j<=3;j++) k[1][j]=f[j];for(j=1;j<=3;j++)/*求K2值*/ y1[j]=Y[j]+0.5*h*k[1][j]; F(x,y1,f); for(j=1;j<=3;j++) k[2][j]=f[j]; for(j=1;j<=3;j++)/*求K3值*/ y1[j]=Y