北航数值分析B作业二
北航数值分析大作业二(纯原创,高分版)
(R_5 ,I_5 )=(-1.493147080915e+000, 0.000000000000e+000)
(R_6 ,I_6 )=(-9.891143464723e-001, 1.084758631502e-001)
-0.8945216982
-0.0993313649
-1.0998317589
0.9132565113
-0.6407977009
0.1946733679
-2.3478783624
2.3720579216
1.8279985523
-1.2630152661
0.6790694668
-0.4672150886
6.220134985374e-001
-1.119962139645e-001
-2.521344456568e+000
-1.306189420531e+000
-3.809101150714e+000
8.132800093357e+000
-1.230295627285e+000
-6.753086301215e-001
而其本质就是
1.令 以及最大迭代步数L；
2.若m≤0，则结束计算，已求出A的全部特征值，判断 或 或m≤2是否成立，成立则转3，否则转4；
3.若 ，则得一个特征值 ，m=m-1，降阶；若 ，则计算矩阵：
的特征值得矩阵A的两个特征值，m=m-2，降阶，转2.；
4.若k≤L，成立则令
k=k+1，转2，否则结束计算，为计算出矩阵A的全部特征值；
北航数值分析计算实习题目二 矩阵QR分解
数值分析实习二院（系）名称航空科学与工程学院专业名称动力工程及工程热物理学号SY0905303学生姓名解立垚1. 题目试用带双步位移QR 的分解法求矩阵A=[a ij ]10*10的全部特征值，并对其中的每一个实特征值求相应的特征向量。
已知()sin 0.50.2,1.5cos 1.2,ij i j i j a i j i j ⎧⎫+≠⎪⎪=⎨⎬+=⎪⎪⎩⎭(),1,2,...,10i j =。
说明：1、求矩阵特征值时，要求迭代的精度水平为1210ε-=。
2、打印以下内容：算法的设计方案；全部源程序(要求注明主程序和每个子程序的功能)； 矩阵A 经过拟上三角话之后所得的矩阵()1n A -；对矩阵()1n A-进行QR 分解方法结束后所得的矩阵；矩阵A 的全部特征值()(),1,2,......10i i iR I i λ=，和A 的相应于实特征值的特征向量；其中()(),.i e i m i R R I I λλ==如果i λ是实数，则令0.i I =3、采用e 型输出数据，并且至少显示12位有效数字。
2. 算法设计方案本题采用带双步位移的QR 分解方法。
为了使程序简洁，自定义类Xmatrix ，其中封装了所需要的函数方法。
在Xmatrix 类中封装了运算符重载的函数，即定义了矩阵的加、减、乘、除、数乘运算及转置运算(T())。
同时为了避免传递数组带来的额外内存开销，使用引用(&)代替值传递，以节省内存空间，避免溢出.（1）此程序的主要部分为Xmatrix 中的doubleQR()方法，具体如下：Step1：使用矩阵拟上三角化的算法将A 化为拟上三角阵A （n-1）(此处调用Xmatrix 中的preQR()方法)Step2：令121,,10k m n ε-===， 其中k 为迭代次数。
Step3：如果,1m m a ε-≤，则得到A 的一个特征值,m m a ，令1m m =-，goto Step4;否则goto Step5.Step4： 如果1m =，则得到A 的一个特征值11a ，goto Step11;如果0m =，则goto Step11;如果1m >，则goto Step3;Step5(Step6)：如果2m =，则得到A 的两个特征值12s s 和（12s s 和为右下角两阶子阵对应的特征方程21,1,()det 0m m m m a a D λλ---++=的两个根。
北航数理统计大作业2-聚类与判别分析
应用数理统计作业二学号：姓名：电话：二〇一四年十二月对NBA球队的聚类分析和判别分析摘要：NBA联盟作为篮球的最高殿堂深受广大球迷的喜爱，联盟的30支球队大家也耳熟能详，本文选取NBA联盟30支球队2013-2014常规赛赛季场均数据。
利用spss软件通过聚类分析对27个地区进行实力类型分类，并利用判断分析对其余3支球队对分类结果进行验证。
可以看出各球队实力类型与赛季实际结果相吻合。
关键词：聚类分析，判别分析，NBA目录1. 引言 (4)2、相关统计基础理论 (5)2.1、聚类分析 (5)2.2，判别分析 (6)3.聚类分析 (7)3.1数据文件 (7)3.2聚类分析过程 (9)3.3 聚类结果分析 (11)4、判别分析 (12)4.1 判别分析过程 (12)4.2判别检验 (17)5、结论 (20)参考文献 (21)致谢 (22)1. 引言1896年，美国第一个篮球组织"全国篮球联盟（简称NBL）"成立，但当时篮球规则还不完善，组织机构也不健全，经过几个赛季后，该组织就名存实亡了。
1946年4月6日，由美国波士顿花园老板沃尔特.阿.布朗发起成立了“美国篮球协会”（简称BAA）。
1949年在布朗的努力下，美国两大篮球组织BAA和NBL合并为“全国篮球协会”（简称NBA）。
NBA季前赛是 NBA各支队伍的热身赛，因为在每个赛季结束后，每支球队在阵容上都有相当大的变化，为了让各队磨合阵容，熟悉各自球队的打法，确定各队新赛季的比赛阵容、同时也能增进队员、教练员之间的沟通，所以在每个赛季开始之前，NBA就举办若干场季前赛，使他们能以比较好的状态投入到漫长的常规赛的比赛当中。
为了扩大NBA在全球的影响，季前赛有约三分之一的球队在美国以外的国家举办。
从总体上看，NBA的赛程安排分为常规赛、季后赛和总决赛。
常规赛采用主客场制，季后赛和总决赛采用七场四胜制的淘汰制。
[31]NBA常规赛从每年的11月的第一个星期二开罗，到次年的4月20日左右结束。
北航数值分析第二次大作业--QR分解
《数值分析A》计算实习题目二姓名学号联系方式班级指导教师2012年10月一、算法设计方案整个程序主要分为四个函数，主函数，拟上三角化函数，QR分解函数以及使用双步位移求解矩阵特征值、特征向量的函数。
因为在最后一个函数中也存在QR分解，所以我没有采用参考书上把矩阵M进行的QR分解与矩阵Ak的迭代合并的方法，而是在该函数中调用了QR分解函数，这样增强了代码的复用性，减少了程序长度；但由于时间关系，对阵中方法的运算速度没有进行深入研究。
1.为了减少QR分解法应用时的迭代次数，首先对给定矩阵进行拟上三角化处理。
2.对经过拟上三角化处理的矩阵进行QR分解。
3.注意到计算特征值与特征向量的过程首先要应用前面两个函数，于是在拟上三角化矩阵的基础上对QR分解函数进行了调用。
计算过程中，没有采用goto语句，而是根据流程图采用其他循环方式完成了设计，通过对迭代过程的合并，简化了程序的循环次数，最后在计算特征向量的时候采用了列主元高斯消去法。
二、源程序代码#include<stdio.h>#include<math.h>#include<string.h>int i,j,k,l,m; //定义外部变量double d,h,b,c,t,s;double A[10][10],AA[10][10],R[10][10],Q[10][10],RQ[10][10]; double X[10][10],Y[10][10],Qt[10][10],M[10][10];double U[10],P[10],T[10],W[10],Re[10]={0},Im[10]={0}; double epsilon=1e-12;void main(){void Quasiuppertriangular(double A[][10]);void QRdecomposition(double A[][10]);void DoublestepsQR(double A[][10]);int i,j;for(i=0;i<10;i++){for(j=0;j<10;j++){A[i][j]=sin(0.5*(i+1)+0.2*(j+1));Q[i][j]=0;AA[i][j]=A[i][j];}A[i][i]=1.5*cos(2.2*(i+1));AA[i][i]=A[i][i];}Quasiuppertriangular(A); //调用拟上三角化函数printf( "\n A经过拟上三角化矩阵为：\n\n");for(i=0;i<10;i++) //输出拟上三角化矩阵{for(j=0;j<10;j++){printf("%.12e ",A[i][j]); //输出拟上三角化矩阵}printf( "\n\n");}QRdecomposition(A); //调用QR分解函数printf( " 进行QR分解后，R矩阵为：\n\n"); //输出R矩阵for(i=0;i<10;i++){for(j=0;j<10;j++){printf("%.12e ",R[i][j]);}printf( "\n\n");}printf( " Q矩阵为:\n\n"); //输出Q矩阵for(i=0;i<10;i++){for(j=0;j<10;j++){printf("%.12e ",Q[i][j]);}printf( "\n\n");}printf( " RQ矩阵为：\n\n"); //输出RQ矩阵for(i=0;i<10;i++){for(j=0;j<10;j++){printf("%.12e ",RQ[i][j]);}printf( "\n\n");}DoublestepsQR(A); //调用双步位移函数printf( "\n\n 特征值实部依次为：\n\n"); //输出特征值实部for(j=0;j<10;j++){printf("%.12e ",Re[j]);}printf("\n\n 特征值虚部依次为：\n\n "); //输出特征值虚部for(j=0;j<10;j++){printf("%.12e ",Im[j]);}//按行输出特征向量printf( "\n\n 按行输出实特征根相应特征向量为：\n\n");for(i=0;i<10;i++){if(i==1||i==2||i==5||i==6){continue;}for(j=0;j<10;j++){printf("%.12e ",X[i][j]);}printf( "\n\n");}getchar();}//拟上三角化函数void Quasiuppertriangular(double A[][10]) {for(j=0;j<8;j++){for(i=0;i<10;i++){U[i]=0;P[i]=0;T[i]=0;W[i]=0;}m=0;for(i=j+2;i<10;i++){if(A[i][j]!=0){m=m+1;}}if(m==0){continue;}d=0;for(i=j+1;i<10;i++){d=d+pow(A[i][j],2);}d=sqrt(d);c=-d;if(A[j+1][j]<=0){c=d;}h=c*(c-A[j+1][j]);U[j+1]=A[j+1][j]-c;for(i=j+2;i<10;i++){U[i]=A[i][j];}for(i=0;i<10;i++){for(k=0;k<10;k++){P[i]=P[i]+U[k]*A[k][i];}P[i]=P[i]/h;}t=0;for(i=0;i<10;i++){for(k=0;k<10;k++){T[i]=T[i]+U[k]*A[i][k];}T[i]=T[i]/h;t=t+P[i]*U[i];}t=t/h;for(i=0;i<10;i++){W[i]=T[i]-t*U[i];for(k=0;k<10;k++){A[i][k]=A[i][k]-W[i]*U[k]-U[i]*P[k];if(abs(A[i][k])<1e-12){A[i][k]=0;}}}}}//QR分解函数void QRdecomposition(double A[][10]) {for(i=0;i<10;i++){for(j=0;j<10;j++){RQ[i][j]=0;Q[i][j]=0;R[i][j]=A[i][j];}Q[i][i]=1;}for(j=0;j<9;j++){for(i=0;i<10;i++){U[i]=0;P[i]=0;W[i]=0;}m=0;for(i=j+1;i<10;i++){if(R[i][j]!=0){m=m+1;}}if(m==0){continue;}d=0;for(i=j;i<10;i++){d=d+pow(R[i][j],2);}d=sqrt(d);c=-d;if(R[j][j]<=0){c=d;}h=c*(c-R[j][j]);U[j]=R[j][j]-c;for(i=j+1;i<10;i++){U[i]=R[i][j];}for(i=0;i<10;i++){for(k=0;k<10;k++){W[i]=W[i]+U[k]*Q[i][k];}}for(i=0;i<10;i++){for(k=0;k<10;k++){Q[i][k]=Q[i][k]-((W[i]*U[k])/h);}}for(i=0;i<10;i++){for(k=0;k<10;k++){P[i]=P[i]+U[k]*R[k][i];}P[i]=P[i]/h;}for(i=0;i<10;i++){for(k=0;k<10;k++){R[i][k]=R[i][k]-U[i]*P[k];if(abs(R[i][k])<epsilon){R[i][k]=0;}}}}for(i=0;i<10;i++) //计算A(n+1)=RQ {for(j=0;j<10;j++){for(k=0;k<10;k++){RQ[i][j]=RQ[i][j]+R[i][k]*Q[k][j];}}}}//双步位移法计算特征值特征向量函数void DoublestepsQR(double A[][10]){int L=1000,m=9; //定义最大循环次数for(i=0;i<L;i++){for(;m>-1;){if(abs(A[m][m-1])<=epsilon){Re[m]=A[m][m];m=m-1; //降阶if(m==0) //4{Re[0]=A[0][0];break;}if(m==-1){break;}if(m>1){continue;}}b=-A[m][m]-A[m-1][m-1]; //5c=A[m][m]*A[m-1][m-1]-A[m][m-1]*A[m-1][m];if(m==1) //6{if((b*b-4*c)>=0){Re[m]=(-b+sqrt(b*b-4*c))/2;Re[m-1]=(-b-sqrt(b*b-4*c))/2;}if((b*b-4*c)<0){Re[m]=-b/2; Im[m]=sqrt(4*c-b*b)/2;Re[m-1]=-b/2; Im[m-1]=-sqrt(4*c-b*b)/2;}m=m-1; //循环出口条件break;}if((m>1)&&(abs(A[m-1][m-2])>epsilon)) //8{if(i==L-1){printf("No results! \n");m=0; //循环出口条件break;}break;}if((m>1)&&(abs(A[m-1][m-2])<=epsilon)) //7 {if((b*b-4*c)>0){Re[m]=(-b+sqrt(b*b-4*c))/2;Re[m-1]=(-b-sqrt(b*b-4*c))/2;}if((b*b-4*c)<0){Re[m]=-b/2; Im[m]=sqrt(4*c-b*b)/2;Re[m-1]=-b/2; Im[m-1]=-sqrt(4*c-b*b)/2;}m=m-2; //降阶if(m>0){continue;}if(m==0){Re[0]=A[0][0];break;}}}if(m<=0){break;}s=A[m-1][m-1]+A[m][m]; //9t=A[m][m]*A[m-1][m-1]-A[m][m-1]*A[m-1][m];for(j=0;j<10;j++){for(k=0;k<10;k++){Qt[j][k]=0;Q[j][k]=0;M[j][k]=0;X[j][k]=0;Y[j][k]=0;}}for(j=0;j<m+1;j++){for(k=0;k<m+1;k++){for(l=0;l<m+1;l++){M[j][k]=M[j][k]+A[j][l]*A[l][k];}}}for(j=0;j<m+1;j++){for(k=0;k<m+1;k++){M[j][k]=M[j][k]-s*A[j][k];}M[j][j]=M[j][j]+t;}//调用QR分解函数对M矩阵进行分解并传递参数矩阵QQRdecomposition(M);for(j=0;j<10;j++){for(k=0;k<10;k++){Qt[j][k]=Q[k][j];}}for(j=0;j<m+1;j++){for(k=0;k<m+1;k++){for(l=0;l<m+1;l++){X[j][k]=X[j][k]+Qt[j][l]*A[l][k];}}}for(j=0;j<m+1;j++){for(k=0;k<m+1;k++){for(l=0;l<m+1;l++){Y[j][k]=Y[j][k]+X[j][l]*Q[l][k];}}}for(j=0;j<10;j++){{A[j][k]=Y[j][k];}}}//应用列主元高斯消元法计算实部特征向量for(l=0;l<10;l++){if(l==1||l==2||l==5||l==6){continue;}for(k=0;k<10;k++){for(m=0;m<10;m++){A[k][m]=AA[k][m];}A[k][k]=A[k][k]-Re[l];}for(j=0;j<9;j++){m=j;for(i=j+1;i<10;i++){if(abs(A[i][j])>abs(A[m][j])){m=i;}}{Y[j][k]=A[j][k];A[j][k]=A[m][k];A[m][k]=Y[j][k];}for(k=j+1;k<10;k++){b=A[k][j]/A[j][j];for(i=j;i<10;i++){A[k][i]=A[k][i]-A[j][i]*b;}}}X[l][9]=1;for(i=8;i>=0;i--){c=0;for(j=i+1;j<10;j++){c=c+A[i][j]*X[l][j];}X[l][i]=-c/A[i][i];}}}三、程序输出结果1819。
北航数值分析大作业二
数值分析第二题1 算法设计方案要想得出该题的答案首先要将矩阵A 进行拟上三角化，把矩阵A 进行QR 分解。
要得出矩阵A 的全部特征值先对A 进行QR 的双步位移得出特征值。
采用列主元的高斯消元法求解特征向量。
1.1 A 的拟上三角化因为对矩阵进行QR 分解并不改变矩阵的结构，因此在进行QR 分解前对矩阵A 进行拟上三角化可以大大减少计算机的计算量，提高程序的运行效率。
具体算法如下所示，记A A =)1(，并记)(r A 的第r 列至第n 列的元素为()n r r j n i a r ij,,1,;,,2,1)( +==。
对于2,,2,1-=n r 执行1. 若()n r r i a r ir,,3,2)( ++=全为零，则令)()1(r r A A =+,转5；否则转2。
2. 计算()∑+==nr i r irr a d 12)(()()r r r r r r r r r r d c a d a c ==-=++则取,0sgn )(,1)(,1若 )(,12r rr r r r a c c h +-=3. 令()n Tr nrr r r r r r r r R a a c a u ∈-=++)()(,2)(,1,,,,0,,0 。
4. 计算r r T r r h u A p /)(= r r r r h u A q /)(=r r Tr r h u p t /=r r r r u t q -=ωTrr T r r r r p u u A A --=+ω)()1(5. 继续。
1.2 A 的QR 分解具体算法如下所示，记)1(1-=n A A ，并记[]nn r ij r a A ⨯=)(，令I Q =1 对于1,,2,1-=n r 执行1.若()n r r i a r ir,,3,1)( ++=全为零，则令r r Q Q =+1r r A A =+1,转5；否则转2。
2.计算()∑==nri r irr a d 2)(()()r r r r r r r r r r d c a d a c ==-=++则取,0sgn )(,1)(,1若)(,2r r r r r r a c c h -=3.令()n Tr nrr r r r r r r r R a a c a u ∈-=+)()(,2)(,,,,,0,,0 。
北航数值分析全部三次大作业
北航数值分析全部三次大作业第一次大作业是关于解线性方程组的数值方法。
我们被要求实现各种常用的线性方程组求解算法，例如高斯消元法、LU分解法和迭代法等。
我首先学习了这些算法的原理和实现方法，并借助Python编程语言编写了这些算法的代码。
在实验中，我们使用了不同规模和条件的线性方程组进行测试，并比较了不同算法的性能和精度。
通过这个作业，我深入了解了线性方程组求解的原理和方法，提高了我的编程和数值计算能力。
第二次大作业是关于数值积分的方法。
数值积分是数值分析中的重要内容，它可以用于计算曲线的长度、函数的面积以及求解微分方程等问题。
在这个作业中，我们需要实现不同的数值积分算法，例如矩形法、梯形法和辛普森法等。
我学习了这些算法的原理和实现方法，并使用Python编写了它们的代码。
在实验中，我们计算了不同函数的积分值，并对比了不同算法的精度和效率。
通过这个作业，我深入了解了数值积分的原理和方法，提高了我的编程和数学建模能力。
第三次大作业是关于常微分方程的数值解法。
常微分方程是数值分析中的核心内容之一，它可以用于描述众多物理、化学和生物现象。
在这个作业中，我们需要实现不同的常微分方程求解算法，例如欧拉法、龙格-库塔法和Adams法等。
我学习了这些算法的原理和实现方法，并使用Python编写了它们的代码。
在实验中，我们解决了一些具体的常微分方程问题，并比较了不同算法的精度和效率。
通过这个作业，我深入了解了常微分方程的原理和方法，提高了我的编程和问题求解能力。
总的来说，北航数值分析课程的三次大作业非常有挑战性，但也非常有意义。
通过这些作业，我在数值计算和编程方面得到了很大的提升，也更加深入地了解了数值分析的理论和方法。
虽然这些作业需要大量的时间和精力，但我相信这些努力将会对我未来的学习和工作产生积极的影响。
北航数值分析大作业 第二题 QR分解
《数值分析B》课计算实习第一题设计文档与源程序姓名：杨彦杰学号：SY10171341 算法的设计方案（1）运行平台操作系统：Windows XP；开发平台：VC6.0++；工程类型：文档视图类；工程名：Numanalysis；（2）开发描述首先新建类CMetrix，该类完成矩阵之间的相关运算，包括相乘、加减等，以主程序方便调用；题目的解算过程在视图类CNumanalysisView中实现，解算结果在视图界面中显示；（3）运行流程（4）运行界面2、全部源代码（1）类CMetrixMetrix.h文件：class CMetrix{public:double** MetrixMultiplyConst(double**A,int nRow,int nCol,double nConst);//矩阵乘常数double** MetrixMultiplyMetrix(double**A,double**mA,int nRow,int nCol);//矩阵相乘double** MetrixSubtractMetrix(double **A, double **subA, int nRow,int nCol);//矩阵减矩阵double VectorMultiplyVector(double*V,double*mulV,int nV);//向量点积double** VectorMultiplyVectortoMetrix(double*V,double*VT,int nV);//向量相乘为矩阵double* VectorSubtractVector(double*V,double*subV,int nV);//向量相减double* VectorMultiplyConst(double *V, int nV, double nConst);//向量乘常数double LengthofVector(double *V,int nV);//求向量的长度double* MetrixMultiplyVector(double**A,int nRow,int nCol,double*V,int nV);//矩阵与向量相乘double** AtoAT(double **A,int Row,int Col);//矩阵转置运算void FreeMem();CMetrix(int nRow,int nCol);uCMetrix();virtual ~CMetrix();double* vector; //过渡向量double** B; //过渡矩阵};Metrix.cpp文件：CMetrix::CMetrix(int nRow, int nCol){B = new double*[nRow];for (int i = 0;i < nCol;i++){B[i] = new double[nCol];}vector = new double[nRow];}CMetrix::~CMetrix(){delete vector;B = NULL;delete B;}double** CMetrix::AtoAT(double **A, int nRow, int nCol){for (int row = 0;row < nRow;row++){for (int col = 0;col < nCol;col++){B[col][row] = A[row][col];}}return B;}double* CMetrix::MetrixMultiplyVector(double **A, int nRow, int nCol, double *V, int nV) {if (nCol != nV){AfxMessageBox("矩阵列数和向量维数不等,不能相乘!");return 0;}double sum = 0.0;for (int row = 0;row < nRow;row++){for (int col = 0;col < nCol;col++){sum += A[row][col]*V[col];}vector[row] = sum;sum = 0.0;}return vector;}double CMetrix::LengthofVector(double *V, int nV){double length = 0.0;for (int col = 0;col < nV;col++){length += V[col]*V[col];}return length;}double* CMetrix::VectorMultiplyConst(double *V, int nV, double nConst){for (int col = 0;col < nV;col++){vector[col] = V[col]*nConst;}return vector;}double* CMetrix::VectorSubtractVector(double *V, double *subV, int nV){for (int col = 0;col < nV;col++){vector[col] = V[col]-subV[col];}return vector;}double** CMetrix::VectorMultiplyVectortoMetrix(double*V, double *VT, int nV){for (int row = 0;row < nV;row++){for (int col = 0;col < nV;col++){B[row][col] = V[row]*VT[col];}}return B;}double CMetrix::VectorMultiplyVector(double *V, double *mulV, int nV){double length = 0.0;for (int col = 0;col < nV;col++){length += V[col]*mulV[col];}return length;}double** CMetrix::MetrixSubtractMetrix(double **A, double **subA, int nRow, int nCol) {for (int row = 0;row < nRow;row++){for (int col = 0;col < nCol;col++){B[row][col] = A[row][col]-subA[row][col];}}return B;}double** CMetrix::MetrixMultiplyMetrix(double **A, double **mA, int nRow, int nCol) {double sum = 0.0;for (int row = 0;row < nRow;row++){for (int col = 0;col < nCol;col++){for(int n = 0;n < nCol;n++){sum += A[row][n]*mA[n][col];}B[row][col] = sum;sum = 0.0;}}return B;}double** CMetrix::MetrixMultiplyConst(double **A, int nRow, int nCol, double nConst) {for (int row = 0;row < nRow;row++){for (int col = 0;col < nCol;col++){B[row][col] = A[row][col]*nConst;}}return B;}（2）类CNumanalysisViewNumanalysisview.hclass CNumanalysisView : public CEditView{…………public:double Sign(double x);void DisplayVector(double*V,int nV); // 显示向量数据void DisplayMetrix(double **A,int Row,int Col); //显示矩阵void DisplayText(CString str); //显示文本protected://{{AFX_MSG(CNumanalysisView)afx_msg void OnQRanalyze(); //运行主函数…………};Numanalysisview.cppvoid CNumanalysisView::OnQRanalyze(){//开辟空间int nRow = 10;int nCol = 10;CString str;CMetrix Metrix(nRow,nCol);double tempa = 0.0;double *V = new double[nCol]; //分配10*10矩阵空间double *ur = new double[nCol];double *pr = new double[nCol];double *qr = new double[nCol];double *wr = new double[nCol];double *tempV = new double[nCol];double **Ar = new double*[nRow];double **C = new double*[nRow];double **Cr = new double*[nRow];double **tempA = new double*[nRow];double **A = new double*[nRow];double **R = new double*[nRow];for (int col = 0;col < nRow;col++){A[col] = new double[nCol];Ar[col] = new double[nCol];C[col] = new double[nCol];Cr[col] = new double[nCol];tempA[col] = new double[nCol];R[col] = new double[nCol];}//矩阵A求解for (int i = 0;i < nRow;i++){for (int j = 0;j < nCol;j++){if(i == j)A[i][j] = 1.5*cos((i+1.0)+1.2*(j+1.0));elseA[i][j] = sin(0.5*(i+1.0)+0.2*(j+1.0));}}//--------------------拟上三角化-------------------------// double dr = 0.0,cr = 0.0,hr = 0.0,tr = 0.0;for (int r = 0;r < nCol - 2;r++){dr = 0.0;for (i = r+1;i < nCol;i++) //dr{dr += A[i][r]*A[i][r];}dr = sqrt(dr);for (i = r+2;i < nCol;i++) //判断air是否全为零tempa += fabs(A[i][r]);if (tempa <= IPSLEN)continue;if (A[r+1][r] == 0.0) //crcr = dr;elsecr = -1*Sign(A[r+1][r])*dr;hr = cr*cr - cr*A[r+1][r]; //hrstr.Format("dr = %.6e, cr = %.6e, hr = %.6e",dr,cr,hr);for (int row = 0;row < nRow;row++) //ur{if (row < r+1)ur[row] = 0.0;else if (row == r+1)ur[row] = A[row][r]-cr;elseur[row] = A[row][r];}tempA = Metrix.AtoAT(A,nRow,nCol);for (row = 0;row < nRow;row++){for (col = 0;col < nCol;col++)Ar[row][col] = tempA[row][col];}tempV = Metrix.MetrixMultiplyVector(Ar,nRow,nCol,ur,nCol); //pr memcpy(pr,tempV,nCol*8);tempV = Metrix.VectorMultiplyConst(pr,nCol,1.0/hr);memcpy(pr,tempV,nCol*8);tempV = Metrix.MetrixMultiplyVector(A,nRow,nCol,ur,nCol); //qr memcpy(qr,tempV,nCol*8);tempV = Metrix.VectorMultiplyConst(qr,nCol,1.0/hr);memcpy(qr,tempV,nCol*8);tr = Metrix.VectorMultiplyVector(pr,ur,nCol)/hr; //trtempV = Metrix.VectorMultiplyConst(ur,nCol,tr); //wr memcpy(wr,tempV,nCol*8);tempV = Metrix.VectorSubtractVector(qr,wr,nCol);memcpy(wr,tempV,nCol*8);tempA = Metrix.VectorMultiplyVectortoMetrix(wr,ur,nCol); //Arfor (row = 0;row < nRow;row++){for (col = 0;col < nCol;col++)Ar[row][col] = tempA[row][col];}tempA = Metrix.MetrixSubtractMetrix(A,Ar,nRow,nCol);for (row = 0;row < nRow;row++){for (col = 0;col < nCol;col++)A[row][col] = tempA[row][col];}tempA = Metrix.VectorMultiplyVectortoMetrix(ur,pr,nCol);for (row = 0;row < nRow;row++){for (col = 0;col < nCol;col++)Ar[row][col] = tempA[row][col];}tempA = Metrix.MetrixSubtractMetrix(A,Ar,nRow,nCol);for (row = 0;row < nRow;row++){for (col = 0;col < nCol;col++){A[row][col] = tempA[row][col];if (fabs(A[row][col]) < IPSLEN){A[row][col] = 0.0;}}}}DisplayText("矩阵A拟上三角化后所得的矩阵为：");DisplayMetrix(A,nRow,nCol);for (int row = 0;row < nRow;row++) //用于计算特征向量{for (col = 0;col < nCol;col++)R[row][col] = A[row][col];}// -------------------------------------------------////--------------------带双步位移的QR分解-------------------------// int m = nCol;struct EigenVal //定义特征值结构，实数和虚数{double Realnum;double Imagnum;};EigenVal *eigenvalue = new EigenVal[m];EigenVal tmpEigen1,tmpEigen2;double b = 0.0,c = 0.0,delta = 0.0,s = 0.0,t = 0.0;double *vr = new double[m];for (int k = 1;k < 100; k++){//m代表矩阵阶数，判断式中直接用，运算中需要-1while (m > 1 && fabs(A[m-1][m-2]) <= IPSLEN)//第三步和第四步{eigenvalue[m-1].Realnum = A[m-1][m-1];eigenvalue[m-1].Imagnum = 0.0;m = m - 1;}if (m == 1){eigenvalue[m-1].Realnum = A[m-1][m-1];eigenvalue[m-1].Imagnum = 0.0;DisplayText("已求出A的全部特征值：");break;}b = -(A[m-2][m-2]+A[m-1][m-1]); //第五步求一元二次方程式的根s1,s2c = A[m-2][m-2]*A[m-1][m-1]-A[m-2][m-1]*A[m-1][m-2];delta =b*b - 4*c;if (delta >= 0.0){tmpEigen1.Realnum = (-b-sqrt(delta))/2;tmpEigen1.Imagnum = 0.0;tmpEigen2.Realnum = (-b+sqrt(delta))/2;tmpEigen2.Imagnum = 0.0;}else{tmpEigen1.Realnum = -b/2;tmpEigen1.Imagnum = -sqrt(fabs(delta))/2 ;tmpEigen2.Realnum = -b/2;tmpEigen2.Imagnum = sqrt(fabs(delta))/2;}if (m == 2) //第六步 m=2时结束运算{eigenvalue[m-1] = tmpEigen1;eigenvalue[m-2] = tmpEigen2;DisplayText("已求出A的全部特征值：");break;}else //第七步 m > 1{if (fabs(A[m-2][m-3]) <= IPSLEN){eigenvalue[m-1] = tmpEigen1;eigenvalue[m-2] = tmpEigen2;m = m - 2;continue;}}for (int row = 0;row < m;row++) //Mk求之前需要把A付给C{for (int col = 0;col < m;col++)C[row][col] = A[row][col];}double **I = new double*[m]; //第九步求Mk和Mk的QR分解for (int i = 0;i < m;i++) //求单位矩阵I,分配m*m矩阵空间{I[i] = new double[m];}for (i = 0;i < m;i++){for (int j = 0;j < m;j++){if(i == j)I[i][j] = 1;else I[i][j] = 0;}}s = A[m-2][m-2]+A[m-1][m-1];t = A[m-2][m-2]*A[m-1][m-1] - A[m-2][m-1]*A[m-2][m-1];tempA = Metrix.MetrixMultiplyMetrix(A,A,m,m);//A*Afor (row = 0;row < m;row++){for (col = 0;col < m;col++)Ar[row][col] = tempA[row][col];}tempA = Metrix.MetrixMultiplyConst(A,m,m,s);//s*Afor (row = 0;row < m;row++){for (col = 0;col < m;col++)A[row][col] = tempA[row][col];}tempA = Metrix.MetrixSubtractMetrix(Ar,A,m,m);//A*A-s*Afor (row = 0;row < m;row++){for (col = 0;col < m;col++)A[row][col] = tempA[row][col]; }tempA = Metrix.MetrixMultiplyConst(I,m,m,-1*t);//-t*Ifor (row = 0;row < m;row++){for (col = 0;col < m;col++)Ar[row][col] = tempA[row][col]; }tempA = Metrix.MetrixSubtractMetrix(A,Ar,m,m);//A*A - s*A + r*I for (row = 0;row < m;row++){for (col = 0;col < m;col++){A[row][col] = tempA[row][col];if (fabs(A[row][col]) < IPSLEN){A[row][col] = 0.0;}}}delete I;//Mk的QR分解for (int r = 0;r < m - 1;r++){dr = 0.0;for (i = r;i < m;i++) //dr{dr += A[i][r]*A[i][r];}dr = sqrt(dr);for (i = r+1;i < m;i++) //判断air是否全为零tempa += fabs(A[i][r]);if (tempa <= IPSLEN)continue;if (A[r][r] == 0.0) //crcr = dr;elsecr = -1*Sign(A[r][r])*dr;hr = cr*cr - cr*A[r][r]; //hrfor (int row = 0;row < m;row++) //ur{if (row < r)ur[row] = 0.0;else if (row == r)ur[row] = A[row][r]-cr;elseur[row] = A[row][r];}tempA = Metrix.AtoAT(A,m,m); //Btfor (row = 0;row < m;row++){for (col = 0;col < m;col++)Ar[row][col] = tempA[row][col];}tempV = Metrix.MetrixMultiplyVector(Ar,m,m,ur,m); //Bt*ur memcpy(vr,tempV,m*8);tempV = Metrix.VectorMultiplyConst(vr,m,1.0/hr); //vr = Bt*ur/hr memcpy(vr,tempV,m*8);tempA = Metrix.VectorMultiplyVectortoMetrix(ur,vr,m);//Ur*vrfor (row = 0;row < m;row++){for (col = 0;col < m;col++)Ar[row][col] = tempA[row][col];}tempA = Metrix.MetrixSubtractMetrix(A,Ar,m,m); //Br-ur*vrfor (row = 0;row < m;row++){for (col = 0;col < m;col++){A[row][col] = tempA[row][col];if (fabs(A[row][col]) < IPSLEN){A[row][col] = 0.0;}}}tempA = Metrix.AtoAT(C,m,m); //Ctfor (row = 0;row < m;row++){for (col = 0;col < m;col++)Cr[row][col] = tempA[row][col]; }tempV = Metrix.MetrixMultiplyVector(Cr,m,m,ur,m); //pr memcpy(pr,tempV,m*8);tempV = Metrix.VectorMultiplyConst(pr,m,1.0/hr);memcpy(pr,tempV,m*8);tempV = Metrix.MetrixMultiplyVector(C,m,m,ur,m); //qr memcpy(qr,tempV,m*8);tempV = Metrix.VectorMultiplyConst(qr,m,1.0/hr);memcpy(qr,tempV,m*8);tr = Metrix.VectorMultiplyVector(pr,ur,m)/hr; //trtempV = Metrix.VectorMultiplyConst(ur,m,tr); //wr memcpy(wr,tempV,m*8);tempV = Metrix.VectorSubtractVector(qr,wr,m);memcpy(wr,tempV,m*8);tempA = Metrix.VectorMultiplyVectortoMetrix(wr,ur,m);//Cr+1for (row = 0;row < m;row++){for (col = 0;col < m;col++)Cr[row][col] = tempA[row][col]; }tempA = Metrix.MetrixSubtractMetrix(C,Cr,m,m);for (row = 0;row < m;row++){for (col = 0;col < m;col++)C[row][col] = tempA[row][col]; }tempA = Metrix.VectorMultiplyVectortoMetrix(ur,pr,m);for (row = 0;row < m;row++){for (col = 0;col < m;col++)Cr[row][col] = tempA[row][col]; }tempA = Metrix.MetrixSubtractMetrix(C,Cr,m,m);for (row = 0;row < m;row++){for (col = 0;col < m;col++){C[row][col] = tempA[row][col];if (fabs(C[row][col]) < IPSLEN){C[row][col] = 0.0;}}}}str.Format("矩阵A%d QR分解结束后所得到的矩阵为：",m);//计算结果输出DisplayText(str);DisplayMetrix(A,m,m);for (row = 0;row < m;row++) //Mk的QR分解后需要把C付给A{for (col = 0;col < m;col++)A[row][col] = C[row][col];}str.Format("迭代完成后的矩阵A%d = ",k);DisplayText(str);DisplayMetrix(A,m,m);}DisplayText("矩阵A的全体特征值如下: ");for (i = 0;i<nCol;i++){str.Format("%.6e + j%.6e",eigenvalue[i].Realnum,eigenvalue[i].Imagnum);DisplayText(str);}// -------------------------------------------------//求实特征值的特征向量，在拟上三角矩阵基础上直接求解即可////(A-egiI)X = 0.0;m = nRow;for (row = 0;row < nRow;row++) //用于计算特征向量{for (col = 0;col < nCol;col++)A[row][col] = R[row][col];}double **I = new double*[m]; //求单位矩阵I,分配m*m矩阵空间double sum = 0.0;for (i = 0;i < m;i++){I[i] = new double[m];}for (i = 0;i < m;i++){for (int j = 0;j < m;j++){if(i == j)I[i][j] = 1;else I[i][j] = 0;}}for (i = 0;i < nRow;i++){if (eigenvalue[i].Imagnum != 0.0){str.Format("特征值%.6e+j%.6e为虚数，不需要求特征向量。
北航2010-2015年研究生数值分析报告期末模拟试卷与真题
北航2010-2015年研究生数值分析报告期末模拟试卷与真题数值分析模拟卷A一、填空（共30分，每空3分）1 设-=1511A ，则A 的谱半径=)(a ρ______，A 的条件数)(1A cond =________. 2 设 ,2,1,0,,53)(2==+=k kh x x x f k ，则],,[21++n n n x x x f ＝________, ],,[321+++n n n n x x x x f ，＝________.3 设≤≤-++≤≤+=21,1210,)(2323x cx bx x x x x x S ，是以0，1，2为节点的三次样条函数，则b=________,c=________.4 设∞=0)]([k k x q 是区间［0，1］上权函数为x x =)(ρ的最高项系数为1的正交多项式族，其中1)(0=x q ，则?=10)(dx x xq k ________，=)(2x q ________.5 设=11001a a a a A ，当∈a ________时，必有分解式，其中L 为下三角阵，当其对角线元素)3,2,1(=i L ii 满足条件________时，这种分解是唯一的.二、（14分）设49,1,41,)(21023====x x x x x f , （1）试求)(x f 在]49,41[上的三次Hermite 插值多项式)(x H 使满足2,1,0),()(==i x f x H i i ，)()(11x f x H '='.（2）写出余项)()()(x H x f x R -=的表达式.三、（14分）设有解方程0cos 2312=+-x x 的迭代公式为n n x x cos 3241+=+，（1）证明R x ∈?0均有?∞→=x x n x lim （?x 为方程的根）；（2）取40=x ，用此迭代法求方程根的近似值，误差不超过，列出各次迭代值；（3）此迭代的收敛阶是多少？证明你的结论.四、(16分) 试确定常数A ，B ，C 和，使得数值积分公式有尽可能高的代数精度. 试问所得的数值积分公式代数精度是多少？它是否为Gauss 型的？五、（15分）设有常微分方程的初值问题=='00)(),(y x y y x f y ，试用Taylor 展开原理构造形如)()(11011--++++=n n n n n f f h y y y ββα的方法，使其具有二阶精度，并推导其局部截断误差主项.六、（15分）已知方程组b Ax ＝，其中= ??=21,13.021b A ，（1）试讨论用Jacobi 迭代法和Gauss-Seidel 迭代法求解此方程组的收敛性.（2）若有迭代公式)()()()1(b Ax a x x k k k ++=+，试确定一个的取值围，在这个围任取一个值均能使该迭代公式收敛.七、（8分）方程组，其中，A 是对称的且非奇异.设A 有误差，则原方程组变化为，其中为解的误差向量，试证明 .其中1λ和2λ分别为A 的按模最大和最小的特征值.数值分析模拟卷B填空题（每空2分，共30分）1. 近似数231.0=*x 关于真值229.0=x 有____________位有效数字；2. 设)(x f 可微，求方程)(x f x =根的牛顿迭代格式是_______________________________________________；3. 对1)(3++=x x x f ，差商=]3,2,1,0[f _________________；=]4,3,2,1,0[f ________；4. 已知???? ??-='-=1223,)3,2(A x ，则=∞||||Ax ________________，=)(1A Cond ______________________ ；5. 用二分法求方程01)(3=-+=x x x f 在区间[0,1]的根，进行一步后根所在区间为_________，进行二步后根所在区间为_________________；6. 求解线性方程组=+=+04511532121x x x x 的高斯—赛德尔迭代格式为_______________________________________；该迭代格式迭代矩阵的谱半径=)(G ρ_______________；7. 为使两点数值求积公式：?-+≈111100)()()(x f x f dx x f ωω具有最高的代数精确度，其求积节点应为=0x _____ , =1x _____,==10ωω__________.8. 求积公式)]2()1([23)(30f f dx x f +≈?是否是插值型的__________，其代数精度为___________。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
目录目录 (1)1、题目 (2)2、算法设计方案 (2)3、源程序 (3)4、输出结果 (16)1、题目关于x 、y 、t 、u 、w 的下列方程组0.5cos 2.670.5sin 1.07(1)0.5cos 3.740.5sin 0.79t u v w x t u v w y t u v w x t u v w y +++-=⎧⎪+++-=⎪⇒⎨+++-=⎪⎪+++-=⎩确定了一个二元函数z (,)f x y =。
1、试用数值方法求出(,)f x y 在区域{(,)|00.8,0.54 1.5}D x y x y =≤≤≤≤上的一个近似表达式,0(,)kr s rsr s p x y cx y ==∑要求(,)p x y 一最小的k 值达到以下精度10202700((,)(,))10i j i j i j f x y p x y σ-===-≤∑∑其中0.08,0.50.05i j x i y j ==+。
2、计算****(,),(,)(1,2,...,8,1,2,...,5)i j i j f x y p x y i j ==的值，以观察(,)p x y 逼近(,)f x y 的效果，其中0.1,0.50.2i j x i y j ==+2、算法设计方案根据题目可知，要得到(,)p x y ，必须首先确定(,)f x y 的值，通过一系列(,,)x y z 点进行曲面拟合得到(,)p x y 。
所以之前必须得到对应于x 、y 的z 值。
由已知条件可知，当给定一组(,)x y 的值，方程组（1）就唯一确定，通过对方程组求解就可以解得一组(,,,)Tt u v w ，故而可得到对应于0.08,0.50.05i j x i y j ==+的一系列，(,)t u 值。
然后根据给出的二位数表进行二元函数插值得到对应于(,)x y 的z 值。
1、采用Newton 迭代法进行方程组（1）求解；2、做二元插值时采用分片二次插值，进行分片二次插值需要九个插值节点。
跟定Newton 法解得的(,)t u 值，求得与其最近的插值节点，以该插值节点为中心找八个插值节点即可满足要求。
3、根据求出的一系列z 值，进行曲面拟合，得到近似表达式,0(,)kr s rsr s p x y cx y ==∑4、**(,)i j f x y 的值的求法见1、2，而对于(,)p x y 而言，只需将对应的**(,)i j x y 代入近似表达式中计算即可。
5、在进行曲面拟合时，涉及到矩阵的求逆问题。
由逆矩阵的性质1AA E -=可知，只需求解n 次i i Ax e =，其中i e 表示第（i +1）个元素为1，其余元素为0的 n 维单位坐标向量，=0,1,2...(1)i n -，n 为矩阵A 的阶数，则{i x }组成的向量组即为A 的逆矩阵。
6、求解过程中涉及到的线性方程组求解问题，本次设计采用Gauss 选主元求解方程组。
3、源程序源程序如下： #include "stdafx.h" #include "math.h"#define M 6 #define N 4#define preci 1.0e-12 //Newton 迭代的精度水平； #define Sigma 1.0e-7 //曲面拟合的精度水平； #define X_NUM 11 #define Y_NUM 21FILE *p;double u[M] = {0, 0.4, 0.8, 1.2, 1.6, 2.0}, // u 向量；t[M] = {0, 0.2, 0.4, 0.6, 0.8, 1.0}, // t 向量；ztu[M][M] = { // z 关于t 、u的二维数表；{-0.5, -0.34, 0.14, 0.94, 2.06, 3.5},{-0.42, -0.5, -0.26, 0.3, 1.18, 2.38},{-0.18, -0.5, -0.5, -0.18, 0.46, 1.42},{0.22, -0.34, -0.58, -0.5, -0.1, 0.62},{0.78, -0.02, -0.5, -0.66, -0.5, -0.02},{1.5, 0.46, -0.26, -0.66, -0.74, -0.5},},z[X_NUM][Y_NUM] = {0}, //z的值z_x[X_NUM][Y_NUM] = {0}, //z*的值tuvw[4] = {1,2,1,2}; //Newton迭代法迭代初始向量/************************************************************用Gauss (选主元)方法解带状方程组************************************************************/void Gauss(double F[K][K], double x[K], double _b[K], int rank){int i,j,k,t;int ik;double temp, mik, b[K];for (k = 0; k < rank; k++){b[k] = _b[k];}for(k = 0; k < rank-1; k++){temp = 0;ik = k;//选行号for (i = k; i <= rank; i++){if (fabs(F[i - 1][k - 1]) > fabs(temp)){temp = F[i - 1][k - 1];ik = i;}}//交换if (ik != k){for (j = k; j <= rank; j++){temp = F[k - 1][j - 1];F[k - 1][j - 1] = F[ik - 1][j - 1];F[ik - 1][j - 1] = temp;}temp = b[k];b[k] = b[ik - 1];b[ik] = temp;}for(i = k + 1; i < rank; i++){mik = F[i][k] / F[k][k];for(j = k; j < rank; j++)F[i][j] = F[i][j] - mik * F[k][j];b[i] -= mik * b[k];}}x[rank - 1] = b[rank - 1] / F[rank - 1][rank - 1];for(k = rank - 2; k >= 0; k--){temp = 0;for(j = k+1; j < rank; j++){temp += F[k][j] * x[j];}x[k] = (b[k] - temp) / F[k][k];}}/************************************************************求解矩阵的逆************************************************************/ void inv_N(double A[K][K], int k){int i, j, m, n, t = -1;double temp[K][K] = {0},bak[K][K],e[K] = {0},x[K] = {0},A1[K][K] = {0};for (i = 0; i < k; i++){//复制矩阵A ,以便循环求解方程组for (m = 0; m < k; m++){for (n = 0; n < k; n++){temp[m][n] = A[m][n];}}//设置单位向量if(i == 0)e[i] = 1;else{e[i - 1] = 0;e[i] = 1;}//求解逆矩阵的第i 列向量Gauss(temp, x, e, k);for (j = 0; j < k; j++){A1[j][i] = x[j];}}//改变矩阵A ,使A 存储其逆矩阵for (m = 0; m < k; m++){for (n = 0; n < k; n++){A[m][n] = A1[m][n];}}}/************************************************************求解向量的无穷范数************************************************************/ double Nome(double x[N]){int i;double temp = fabs(x[0]);for (i = 0; i < N; i++){if (fabs(x[i]) > fabs(temp))temp = fabs(x[i]);}return temp;}/************************************************************牛顿迭代法求方程组的解************************************************************/int Newton(double F[K], double F1[K][K], double Answer[K], double x, double y) {int i, j, k;double D[K][K];double rx[K]; //偏差/*复制矩阵F ，使矩阵A 在计算过程中保持不变，使D = A*/for (i = 0;i < N; i++){for (j = 0;j < N; j++){D[i][j] = F1[i][j];}}//给定初始值for (i = 0; i < N; i++){Answer[i] = tuvw[i];}for (k = 0; ; k++) //不限制最大迭代次数{/* 求F（xk）和F1（xk）*/F[0] = -1 * (0.5 * cos(Answer[0]) + Answer[1] + Answer[2] + Answer[3] - x - 2.67);F[1] = -1 * (Answer[0] + 0.5 * sin(Answer[1]) + Answer[2] + Answer[3] - y - 1.07);F[2] = -1 * (0.5 * Answer[0] + Answer[1] + cos(Answer[2]) + Answer[3] - x - 3.74);F[3] = -1 * (Answer[0] + 0.5 * Answer[1] + Answer[2] + sin(Answer[3]) - y - 0.79);F1[0][0] = -0.5 * sin(Answer[0]);F1[0][1] = 1;F1[0][2] = 1;F1[0][3] = 1;F1[1][0] = 1;F1[1][1] = 0.5 * cos(Answer[1]);F1[1][2] = 1;F1[1][3] = 1;F1[2][0] = 0.5;F1[2][1] = 1;F1[2][2] = -1 * sin(Answer[2]);F1[2][3] = 1;F1[3][0] = 1;F1[3][1] = 0.5;F1[3][2] = 1;F1[3][3] = cos(Answer[3]);Gauss(F1, rx, F, 4);for (i = 0; i < N; i++)Answer[i] = Answer[i] + rx[i];if (Nome(rx) / Nome(Answer) < preci)break;}return 1;}/************************************************************二次插值************************************************************/double Qua_Interpolation(double _t, double _u){int r, i, j, k;double I_u[3], I_t[3]; //进行二次插值需要三个点double z;z = 0;r = (int)(fabs((_t / 0.2) + 0.5)); //寻找与t，u 最接近的点，以此为中心找插值节点k = (int)(fabs((_u / 0.4) + 0.5));if(r = 0)r = 1;if(r = 5)r = 4;if(k = 0)k = 1;if(k = 5)k = 4;//关于u 的拉格朗日插值基函数I_u[0] = (_u - u[k]) * (_u - u[k + 1]) / (u[k - 1] - u[k]) / (u[k - 1] - u[k + 1]);I_u[1] = (_u - u[k - 1]) * (_u - u[k + 1]) / (u[k] - u[k - 1]) / (u[k] - u[k + 1]);I_u[2] = (_u - u[k - 1]) * (_u - u[k]) / (u[k + 1] - u[k - 1]) / (u[k + 1] - u[k]);//关于t 的拉格朗日插值基函数I_t[0] = (_t - t[r]) * (_t - t[r + 1]) / (t[r - 1] - t[r]) / (t[r - 1] - t[r + 1]);I_t[1] = (_t - t[r - 1]) * (_t - t[r + 1]) / (t[r] - t[r - 1]) / (t[r] - t[r + 1]);I_t[2] = (_t - t[r - 1]) * (_t - t[r]) / (t[r + 1] - t[r - 1]) / (t[r + 1] - t[r]);for(i = r - 1; i <= r + 1; i++){for(j = k - 1; j <= k + 1; j++){z += I_t[i - r + 1] * I_u[j - k + 1] * ztu[i][j];}}return z;}/************************************************************曲面拟合************************************************************/void Sur_fitting(double A[K][K], int k_now){int i, j, k, m, n, k_max, k_inv;double x[X_NUM],y[Y_NUM],BTB[K][K] = {0},GTG[K][K] = {0},BTBBT[K][X_NUM],BTBBTU[K][Y_NUM],temp, sigma;for(i = 0; i < X_NUM; i++){x[i] = 0.08 * i;}for(j = 0; j < Y_NUM; j++){y[j] = 0.5 + 0.05 * j;}for(i = 0; i < k_now; i++){for(j = 0; j < k_now; j++){temp = 0;for(k = 0; k < X_NUM; k++)temp += pow(x[k], i) * pow(x[k], j);BTB[i][j] = temp;}}inv_N(BTB, k_now);for(i = 0; i < k_now; i++){for(j = 0; j < k_now; j++){temp = 0;for(k = 0; k < Y_NUM; k++)temp+= pow(y[k], i) * pow(y[k], j);GTG[i][j] = temp;}}inv_N(GTG, k_now);for(i = 0; i < k_now; i++){for(j = 0; j < X_NUM; j++)temp = 0;for(k = 0; k < k_now; k++)temp += BTB[i][k] * pow(x[j],k);BTBBT[i][j] = temp;}}for(i = 0; i < k_now; i++){for(j = 0; j < Y_NUM; j++){temp = 0;for(k = 0; k < X_NUM; k++)temp += BTBBT[i][k] * z[k][j];BTBBTU[i][j] = temp;}}for(i = 0; i < k_now; i++){for(j = 0; j < k_now; j++){temp = 0;for(k = 0; k < Y_NUM; k++)temp += BTBBTU[i][k] * pow(y[k],j);BTB[i][j] = temp;}}for(i = 0; i < k_now; i++){for(j = 0; j < k_now; j++){temp = 0;for(k = 0; k < k_now; k++)temp += BTB[i][k] * GTG[k][j];A[i][j] = temp;}}/************************************************************主函数************************************************************/int main(int argc, char* argv[]){double temp, temp1;double F[K], F1[K][K];/*t, u, v, w的向量*/double Answer[K], Answer_x[K];double Crs[K][K];int i, j, k, k_inv;int m, n;p = fopen("e:\\Final_Result.txt", "w");fprintf(p, "1、插值求z 如下：\n\n");fclose(p);for (i = 0; i <= 10; i++){for (j = 0; j <= 20; j++){Newton(F, F1, Answer, 0.08 * i, 0.5 + 0.05 * j);z[i][j] = Qua_Interpolation(Answer[0], Answer[1]);p = fopen("e:\\Final_Result.txt", "at+");fprintf(p, " x%-2d=%f, y%-2d=%f, t=%f, u=%f, z=f(x%-2d,y%-2d)=%13.12E\n", i, 0.08 * i, j, 0.5 + 0.05 * j, Answer[0], Answer[1], i, j, z[i][j]);fclose(p);}}p = fopen("e:\\Final_Result.txt", "at+");fprintf(p, "\n2、曲面拟合结果如下：\n\n");fclose(p);/*不限制次数*/for (k = 1; ; k++){temp = 0;temp1 = 0;k_inv = k;printf("k = %d: \n", k);Sur_fitting(Crs, k + 1);temp = 0;for(i = 0; i < X_NUM; i++){for(j = 0; j < Y_NUM; j++){temp1 = 0;for(m = 0; m < k + 1; m++){for(n = 0; n < k + 1; n++){temp1 += Crs[m][n] * pow(0.08 * i, m) * pow(0.5 + 0.05 * j, n);}}temp += (temp1 - z[i][j]) * (temp1 - z[i][j]);}}p = fopen("e:\\Final_Result.txt", "at+");fprintf(p, " k = %d σ = %13.12E\n\n", k, temp);fclose(p);if (temp < Sigma){p = fopen("e:\\Final_Result.txt", "at+");fprintf(p,"系数矩阵:\n Crs[%d][%d] = \n {\n", k + 1, k + 1);for(i = 0; i < k + 1; i++){fprintf(p," {");for(j = 0; j < k; j++)fprintf(p,"%.12lE ", Crs[i][j]);fprintf(p,"%.12lE", Crs[i][k]);fprintf(p,"},\n");}fprintf(p," }\n\n拟合最小k 值已求出：");fprintf(p,"k_min = %d\n\n", k_inv);fclose(p);break;}}p = fopen("e:\\Final_Result.txt", "at+");fprintf(p, "3、f(x*, y*)和p(x*, y*)结果如下:\n\n");fclose(p);for (i = 1; i <= 8; i++){for (j = 1; j <= 5; j++){Newton(F, F1, Answer_x, 0.1 * i, 0.5 + 0.2 * j);z_x[i][j] = Qua_Interpolation(Answer_x[0], Answer_x[1]);//求p(x*, y*)temp1 = 0;for(m = 0; m <= k_inv; m++){for(n = 0; n <= k_inv; n++){temp1 += Crs[m][n] * pow(0.1 * i, m) * pow(0.5 + 0.2 * j, n);}}p = fopen("e:\\Final_Result.txt", "at+");fprintf(p, " x*%d=%f, y*%d=%f, f(x*%d,y*%d)=%13.12E, p(x*%d,y*%d)=%13.12E\n", i, 0.1 * i, j, 0.5 + 0.2 * j, i, j, z_x[i][j], i, j, temp1);fclose(p);}}printf("结果已求出，详见E:\\Final_Result.txt!\n\n");return 0;}4、输出结果1、插值求z 如下：x0 =0.000000, y0 =0.500000, t=0.243185, u=1.345231, z=f(x0 ,y0 )=4.465040184799E-001 x0 =0.000000, y1 =0.550000, t=0.269004, u=1.321533, z=f(x0 ,y1 )=3.246832629274E-001 x0 =0.000000, y2 =0.600000, t=0.295402, u=1.298659, z=f(x0 ,y2 )=2.101596866825E-001 x0 =0.000000, y3 =0.650000, t=0.322314, u=1.276576, z=f(x0 ,y3 )=1.030436083159E-001 x0 =0.000000, y4 =0.700000, t=0.349678, u=1.255250, z=f(x0 ,y4 )=3.401895562659E-003 x0 =0.000000, y5 =0.750000, t=0.377443, u=1.234655, z=f(x0 ,y5 )=-8.873581363801E-002 x0 =0.000000, y6 =0.800000, t=0.405561, u=1.214764, z=f(x0 ,y6 )=-1.733716327497E-001 x0 =0.000000, y7 =0.850000, t=0.433991, u=1.195553, z=f(x0 ,y7 )=-2.505346114666E-001 x0 =0.000000, y8 =0.900000, t=0.462698, u=1.177002, z=f(x0 ,y8 )=-3.202765063876E-001 x0 =0.000000, y9 =0.950000, t=0.491648, u=1.159089, z=f(x0 ,y9 )=-3.826680661097E-001 x0 =0.000000, y10=1.000000, t=0.520814, u=1.141798, z=f(x0 ,y10)=-4.377957667384E-001 x0 =0.000000, y11=1.050000, t=0.550170, u=1.125109, z=f(x0 ,y11)=-4.857589414438E-001 x0 =0.000000, y12=1.100000, t=0.579694, u=1.109008, z=f(x0 ,y12)=-5.266672548835E-001 x0 =0.000000, y13=1.150000, t=0.609367, u=1.093477, z=f(x0 ,y13)=-5.606384797965E-001 x0 =0.000000, y14=1.200000, t=0.639172, u=1.078503, z=f(x0 ,y14)=-5.877965387677E-001 x0 =0.000000, y15=1.250000, t=0.669094, u=1.064071, z=f(x0 ,y15)=-6.0826********E-001 x0 =0.000000, y16=1.300000, t=0.699119, u=1.050167, z=f(x0 ,y16)=-6.221894528764E-001 x0 =0.000000, y17=1.350000, t=0.729235, u=1.036777, z=f(x0 ,y17)=-6.296883781856E-001 x0 =0.000000, y18=1.400000, t=0.759434, u=1.023889, z=f(x0 ,y18)=-6.308997600028E-001 x0 =0.000000, y19=1.450000, t=0.789706, u=1.011490, z=f(x0 ,y19)=-6.259561525454E-001 x0 =0.000000, y20=1.500000, t=0.820043, u=0.999567, z=f(x0 ,y20)=-6.149885466094E-001 x1 =0.080000, y0 =0.500000, t=0.228064, u=1.414952, z=f(x1 ,y0 )=6.380152265102E-001 x1 =0.080000, y1 =0.550000, t=0.253304, u=1.391219, z=f(x1 ,y1 )=5.0661********E-001 x1 =0.080000, y2 =0.600000, t=0.279170, u=1.368312, z=f(x1 ,y2 )=3.821763692772E-001 x1 =0.080000, y3 =0.650000, t=0.305590, u=1.346197, z=f(x1 ,y3 )=2.648634911536E-001 x1 =0.080000, y4 =0.700000, t=0.332503, u=1.324841, z=f(x1 ,y4 )=1.547802002848E-001 x1 =0.080000, y5 =0.750000, t=0.359851, u=1.304216, z=f(x1 ,y5 )=5.199268349093E-002 x1 =0.080000, y6 =0.800000, t=0.387587, u=1.284294, z=f(x1 ,y6 )=-4.346804020491E-002x1 =0.080000, y7 =0.850000, t=0.415664, u=1.265052, z=f(x1 ,y7 )=-1.316010567885E-001 x1 =0.080000, y8 =0.900000, t=0.444046, u=1.246468, z=f(x1 ,y8 )=-2.124310883088E-001 x1 =0.080000, y9 =0.950000, t=0.472697, u=1.228521, z=f(x1 ,y9 )=-2.860045510580E-001 x1 =0.080000, y10=1.000000, t=0.501586, u=1.211194, z=f(x1 ,y10)=-3.523860789794E-001 x1 =0.080000, y11=1.050000, t=0.530687, u=1.194468, z=f(x1 ,y11)=-4.116554565222E-001 x1 =0.080000, y12=1.100000, t=0.559977, u=1.178326, z=f(x1 ,y12)=-4.639049115188E-001 x1 =0.080000, y13=1.150000, t=0.589432, u=1.162754, z=f(x1 ,y13)=-5.092367247005E-001 x1 =0.080000, y14=1.200000, t=0.619036, u=1.147736, z=f(x1 ,y14)=-5.477611179623E-001 x1 =0.080000, y15=1.250000, t=0.648772, u=1.133257, z=f(x1 ,y15)=-5.795943883391E-001 x1 =0.080000, y16=1.300000, t=0.678624, u=1.119304, z=f(x1 ,y16)=-6.048572588895E-001 x1 =0.080000, y17=1.350000, t=0.708581, u=1.105864, z=f(x1 ,y17)=-6.236734213318E-001 x1 =0.080000, y18=1.400000, t=0.738631, u=1.092923, z=f(x1 ,y18)=-6.361682484133E-001 x1 =0.080000, y19=1.450000, t=0.768765, u=1.080469, z=f(x1 ,y19)=-6.424676566901E-001 x1 =0.080000, y20=1.500000, t=0.798973, u=1.068489, z=f(x1 ,y20)=-6.426971026996E-001 x2 =0.160000, y0 =0.500000, t=0.213867, u=1.483347, z=f(x2 ,y0 )=8.400813957651E-001 x2 =0.160000, y1 =0.550000, t=0.238510, u=1.459575, z=f(x2 ,y1 )=6.997641656726E-001 x2 =0.160000, y2 =0.600000, t=0.263823, u=1.436630, z=f(x2 ,y2 )=5.660614423514E-001 x2 =0.160000, y3 =0.650000, t=0.289734, u=1.414478, z=f(x2 ,y3 )=4.391716081175E-001 x2 =0.160000, y4 =0.700000, t=0.316175, u=1.393085, z=f(x2 ,y4 )=3.192421380408E-001 x2 =0.160000, y5 =0.750000, t=0.343088, u=1.372422, z=f(x2 ,y5 )=2.063761923874E-001 x2 =0.160000, y6 =0.800000, t=0.370420, u=1.352463, z=f(x2 ,y6 )=1.006385238914E-001 x2 =0.160000, y7 =0.850000, t=0.398126, u=1.333182, z=f(x2 ,y7 )=2.060740067836E-003 x2 =0.160000, y8 =0.900000, t=0.426163, u=1.314558, z=f(x2 ,y8 )=-8.935402476698E-002 x2 =0.160000, y9 =0.950000, t=0.454495, u=1.296571, z=f(x2 ,y9 )=-1.736269688649E-001 x2 =0.160000, y10=1.000000, t=0.483090, u=1.279201, z=f(x2 ,y10)=-2.507999561599E-001 x2 =0.160000, y11=1.050000, t=0.511918, u=1.262430, z=f(x2 ,y11)=-3.209322694446E-001 x2 =0.160000, y12=1.100000, t=0.540954, u=1.246242, z=f(x2 ,y12)=-3.840977350046E-001 x2 =0.160000, y13=1.150000, t=0.570176, u=1.230621, z=f(x2 ,y13)=-4.403821754175E-001 x2 =0.160000, y14=1.200000, t=0.599562, u=1.215552, z=f(x2 ,y14)=-4.898811523126E-001 x2 =0.160000, y15=1.250000, t=0.629096, u=1.201021, z=f(x2 ,y15)=-5.326979655338E-001 x2 =0.160000, y16=1.300000, t=0.658761, u=1.187013, z=f(x2 ,y16)=-5.689418792921E-001 x2 =0.160000, y17=1.350000, t=0.688544, u=1.173516, z=f(x2 ,y17)=-5.987265495151E-001 x2 =0.160000, y18=1.400000, t=0.718432, u=1.160517, z=f(x2 ,y18)=-6.221686297503E-001 x2 =0.160000, y19=1.450000, t=0.748414, u=1.148002, z=f(x2 ,y19)=-6.393865356972E-001 x2 =0.160000, y20=1.500000, t=0.778481, u=1.135960, z=f(x2 ,y20)=-6.504993507878E-001 x3 =0.240000, y0 =0.500000, t=0.200585, u=1.550507, z=f(x3 ,y0 )=1.0515********E+000 x3 =0.240000, y1 =0.550000, t=0.224618, u=1.526693, z=f(x3 ,y1 )=9.029*********E-001 x3 =0.240000, y2 =0.600000, t=0.249366, u=1.503707, z=f(x3 ,y2 )=7.605802668593E-001x3 =0.240000, y3 =0.650000, t=0.274751, u=1.481514, z=f(x3 ,y3 )=6.247151981455E-001 x3 =0.240000, y4 =0.700000, t=0.300705, u=1.460080, z=f(x3 ,y4 )=4.955197560009E-001 x3 =0.240000, y5 =0.750000, t=0.327165, u=1.439375, z=f(x3 ,y5 )=3.731340427746E-001 x3 =0.240000, y6 =0.800000, t=0.354078, u=1.419373, z=f(x3 ,y6 )=2.576567488723E-001 x3 =0.240000, y7 =0.850000, t=0.381394, u=1.400049, z=f(x3 ,y7 )=1.491505594102E-001 x3 =0.240000, y8 =0.900000, t=0.409069, u=1.381381, z=f(x3 ,y8 )=4.764698677336E-002 x3 =0.240000, y9 =0.950000, t=0.437066, u=1.363347, z=f(x3 ,y9 )=-4.684932320146E-002 x3 =0.240000, y10=1.000000, t=0.465349, u=1.345928, z=f(x3 ,y10)=-1.343567603849E-001 x3 =0.240000, y11=1.050000, t=0.493888, u=1.329107, z=f(x3 ,y11)=-2.149133449274E-001 x3 =0.240000, y12=1.100000, t=0.522655, u=1.312867, z=f(x3 ,y12)=-2.885737006348E-001 x3 =0.240000, y13=1.150000, t=0.551627, u=1.297192, z=f(x3 ,y13)=-3.554063647857E-001 x3 =0.240000, y14=1.200000, t=0.580781, u=1.282067, z=f(x3 ,y14)=-4.154913964886E-001 x3 =0.240000, y15=1.250000, t=0.610098, u=1.267477, z=f(x3 ,y15)=-4.689182499695E-001 x3 =0.240000, y16=1.300000, t=0.639561, u=1.253409, z=f(x3 ,y16)=-5.157838831247E-001 x3 =0.240000, y17=1.350000, t=0.669156, u=1.239851, z=f(x3 ,y17)=-5.561910752001E-001 x3 =0.240000, y18=1.400000, t=0.698868, u=1.226787, z=f(x3 ,y18)=-5.902469305629E-001 x3 =0.240000, y19=1.450000, t=0.728686, u=1.214207, z=f(x3 ,y19)=-6.180615482412E-001 x3 =0.240000, y20=1.500000, t=0.758600, u=1.202098, z=f(x3 ,y20)=-6.397468392579E-001 x4 =0.320000, y0 =0.500000, t=0.188204, u=1.616511, z=f(x4 ,y0 )=1.271246751481E+000 x4 =0.320000, y1 =0.550000, t=0.211620, u=1.592655, z=f(x4 ,y1 )=1.115002018146E+000 x4 =0.320000, y2 =0.600000, t=0.235792, u=1.569627, z=f(x4 ,y2 )=9.646077272154E-001 x4 =0.320000, y3 =0.650000, t=0.260641, u=1.547391, z=f(x4 ,y3 )=8.203473694749E-001 x4 =0.320000, y4 =0.700000, t=0.286094, u=1.525914, z=f(x4 ,y4 )=6.824476781794E-001 x4 =0.320000, y5 =0.750000, t=0.312089, u=1.505166, z=f(x4 ,y5 )=5.510852085975E-001 x4 =0.320000, y6 =0.800000, t=0.338568, u=1.485118, z=f(x4 ,y6 )=4.263923859018E-001 x4 =0.320000, y7 =0.850000, t=0.365480, u=1.465747, z=f(x4 ,y7 )=3.084629956332E-001 x4 =0.320000, y8 =0.900000, t=0.392779, u=1.447030, z=f(x4 ,y8 )=1.973571296919E-001 x4 =0.320000, y9 =0.950000, t=0.420426, u=1.428946, z=f(x4 ,y9 )=9.310562085941E-002 x4 =0.320000, y10=1.000000, t=0.448382, u=1.411475, z=f(x4 ,y10)=-4.285992234033E-003 x4 =0.320000, y11=1.050000, t=0.476617, u=1.394600, z=f(x4 ,y11)=-9.483392529689E-002 x4 =0.320000, y12=1.100000, t=0.505101, u=1.378303, z=f(x4 ,y12)=-1.785729903640E-001 x4 =0.320000, y13=1.150000, t=0.533808, u=1.362570, z=f(x4 ,y13)=-2.555537790546E-001 x4 =0.320000, y14=1.200000, t=0.562715, u=1.347384, z=f(x4 ,y14)=-3.258401501575E-001 x4 =0.320000, y15=1.250000, t=0.591802, u=1.332732, z=f(x4 ,y15)=-3.895069883634E-001 x4 =0.320000, y16=1.300000, t=0.621051, u=1.318600, z=f(x4 ,y16)=-4.466382045995E-001 x4 =0.320000, y17=1.350000, t=0.650444, u=1.304975, z=f(x4 ,y17)=-4.973249517677E-001 x4 =0.320000, y18=1.400000, t=0.679968, u=1.291844, z=f(x4 ,y18)=-5.416640326994E-001 x4 =0.320000, y19=1.450000, t=0.709610, u=1.279194, z=f(x4 ,y19)=-5.797564797951E-001x4 =0.320000, y20=1.500000, t=0.739359, u=1.267014, z=f(x4 ,y20)=-6.117062881476E-001 x5 =0.400000, y0 =0.500000, t=0.176701, u=1.681430, z=f(x5 ,y0 )=1.498321052481E+000 x5 =0.400000, y1 =0.550000, t=0.199497, u=1.657534, z=f(x5 ,y1 )=1.334998632066E+000 x5 =0.400000, y2 =0.600000, t=0.223088, u=1.634465, z=f(x5 ,y2 )=1.177125123739E+000 x5 =0.400000, y3 =0.650000, t=0.247393, u=1.612188, z=f(x5 ,y3 )=1.025*********E+000 x5 =0.400000, y4 =0.700000, t=0.272338, u=1.590667, z=f(x5 ,y4 )=8.789600231743E-001 x5 =0.400000, y5 =0.750000, t=0.297858, u=1.569874, z=f(x5 ,y5 )=7.391451087035E-001 x5 =0.400000, y6 =0.800000, t=0.323892, u=1.549780, z=f(x5 ,y6 )=6.0574********E-001 x5 =0.400000, y7 =0.850000, t=0.350389, u=1.530361, z=f(x5 ,y7 )=4.788838610666E-001 x5 =0.400000, y8 =0.900000, t=0.377301, u=1.511593, z=f(x5 ,y8 )=3.586506258818E-001 x5 =0.400000, y9 =0.950000, t=0.404584, u=1.493457, z=f(x5 ,y9 )=2.451022361964E-001 x5 =0.400000, y10=1.000000, t=0.432202, u=1.475931, z=f(x5 ,y10)=1.382683509285E-001 x5 =0.400000, y11=1.050000, t=0.460120, u=1.458999, z=f(x5 ,y11)=3.815486540699E-002 x5 =0.400000, y12=1.100000, t=0.488308, u=1.442644, z=f(x5 ,y12)=-5.525282116814E-002 x5 =0.400000, y13=1.150000, t=0.516738, u=1.426849, z=f(x5 ,y13)=-1.419868808137E-001 x5 =0.400000, y14=1.200000, t=0.545386, u=1.411600, z=f(x5 ,y14)=-2.220944390959E-001 x5 =0.400000, y15=1.250000, t=0.574230, u=1.396883, z=f(x5 ,y15)=-2.956352324598E-001 x5 =0.400000, y16=1.300000, t=0.603252, u=1.382683, z=f(x5 ,y16)=-3.626795115028E-001 x5 =0.400000, y17=1.350000, t=0.632432, u=1.368989, z=f(x5 ,y17)=-4.233061642240E-001 x5 =0.400000, y18=1.400000, t=0.661756, u=1.355786, z=f(x5 ,y18)=-4.776010361325E-001 x5 =0.400000, y19=1.450000, t=0.691211, u=1.343064, z=f(x5 ,y19)=-5.256554266672E-001 x5 =0.400000, y20=1.500000, t=0.720783, u=1.330809, z=f(x5 ,y20)=-5.675647436551E-001 x6 =0.480000, y0 =0.500000, t=0.166048, u=1.745329, z=f(x6 ,y0 )=1.731892740382E+000 x6 =0.480000, y1 =0.550000, t=0.188225, u=1.721396, z=f(x6 ,y1 )=1.562034577208E+000 x6 =0.480000, y2 =0.600000, t=0.211234, u=1.698290, z=f(x6 ,y2 )=1.397216918208E+000 x6 =0.480000, y3 =0.650000, t=0.234993, u=1.675972, z=f(x6 ,y3 )=1.237801006739E+000 x6 =0.480000, y4 =0.700000, t=0.259425, u=1.654410, z=f(x6 ,y4 )=1.0840********E+000 x6 =0.480000, y5 =0.750000, t=0.284463, u=1.633572, z=f(x6 ,y5 )=9.363227723149E-001 x6 =0.480000, y6 =0.800000, t=0.310047, u=1.613433, z=f(x6 ,y6 )=7.947044490537E-001 x6 =0.480000, y7 =0.850000, t=0.336120, u=1.593965, z=f(x6 ,y7 )=6.593871980282E-001 x6 =0.480000, y8 =0.900000, t=0.362635, u=1.575147, z=f(x6 ,y8 )=5.304875868399E-001 x6 =0.480000, y9 =0.950000, t=0.389546, u=1.556958, z=f(x6 ,y9 )=4.080886854542E-001 x6 =0.480000, y10=1.000000, t=0.416815, u=1.539377, z=f(x6 ,y10)=2.922442012295E-001 x6 =0.480000, y11=1.050000, t=0.444406, u=1.522387, z=f(x6 ,y11)=1.829822068535E-001 x6 =0.480000, y12=1.100000, t=0.472287, u=1.505971, z=f(x6 ,y12)=8.030849403543E-002 x6 =0.480000, y13=1.150000, t=0.500429, u=1.490114, z=f(x6 ,y13)=-1.579041305164E-002 x6 =0.480000, y14=1.200000, t=0.528807, u=1.474800, z=f(x6 ,y14)=-1.0534********E-001 x6 =0.480000, y15=1.250000, t=0.557397, u=1.460015, z=f(x6 ,y15)=-1.883980906096E-001x6 =0.480000, y16=1.300000, t=0.586180, u=1.445746, z=f(x6 ,y16)=-2.650071493189E-001 x6 =0.480000, y17=1.350000, t=0.615137, u=1.431980, z=f(x6 ,y17)=-3.352378389040E-001 x6 =0.480000, y18=1.400000, t=0.644251, u=1.418704, z=f(x6 ,y18)=-3.991645038868E-001 x6 =0.480000, y19=1.450000, t=0.673508, u=1.405907, z=f(x6 ,y19)=-4.568681433016E-001 x6 =0.480000, y20=1.500000, t=0.702894, u=1.393575, z=f(x6 ,y20)=-5.084349932782E-001 x7 =0.560000, y0 =0.500000, t=0.156210, u=1.808266, z=f(x7 ,y0 )=1.971221786400E+000 x7 =0.560000, y1 =0.550000, t=0.177775, u=1.784302, z=f(x7 ,y1 )=1.795329599501E+000 x7 =0.560000, y2 =0.600000, t=0.200206, u=1.761161, z=f(x7 ,y2 )=1.624067113228E+000 x7 =0.560000, y3 =0.650000, t=0.223420, u=1.738807, z=f(x7 ,y3 )=1.457830582708E+000 x7 =0.560000, y4 =0.700000, t=0.247338, u=1.717205, z=f(x7 ,y4 )=1.296954649752E+000 x7 =0.560000, y5 =0.750000, t=0.271893, u=1.696326, z=f(x7 ,y5 )=1.141718105447E+000 x7 =0.560000, y6 =0.800000, t=0.297020, u=1.676142, z=f(x7 ,y6 )=9.923495333243E-001 x7 =0.560000, y7 =0.850000, t=0.322665, u=1.656628, z=f(x7 ,y7 )=8.490326633294E-001 x7 =0.560000, y8 =0.900000, t=0.348777, u=1.637760, z=f(x7 ,y8 )=7.119113522641E-001 x7 =0.560000, y9 =0.950000, t=0.375309, u=1.619518, z=f(x7 ,y9 )=5.810941589219E-001 x7 =0.560000, y10=1.000000, t=0.402221, u=1.601883, z=f(x7 ,y10)=4.566585132334E-001 x7 =0.560000, y11=1.050000, t=0.429477, u=1.584835, z=f(x7 ,y11)=3.386544961394E-001 x7 =0.560000, y12=1.100000, t=0.457042, u=1.568359, z=f(x7 ,y12)=2.271082557696E-001 x7 =0.560000, y13=1.150000, t=0.484887, u=1.552439, z=f(x7 ,y13)=1.220250891932E-001 x7 =0.560000, y14=1.200000, t=0.512986, u=1.537060, z=f(x7 ,y14)=2.339221963760E-002 x7 =0.560000, y15=1.250000, t=0.541313, u=1.522207, z=f(x7 ,y15)=-6.881870197104E-002 x7 =0.560000, y16=1.300000, t=0.569849, u=1.507868, z=f(x7 ,y16)=-1.546493442129E-001 x7 =0.560000, y17=1.350000, t=0.598573, u=1.494030, z=f(x7 ,y17)=-2.341526664587E-001 x7 =0.560000, y18=1.400000, t=0.627467, u=1.480679, z=f(x7 ,y18)=-3.0739********E-001 x7 =0.560000, y19=1.450000, t=0.656517, u=1.467806, z=f(x7 ,y19)=-3.744348623481E-001 x7 =0.560000, y20=1.500000, t=0.685708, u=1.455396, z=f(x7 ,y20)=-4.353605565359E-001 x8 =0.640000, y0 =0.500000, t=0.147151, u=1.870297, z=f(x8 ,y0 )=2.215667863688E+000 x8 =0.640000, y1 =0.550000, t=0.168114, u=1.846307, z=f(x8 ,y1 )=2.034201133607E+000 x8 =0.640000, y2 =0.600000, t=0.189975, u=1.823136, z=f(x8 ,y2 )=1.856955143619E+000 x8 =0.640000, y3 =0.650000, t=0.212648, u=1.800749, z=f(x8 ,y3 )=1.684358164161E+000 x8 =0.640000, y4 =0.700000, t=0.236055, u=1.779112, z=f(x8 ,y4 )=1.516776352400E+000 x8 =0.640000, y5 =0.750000, t=0.260127, u=1.758194, z=f(x8 ,y5 )=1.354519041151E+000 x8 =0.640000, y6 =0.800000, t=0.284798, u=1.737969, z=f(x8 ,y6 )=1.197844086673E+000 x8 =0.640000, y7 =0.850000, t=0.310012, u=1.718410, z=f(x8 ,y7 )=1.046963049419E+000 x8 =0.640000, y8 =0.900000, t=0.335717, u=1.699495, z=f(x8 ,y8 )=9.020*********E-001 x8 =0.640000, y9 =0.950000, t=0.361866, u=1.681203, z=f(x8 ,y9 )=7.632264776629E-001 x8 =0.640000, y10=1.000000, t=0.388416, u=1.663514, z=f(x8 ,y10)=6.306048219543E-001 x8 =0.640000, y11=1.050000, t=0.415330, u=1.646410, z=f(x8 ,y11)=5.042528145972E-001x8 =0.640000, y12=1.100000, t=0.442573, u=1.629875, z=f(x8 ,y12)=3.842167155457E-001 x8 =0.640000, y13=1.150000, t=0.470115, u=1.613892, z=f(x8 ,y13)=2.705204766410E-001 x8 =0.640000, y14=1.200000, t=0.497927, u=1.598448, z=f(x8 ,y14)=1.631685723996E-001 x8 =0.640000, y15=1.250000, t=0.525984, u=1.583528, z=f(x8 ,y15)=6.214855811676E-002 x8 =0.640000, y16=1.300000, t=0.554264, u=1.569119, z=f(x8 ,y16)=-3.256661939682E-002 x8 =0.640000, y17=1.350000, t=0.582747, u=1.555209, z=f(x8 ,y17)=-1.210165348444E-001 x8 =0.640000, y18=1.400000, t=0.611414, u=1.541784, z=f(x8 ,y18)=-2.032513996228E-001 x8 =0.640000, y19=1.450000, t=0.640248, u=1.528834, z=f(x8 ,y19)=-2.793303595584E-001 x8 =0.640000, y20=1.500000, t=0.669236, u=1.516347, z=f(x8 ,y20)=-3.493199575400E-001 x9 =0.720000, y0 =0.500000, t=0.138831, u=1.931471, z=f(x9 ,y0 )=2.464684222660E+000 x9 =0.720000, y1 =0.550000, t=0.159207, u=1.907461, z=f(x9 ,y1 )=2.278058979399E+000 x9 =0.720000, y2 =0.600000, t=0.180508, u=1.884266, z=f(x9 ,y2 )=2.0952********E+000 x9 =0.720000, y3 =0.650000, t=0.202648, u=1.861851, z=f(x9 ,y3 )=1.916718127997E+000 x9 =0.720000, y4 =0.700000, t=0.225550, u=1.840182, z=f(x9 ,y4 )=1.742854628776E+000 x9 =0.720000, y5 =0.750000, t=0.249143, u=1.819230, z=f(x9 ,y5 )=1.573998427334E+000 x9 =0.720000, y6 =0.800000, t=0.273360, u=1.798966, z=f(x9 ,y6 )=1.410434835231E+000 x9 =0.720000, y7 =0.850000, t=0.298144, u=1.779367, z=f(x9 ,y7 )=1.252401750608E+000 x9 =0.720000, y8 =0.900000, t=0.323442, u=1.760407, z=f(x9 ,y8 )=1.100094409628E+000 x9 =0.720000, y9 =0.950000, t=0.349205, u=1.742067, z=f(x9 ,y9 )=9.536698512613E-001 x9 =0.720000, y10=1.000000, t=0.375391, u=1.724327, z=f(x9 ,y10)=8.132510552489E-001 x9 =0.720000, y11=1.050000, t=0.401960, u=1.707169, z=f(x9 ,y11)=6.789307429659E-001 x9 =0.720000, y12=1.100000, t=0.428876, u=1.690577, z=f(x9 ,y12)=5.507748485043E-001 x9 =0.720000, y13=1.150000, t=0.456109, u=1.674535, z=f(x9 ,y13)=4.288256769731E-001 x9 =0.720000, y14=1.200000, t=0.483629, u=1.659027, z=f(x9 ,y14)=3.131047717398E-001 x9 =0.720000, y15=1.250000, t=0.511410, u=1.644042, z=f(x9 ,y15)=2.036155140327E-001 x9 =0.720000, y16=1.300000, t=0.539429, u=1.629564, z=f(x9 ,y16)=1.003454782409E-001 x9 =0.720000, y17=1.350000, t=0.567664, u=1.615583, z=f(x9 ,y17)=3.268565186572E-003 x9 =0.720000, y18=1.400000, t=0.596096, u=1.602085, z=f(x9 ,y18)=-8.765306591329E-002 x9 =0.720000, y19=1.450000, t=0.624709, u=1.589059, z=f(x9 ,y19)=-1.724672478188E-001 x9 =0.720000, y20=1.500000, t=0.653486, u=1.576494, z=f(x9 ,y20)=-2.512302207523E-001 x10=0.800000, y0 =0.500000, t=0.131209, u=1.991838, z=f(x10,y0 )=2.717811109469E+000 x10=0.800000, y1 =0.550000, t=0.151014, u=1.967813, z=f(x10,y1 )=2.526399501256E+000 x10=0.800000, y2 =0.600000, t=0.171769, u=1.944600, z=f(x10,y2 )=2.338411386860E+000 x10=0.800000, y3 =0.650000, t=0.193389, u=1.922161, z=f(x10,y3 )=2.154329377280E+000 x10=0.800000, y4 =0.700000, t=0.215795, u=1.900466, z=f(x10,y4 )=1.974574556652E+000 x10=0.800000, y5 =0.750000, t=0.238915, u=1.879483, z=f(x10,y5 )=1.799510579099E+000 x10=0.800000, y6 =0.800000, t=0.262683, u=1.859186, z=f(x10,y6 )=1.629448220554E+000 x10=0.800000, y7 =0.850000, t=0.287040, u=1.839549, z=f(x10,y7 )=1.464650043751E+000。