数值分析课程作业

数值分析课程作业（1）解答如下：不稳定的数值计算公式往往会出现“差之毫厘，失之千里”的错误结果。因此，在计算过程中要选用稳定的计算公式。在本题的迭代计算中，由I0往后递推的求解In会使误差快速增大，而先近似计算In然后回代计算In-1…I0会使误差逐渐减小，因此使用此计算方法较为稳定。用C语言实现的代码如下：#includestdio.h/*计算k的阶乘*/intk_factorial(intk){if(0==k)return1;elsereturnk*k_factorial(k-1);}/*计算1/e，k控制计算精度*/doubleget_e_1(intk){if(0==k)return1;elsereturn1.0*((k%2)?-1:1)/k_factorial(k)+get_e_1(k-1);}/**该函数用于计算公式的近似值。对于给定的n值，首先计算(n+k)值对应的*公式的近似值，然后回代计算n值对应的近似值。*/doubleget_value(intn,intk){doubleres;res=(1+get_e_1(7))/((n+k)+1)/2;while(k0){res=(1-res)/(n+k);k--;}returnres;}intmain(intargc,char*argv[]){intn;printf(请输入一个n值：);scanf(%d,&n);if(n0)return0;printf(对于给定的n=%d，公式的近似值为：%lf\n,n,get_value(n,9));return0;}（2）解答如下：考虑拟合函数：)()()()(1100xqaxqaxqaxnn，将数据表xx1x2……xmf(x)y1x2……ym中的数据代入，得超定方程mmnnmmmnnnnyxqaxqaxqaxqayxqaxqaxqaxqayxqaxqaxqaxqa)()()()()()()()()()()()(2211002222221120011122111100（mn）其系数矩阵为)()()()()()()()()()()()(21022221201121110mnmmmnnxqxqxqxqxqxqxqxqxqxqxqxq由于多项式q0(x)，q1(x)，q2(x)，……，qn(x)在点集{x1，x2，……，xm}上的正交，所以超定方程组的系数矩阵中不同列的列向量是相互正交的向量组。于是用这一矩阵的转置矩阵去左乘超定方程组左、右两端得正规方程组),(),(),(),(),(),(11110000yqaqqyqaqqyqaqqnnnn=),/(),(),/(),(),/(),(11110000nnnnqqyqaqqyqaqqyqa其中，miikkkxqqq12)(),(，miiikkyxqyq1)(),(。因为正规方程组中每一个方程都是一元一次方程可以直接写出原超方程组的最小二乘解，所以拟合函数为)(),(),()(),(),()(),(),()(11110000xqqqyqxqqqyqxqqqyqxnnnn这一结果与用次多项式拟合所得结果在理论是完全一样的，只是形式上不同、算法实现上避免了解病态方程组。用C语言实现的代码如下：#includestdio.h#includestdlib.h#includemath.h#includeconio.hconstdoubleEPS=1E-10;//运算精度/***读入节点数组x、函数值数组y、权值数组ω及节点数N*/voidreadFile(double*&x,double*&y,double*&omiga,int&N){FILE*fp;fp=fopen(.\\ZXECF.txt,r);if(fp==NULL){printf(指定位置的文件不存在，请检查！！);getch();exit(0);}for(N=1;;N++){//计算节点个数Nfscanf(fp,%*lf%*lf%*lf);if(feof(fp))break;}x=newdouble[N];y=newdouble[N];omiga=newdouble[N];rewind(fp);for(inti=0;iN;i++)fscanf(fp,%lf%lf%lf,&x[i],&y[i],&omiga[i]);fclose(fp);}/***x数组存放节点的x[i]值，y数组存放节点的y[i]，omiga数组各节点对应权值*a数组为拟和的结果多项式的各系数，m为节点个数，n为多项式最高次数*/voidcurveFitting(double*x,double*y,double*omiga,double*a,intm,intn){double*alpha,*beta,*Q,*b,*t,*s,*d,*q;inti,j,k;doublesum1,sum2;alpha=newdouble[n+1];//α数组beta=newdouble[n+1];//β数组b=newdouble[n+1];//多项式Q(j-1)[x]的系数t=newdouble[n+1];//多项式Q(j)[x]的系数s=newdouble[n+1];//多项式Q(j+1)[x]的系数Q=newdouble[n+1];//临时值d[j]数组d=newdouble[n+1];//临时值d[j]数组q=newdouble[n+1];//临时值q[j]数组for(i=0;i=n;i++)//n次多项式的n+1个系数a[i]=0;if(nm)//拟和多项式的最高次数限定不超过mn=m;//构造Q0(x)Q[0]=b[0]=1;d[0]=sum1=sum2=0;for(i=0;im;i++){sum1+=omiga[i]*x[i]*pow(Q[0],2);d[0]+=omiga[i]*pow(Q[0],2);sum2+=omiga[i]*y[i]*Q[0];}q[0]=sum2/d[0];alpha[0]=sum1/d[0];a[0]=q[0]*b[0];//计算多项式系数a[0]//构造Q1(x)t[0]=-alpha[0];t[1]=1;d[1]=sum1=sum2=0;for(i=0;im;i++){Q[1]=t[0]+t[1]*x[i];sum1+=omiga[i]*x[i]*pow(Q[1],2);d[1]+=omiga[i]*pow(Q[1],2);sum2+=omiga[i]*y[i]*Q[1];}q[1]=sum2/d[1];alpha[1]=sum1/d[1];beta[1]=d[1]/d[0];a[0]+=q[1]*t[0];a[1]=q[1]*t[1];//计算多项式系数a[1]b[1]=t[1];for(j=2;j=n;j++){//构造Qj(x)//计算s[0]--s[j]s[j]=t[j-1];s[j-1]=-alpha[j-1]*t[j-1]+t[j-2];for(k=j-2;k=1;k--)s[k]=-alpha[j-1]*t[k]+t[k-1]-beta[j-1]*b[k];s[0]=-alpha[j-1]*t[0]-beta[j-1]*b[0];sum1=d[j]=sum2=0;for(k=0;km;k++){Q[j]=0;for(i=0;i=j;i++)Q[j]+=s[i]*pow(x[k],i);sum1+=omiga[k]*x[k]*pow(Q[j],2);d[j]+=omiga[k]*pow(Q[j],2);sum2+=omiga[k]*y[k]*Q[j];}q[j]=sum2/d[j];alpha[j]=sum1/d[j];beta[j]=d[j]/d[j-1];for(k=j-1;k=0;k--)//多项式拟和a[k]+=q[j]*s[k];a[j]=q[j]*s[j];for(k=j-1;k=0;k--){//向量赋值：T→B，S→Tb[k]=t[k];t[k]=s[k];}t[j]=s[j];}delete[]alpha;delete[]beta;delete[]s;delete[]t;delete[]b;delete[]Q;delete[]d;delete[]q;}intmain(void){double*x,*y;//x数组存放节点的x[i]值，y数组存放节点的y[i]double*omiga;//各节点对应权值数组double*a;//a数组存放n次多项式的n+1个系数intm,n,i;//m为节点个数，编号为0...m-1；n为多项式拟合次数readFile(x,y,omiga,m);//从文件中读取数据printf(原数据为：\n);printf(%-10s%-10s%-10s\n,x,y,节点权值);for(i=0;im;i++)printf(%-10g%-10g%-10g\n,x[i],y[i],omiga[i]);printf(\n节点个数为：m=%d\n\n,m);n=3;//3次多项式拟合a=newdouble[n+1];curveFitting(x,y,omiga,a,m,n);printf(%d次多项式拟和结果为：\n,n);printf(f(x)=%lg*x^%d,a[n],n);for(i=n-1;i=0;i--){if(fabs(a[i])1E-10)continue;if(a[i]0)printf(+);if(i1)printf(%lg*x^%d,a[i],i);elseif(i==1)printf(%lg*x,a[i]);elseif(i==0)printf(%lg,a[i]);}printf(\n\n);delete[]x;delete[]y;delete[]omiga;delete[]a;return0;}（3）解答如下：LU分解在本质上是高斯消元法的一种表达形式。实质上是将A通过初等行变换变成一个上三角矩阵，其变换矩阵就是一个单位下三角矩阵。这正是所谓的杜尔里特算法（Doolittlealgorithm）：从下至上地对矩阵A做初等行变换，将对角线左下方的元素变成零，然后再证明这些行变换的效果等同于左乘一系列单位下三角矩阵，这一系列单位下三角矩阵的乘积的逆就是L矩阵，它也是一个单位下三角矩阵。这类算法的复杂度一般在(三分之二的n三次方)左右。用C++实现的系数矩阵LU分解代码如下：#includeiostream.h#includestdlib.h#includemath.h#includefstream.h#defineAL100//局阵最大行#defineROW100//局阵最大列intmain(void){intn;//阶doubleA[AL][ROW];//系数矩阵doubleb[AL];//右端项inti,j;charflag='n';//选择标志/**********手工输入**********/cout输入阶:endl;cinn;cout输入系数矩阵A:endl;for(i=0;in;i++)for(j=0;jn;j++)cinA[i][j];cout输入右端项b:endl;for(i=0;in;i++)cinb[i];coutendl;/**********计算L矩阵和U矩阵(都放到A矩阵里)**********/intk;doublesum_1=0;for(k=0;kn-1;k++){//循环//L矩阵算法for(i=k+1;in;i++){sum_1=0;for(intm=0;mk;m++)sum_1+=A[i][m]*A[m][k];A[i][k]=(A[i][k]-sum_1)/A[k][k];}//U矩阵算法for(i=k+1;in;i++