下载文档原格式
数值分析第三次作业1.设计方案对Fredholm积分方程，用迭代法进行求解：()'(())u x A u x=，其中11(())()(,)()A u x g x K x y u y dy-=-⋅⎰对于公式中的积分部分用数值积分方法。
复化梯形积分法，取2601个节点，取迭代次数上限为50次。
实际计算迭代次数为18次，最后算得误差为r= 0.97E-10。
复化Simpson积分法，取迭代次数上限为50次，取2*41+1，即83个节点时能满足精度要求。
实际计算迭代次数为17次，最后的误差为r= 0.97E-10。
Guass积分法选择的Gauss—Legendre法，取迭代次数上限为50次，直接选择8个节点，满足精度要求。
实际计算迭代次数为24次，最后算得误差为r= 0.87E-10。
2.全部源程序module integralimplicit nonecontains!//////////复化梯形subroutine trapezoid(m)implicit noneinteger :: i,j,k,mreal*8 :: x(m+1),u(m+1)real*8 :: sum,sum1,g,r,hreal*8 :: e=1.0e-10h=2./mdo i=1,m+1x(i)=-1.+(i-1)*hend dou=0.02do k=1,50do i=1,m+1sum1=0.g=dexp(x(i)*4.)+(dexp(x(i)+4.)-dexp(-4.-x(i)))/(x(i)+4.)do j=2,msum1=sum1+dexp(x(i)*x(j))*u(j)end dosum=h/2.*(dexp(x(i)*-1.)*u(1)+dexp(x(i)*1.)*u(m+1)+2*sum1)u(i)=g-sumend dor=h/2.*((dexp(x(1)*4)-u(1))**2+(dexp(x(m+1)*4)-u(m+1))**2) do i=2,mr=r+h*(dexp(x(i)*4)-u(i))**2end doif(dabs(r)<=e) exitend dowrite(*,*) kopen(1,file="trapezoid.txt")do i=1,m+1write(1,'(3(f18.12))') x(i),u(i),dexp(x(i)*4.)end dowrite(1,'(4x,a2,e9.2)') "r=",rclose(1)returnend subroutine trapezoid!///////////复化simpsonsubroutine simpson(m)implicit noneinteger :: i,j,k,mreal*8 :: x(2*m+1),u(2*m+1)real*8 :: sum,sum1,sum2,g,r,hreal*8 :: e=1.0e-10h=2./(2.*m)do i=1,2*m+1x(i)=-1.+(i-1)*hend dou=0.02do k=1,50do i=1,2*m+1sum1=0.sum2=0.g=dexp(x(i)*4.)+(dexp(x(i)+4.)-dexp(-4.-x(i)))/(x(i)+4.)do j=1,msum1=sum1+dexp(x(i)*x(2*j))*u(2*j)end dodo j=1,m-1sum2=sum2+dexp(x(i)*x(2*j+1))*u(2*j+1)sum=h/3.*(dexp(x(i)*-1.)*u(1)+dexp(x(i)*1.)*u(2*m+1)+4*sum1+2*sum2) u(i)=g-sumend dor=h/3.*((dexp(x(1)*4)-u(1))**2+(dexp(x(2*m+1)*4)-u(2*m+1))**2)do i=1,mr=r+4.*h/3.*(dexp(x(2*i)*4)-u(2*i))**2end dodo i=1,m-1r=r+2.*h/3.*(dexp(x(2*i+1)*4)-u(2*i+1))**2end doif(dabs(r)<=e) exitend dowrite(*,*) kopen(2,file="simpson.txt")do i=1,2*m+1write(2,'(3(f18.12))') x(i),u(i),dexp(x(i)*4.)end dowrite(2,'(4x,a2,e9.2)') "r=",rclose(2)returnend subroutine simpson!///////////Gauss_Legendre法subroutine Gaussimplicit noneinteger,parameter :: m=8integer :: i,j,kreal*8 :: x(m),u(m),a(m)real*8 :: sum,g,rreal*8 :: e=1.0e-10data x /-0.9602898565,-0.7966664774,-0.5255324099,-0.1834346425,&0.1834346425,0.5255324099,0.7966664774,0.9602898565/data a /0.1012285363,0.2223810345,0.3137066459,0.3626837834,&0.3626837834,0.3137066459,0.2223810345,0.1012285363/u=0.02do k=1,50do i=1,mg=dexp(x(i)*4.)+(dexp(x(i)+4.)-dexp(-4.-x(i)))/(x(i)+4.)do j=1,msum=sum+dexp(x(i)*x(j))*u(j)*a(j)end dou(i)=g-sumend dor=0.do i=1,mr=r+a(i)*(dexp(x(i)*4)-u(i))**2end doif(dabs(r)<=e) exitend dowrite(*,*) kopen(3,file="Gauss.txt")do i=1,mwrite(3,'(3(f18.12))') x(i),u(i),dexp(x(i)*4.)end dowrite(3,'(4x,a2,e9.2)') "r=",rclose(3)returnend subroutine Gaussend module!//////////主程序program mainuse integralimplicit noneinteger :: code1=2600integer :: code2=41call trapezoid(code1)call simpson(code2)call Gaussend program3.各种积分方法的节点和数值解（由于数据太多，在打印时用了较计算时少的有效数字）复化Simpson法4.各方法所得曲线（由于所取节点太多，且精度高，所以图中很难看出各曲线的区别。
北航研究⽣数值分析作业第⼀题北航研究⽣数值分析作业第⼀题：⼀、算法设计⽅案1.要求计算矩阵的最⼤最⼩特征值，通过幂法求得模最⼤的特征值，进⾏⼀定判断即得所求结果；2.求解与给定数值接近的特征值，可以该数做漂移量，新数组特征值倒数的绝对值满⾜反幂法的要求，故通过反幂法即可求得；3.反幂法计算时需要⽅程求解中间过渡向量，需设计Doolite分解求解；4.|A|=|B||C|，故要求解矩阵的秩，只需将Doolite分解后的U矩阵的对⾓线相乘即为矩阵的Det。
算法编译环境：vlsual c++6.0需要编译函数：幂法，反幂法，Doolite分解及⽅程的求解⼆、源程序如下：#include#include#include#includeint Max(int value1,int value2);int Min(int value1,int value2);void Transform(double A[5][501]);double mifa(double A[5][501]);void daizhuangdoolite(double A[5][501],double x[501],double b[501]); double fanmifa(double A[5][501]); double Det(double A[5][501]);/***定义2个判断⼤⼩的函数，便于以后调⽤***/int Max(int value1,int value2){return((value1>value2)?value1:value2);}int Min(int value1,int value2){return ((value1}/*****************************************//***将矩阵值转存在⼀个数组⾥，节省空间***/void Transform(double A[5][501],double b,double c){int i=0,j=0;A[i][j]=0,A[i][j+1]=0;for(j=2;j<=500;j++)A[i][j]=c;i++;j=0;A[i][j]=0;for(j=1;j<=500;j++)A[i][j]=b;i++;for(j=0;j<=500;j++)A[i][j]=(1.64-0.024*(j+1))*sin(0.2*(j+1))-0.64*exp(0.1/(j+1)); i++;for(j=0;j<=499;j++)A[i][j]=b;A[i][j]=0;i++;for(j=0;j<=498;j++)A[i][j]=c;A[i][j]=0,A[i][j+1]=0;}/***转存结束***///⽤于求解模最⼤的特征值，幂法double mifa(double A[5][501]){int s=2,r=2,m=0,i,j;double b2,b1=0,sum,u[501],y[501];for (i=0;i<=500;i++){u[i] = 1.0;}do{sum=0;if(m!=0)b1=b2;m++;for(i=0;i<=500;i++)sum+=u[i]*u[i];for(i=0;i<=500;i++)y[i]=u[i]/sqrt(sum);for(i=0;i<=500;i++){u[i]=0;for(j=Max(i-r,0);j<=Min(i+s,500);j++)u[i]=u[i]+A[i-j+s][j]*y[j];}b2=0;for(i=0;i<=500;i++)b2=b2+y[i]*u[i];}while(fabs(b2-b1)/fabs(b2)>=exp(-12));return b2;}//带状DOOLITE分解，并且求解出⽅程组的解void daizhuangdoolite(double A[5][501],double x[501],double b[501]) { int i,j,k,t,s=2,r=2;double B[5][501],c[501];for(i=0;i<=4;i++){for(j=0;j<=500;j++)B[i][j]=A[i][j];}for(i=0;i<=500;i++)c[i]=b[i];for(k=0;k<=500;k++){for(j=k;j<=Min(k+s,500);j++){for(t=Max(0,Max(k-r,j-s));t<=k-1;t++)B[k-j+s][j]=B[k-j+s][j]-B[k-t+s][t]*B[t-j+s][j]; }for(i=k+1;i<=Min(k+r,500);i++){for(t=Max(0,Max(i-r,k-s));t<=k-1;t++)B[i-k+s][k]=B[i-k+s][k]-B[i-t+s][t]*B[t-k+s][k]; B[i-k+s][k]=B[i-k+s][k]/B[s][k];}}for(i=1;i<=500;i++)for(t=Max(0,i-r);t<=i-1;t++)c[i]=c[i]-B[i-t+s][t]*c[t];x[500]=c[500]/B[s][500];for(i=499;i>=0;i--){x[i]=c[i];for(t=i+1;t<=Min(i+s,500);t++)x[i]=x[i]-B[i-t+s][t]*x[t];x[i]=x[i]/B[s][i];}}//⽤于求解模最⼤的特征值，反幂法double fanmifa(double A[5][501]){int s=2,r=2,m=0,i;double b2,b1=0,sum=0,u[501],y[501];for (i=0;i<=500;i++){u[i] = 1.0;}do{if(m!=0)b1=b2;m++;sum=0;for(i=0;i<=500;i++)sum+=u[i]*u[i];for(i=0;i<=500;i++)y[i]=u[i]/sqrt(sum);daizhuangdoolite(A,u,y);b2=0;for(i=0;i<=500;i++)b2+=y[i]*u[i];}while(fabs(b2-b1)>=fabs(b1)*exp(-12));return 1/b2;}//⾏列式的LU分解，U的主线乘积即位矩阵的DET double Det(double A[5][501]) {int i,j,k,t,s=2,r=2;for(k=0;k<=500;k++){for(j=k;j<=Min(k+s,500);j++){for(t=Max(0,Max(k-r,j-s));t<=k-1;t++)A[k-j+s][j]=A[k-j+s][j]-A[k-t+s][t]*A[t-j+s][j];}for(i=k+1;i<=Min(k+r,500);i++){for(t=Max(0,Max(i-r,k-s));t<=k-1;t++)A[i-k+s][k]=A[i-k+s][k]-A[i-t+s][t]*A[t-k+s][k];A[i-k+s][k]=A[i-k+s][k]/A[s][k];}}double det=1;for(i=0;i<=500;i++)det*=A[s][i];return det;}void main(){double b=0.16,c=-0.064,p,q;int i,j;double A[5][501];Transform(A,b,c); //进⾏A的赋值cout.precision(12); //定义输出精度double lamda1,lamda501,lamdas;double k=mifa(A);if(k>0) //判断求得最⼤以及最⼩的特征值.如果K>0，则它为最⼤特征值值，//并以它为偏移量再⽤⼀次幂法求得新矩阵最⼤特征值，即为最⼤ //与最⼩的特征值的差{lamda501=k;for(i=0;i<=500;i++)A[2][i]=A[2][i]-k;lamda1=mifa(A)+lamda501;for(i=0;i<=500;i++)A[2][i]=A[2][i]+k;}else //如果K<=0，则它为最⼩特征值值，并以它为偏移量再⽤⼀次幂法//求得新矩阵最⼤特征值，即为最⼤与最⼩的特征值的差{lamda1=k;for(i=0;i<=500;i++)A[2][i]=A[2][i]-k;lamda501=mifa(A)+lamda1;for(i=0;i<=500;i++)A[2][i]=A[2][i]+k;}lamdas=fanmifa(A);FILE *fp=fopen("result.txt","w");fprintf(fp,"λ1=%.12e\n",lamda1);fprintf(fp,"λ501=%.12e\n",lamda501);fprintf(fp,"λs=%.12e\n\n",lamdas);fprintf(fp,"\t要求接近的值\t\t\t实际求得的特征值\n");for(i=1;i<=39;i++) //反幂法求得与给定值接近的特征值{p=lamda1+(i+1)*(lamda501-lamda1)/40;for(j=0;j<=500;j++)A[2][j]=A[2][j]-p;q=fanmifa(A)+p;for(j=0;j<=500;j++)A[2][j]=A[2][j]+p;fprintf(fp,"µ%d: %.12e λi%d: %.12e\n",i,p,i,q);}double cond=fabs(mifa(A)/fanmifa(A));double det=Det(A);fprintf(fp,"\ncond(A)=%.12e\n",cond);fprintf(fp,"\ndetA=%.12e\n",det);}三、程序运⾏结果λ1=-1.069936345952e+001λ501=9.722283648681e+000λs=-5.557989086521e-003要求接近的值实际求得的特征值µ1: -9.678281104107e+000 λi1: -9.585702058251e+000µ2: -9.167739926402e+000 λi2: -9.172672423948e+000µ3: -8.657198748697e+000 λi3: -8.652284007885e+000µ4: -8.146657570993e+000 λi4: -8.0934********e+000µ5: -7.636116393288e+000 λi5: -7.659405420574e+000µ6: -7.125575215583e+000 λi6: -7.119684646576e+000µ7: -6.615034037878e+000 λi7: -6.611764337314e+000µ8: -6.104492860173e+000 λi8: -6.0661********e+000µ9: -5.593951682468e+000 λi9: -5.585101045269e+000µ10: -5.0834********e+000 λi10: -5.114083539196e+000µ11: -4.572869327058e+000 λi11: -4.578872177367e+000µ12: -4.062328149353e+000 λi12: -4.096473385708e+000µ13: -3.551786971648e+000 λi13: -3.554211216942e+000µ14: -3.0412********e+000 λi14: -3.0410********e+000µ15: -2.530704616238e+000 λi15: -2.533970334136e+000µ16: -2.020*********e+000 λi16: -2.003230401311e+000µ17: -1.509622260828e+000 λi17: -1.503557606947e+000µ18: -9.990810831232e-001 λi18: -9.935585987809e-001µ19: -4.885399054182e-001 λi19: -4.870426734583e-001µ20: 2.200127228676e-002 λi20: 2.231736249587e-002µ21: 5.325424499917e-001 λi21: 5.324174742068e-001µ22: 1.043083627697e+000 λi22: 1.052898964020e+000µ23: 1.553624805402e+000 λi23: 1.589445977158e+000µ24: 2.064165983107e+000 λi24: 2.060330427561e+000µ25: 2.574707160812e+000 λi25: 2.558075576223e+000µ26: 3.0852********e+000 λi26: 3.080240508465e+000µ27: 3.595789516221e+000 λi27: 3.613620874136e+000µ28: 4.106330693926e+000 λi28: 4.0913********e+000µ29: 4.616871871631e+000 λi29: 4.603035354280e+000µ30: 5.127413049336e+000 λi30: 5.132924284378e+000µ31: 5.637954227041e+000 λi31: 5.594906275501e+000µ32: 6.148495404746e+000 λi32: 6.080933498348e+000µ33: 6.659036582451e+000 λi33: 6.680354121496e+000µ34: 7.169577760156e+000 λi34: 7.293878467852e+000µ35: 7.680118937861e+000 λi35: 7.717111851857e+000µ36: 8.190660115566e+000 λi36: 8.225220016407e+000µ37: 8.701201293271e+000 λi37: 8.648665837870e+000µ38: 9.211742470976e+000 λi38: 9.254200347303e+000µ39: 9.722283648681e+000 λi39: 9.724634099672e+000cond(A)=1.925042185755e+003detA=2.772786141752e+118四、分析如果初始向量选择不当,将导致迭代中X1的系数等于零.但是，由于舍⼊误差的影响,经若⼲步迭代后,.按照基向量展开时,x1的系数可能不等于零。
北航数值分析实验报告篇一：北航数值分析报告第一大题《数值分析》计算实习报告第一大题学号：DY1305姓名：指导老师：一、题目要求已知501*501阶的带状矩阵A，其特征值满足?1?2...?501。
试求：1、?1，?501和?s的值；2、A的与数?k??1?k?501??140最接近的特征值?ik（k=1，2，...，39）；3、A的（谱范数）条件数c nd(A)2和行列式de tA。
二、算法设计方案题目所给的矩阵阶数过大，必须经过去零压缩后进行存储和运算，本算法中压缩后的矩阵A1如下所示。
?0?0?A1??a1??b??c0b a2bcc bb c............c bb ccb a500b0a 3...a499c?b??a501??0?0??由矩阵A的特征值满足的条件可知?1与?501之间必有一个最大，则采用幂法求出的一个特征值必为其中的一个：当所求得的特征值为正数，则为?501；否则为?1。
在求得?1与?501其中的一个后，采用带位移的幂法则可求出它们中的另一个，且位移量即为先求出的特征值的值。
用反幂法求得的特征值必为?s。
由条件数的性质可得，c nd(A)2为模最大的特征值与模最小的特征值之比的模，因此，求出?1，?501和?s的值后，则可以求得c nd(A)2。
一、题目：关于x, y, t, u, v, w 的下列方程组0.5cos 2.670.5sin 1.070.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 +++-=⎧⎪+++-=⎪⎨+++-=⎪⎪+++-=⎩1、试用数值方法求出f(x, y)在区域 {(,)|00.8,0.5 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、计算****(,),(,)i j i j f x y p x y （i = 1, 2, …，8；j = 1, 2,…，5）的值，以观察(,)p x y 逼近(,)f x y 的效果，其中，*i x =0.1i , *j y =0.5+0.2j 。
说明：1、用迭代方法求解非线性方程组时，要求近似解向量()k x 满足()(1)()12||||/||||10k k k x x x --∞∞-≤2、作二元插值时，要使用分片二次代数插值。
3、要由程序自动确定最小的k 值。
4、打印以下内容：●算法的设计方案。
●全部源程序（要求注明主程序和每个子程序的功能）。
●数表：,,i j x y (,)i j f x y （i = 0,1,2,…,10；j = 0,1,2,…,20）。
●选择过程的,k σ值。
●达到精度要求时的,k σ值以及(,)p x y 中的系数rs c （r = 0,1,…,k;s = 0,1,…,k ）。
●数表：**,,i j x y ****(,),(,)i j i j f x y p x y （i = 1, 2, ...,8；j = 1, 2, (5)。
北航数值分析全部三次大作业第一次大作业是关于解线性方程组的数值方法。
我们被要求实现各种常用的线性方程组求解算法，例如高斯消元法、LU分解法和迭代法等。
我首先学习了这些算法的原理和实现方法，并借助Python编程语言编写了这些算法的代码。
在实验中，我们使用了不同规模和条件的线性方程组进行测试，并比较了不同算法的性能和精度。
通过这个作业，我深入了解了线性方程组求解的原理和方法，提高了我的编程和数值计算能力。
第二次大作业是关于数值积分的方法。
数值积分是数值分析中的重要内容，它可以用于计算曲线的长度、函数的面积以及求解微分方程等问题。
在这个作业中，我们需要实现不同的数值积分算法，例如矩形法、梯形法和辛普森法等。
我学习了这些算法的原理和实现方法，并使用Python编写了它们的代码。
在实验中，我们计算了不同函数的积分值，并对比了不同算法的精度和效率。
通过这个作业，我深入了解了数值积分的原理和方法，提高了我的编程和数学建模能力。
第三次大作业是关于常微分方程的数值解法。
常微分方程是数值分析中的核心内容之一，它可以用于描述众多物理、化学和生物现象。
在这个作业中，我们需要实现不同的常微分方程求解算法，例如欧拉法、龙格-库塔法和Adams法等。
我学习了这些算法的原理和实现方法，并使用Python编写了它们的代码。
在实验中，我们解决了一些具体的常微分方程问题，并比较了不同算法的精度和效率。
通过这个作业，我深入了解了常微分方程的原理和方法，提高了我的编程和问题求解能力。
总的来说，北航数值分析课程的三次大作业非常有挑战性，但也非常有意义。
通过这些作业，我在数值计算和编程方面得到了很大的提升，也更加深入地了解了数值分析的理论和方法。
虽然这些作业需要大量的时间和精力，但我相信这些努力将会对我未来的学习和工作产生积极的影响。
北京航空航天大学2020届研究生《数值分析》实验作业第九题院系：xx学院学号：姓名：2020年11月Q9:方程组A.4一、 算法设计方案（一）总体思路1.题目要求∑∑===k i kj s r rsy x cy x p 00),(对f(x, y) 进行拟合，可选用乘积型最小二乘拟合。
),(i i y x 与),(i i y x f 的数表由方程组与表A-1得到。
2.),(**j i y x f 与1使用相同方法求得，),(**j i y x p 由计算得出的p(x,y)直接带入),(**j i y x 求得。
1. ),(i i y x 与),(i i y x f 的数表的获得对区域D ={ （x,y)|1≤x ≤1.24,1.0≤y ≤1.16}上的f (x , y )值可通过xi=1+0.008i ，yj=1+0.008j ，得到),(i i y x 共31×21组。
将每组带入A4方程组，即可获得五个二元函数组，通过简单牛顿迭代法求解这五个二元数组可获得z1~z5有关x,y 的表达式。
再将),(i i y x 分别带入z1~z5表达式即可获得f(x,y)值。
2.乘积型最小二乘曲面拟合2.1使用乘积型最小二乘拟合，根据k 值不用，有基函数矩阵如下：⎪⎪⎪⎭⎫ ⎝⎛=k i i k x x x x B 0000 ， ⎪⎪⎪⎭⎫ ⎝⎛=k j jk y y y y G 0000数表矩阵如下：⎪⎪⎪⎭⎫⎝⎛=),(),(),(),(0000j i i j y x f y x f y x f y x f U记C=[rs c ]，则系数rs c 的表达式矩阵为：11-)(-=G G UG B B B C T TT ）（通过求解如下线性方程，即可得到系数矩阵C 。
UG B G G C B B T T T =)()(2.2计算),(),,(****j i j i y x p y x f (i =1,2,…,31 ; j =1,2,…,21) 的值),(**j i y x f 的计算与),(j i y x f 相同。
数值分析大作业一、算法设计方案1、矩阵初始化矩阵[]501501⨯=ij a A 的下半带宽r=2,上半带宽s=2，设置矩阵[][]5011++s r C ，在矩阵C 中检索矩阵A 中的带内元素ij a 的方法是：j s j i ij c a ,1++-=。
这样所需要的存储单元数大大减少，从而极大提高了运算效率。
2、利用幂法求出5011λλ，幂法迭代格式：0111111nk k k k kk T k k k u R y u u Ay y u ηηβ------⎧∈⎪⎪=⎪=⎨⎪=⎪⎪=⎩非零向量 当1210/-≤-k k βββ时，迭代终止。
首先对于矩阵A 利用幂法迭代求出一个λ，然后求出矩阵B ，其中I A B λ-=(I 为单位矩阵），对矩阵B 进行幂法迭代，求出λ'，之后令λλλ+'=''，比较的大小与λλ''，大者为501λ，小者为1λ。
3、利用反幂法求出iks λλ,反幂法迭代格式：0111111nk k k k kk T k k k u R y u Au y y u ηηβ------⎧∈⎪⎪=⎪=⎨⎪=⎪⎪=⎩非零向量 当1210/-≤-k k βββ时，迭代终止，1s k λβ=。
每迭代一次都要求解一次线性方程组1-=k k y Au ，求解过程为：（1）作分解LU A =对于n k ,...,2,1=执行[][]s k n r k k k i c c c c c n s k k k j c cc c k s ks k t k s k r i t t s t i k s k i k s k i js j t k s j r k t t s t k j s j k j s j k <+++=-=++=-=+++----=++-++-++-++----=++-++-++-∑∑);,min(,...,2,1/)(:),min(,...,1,:,1,11),,1max(,1,1,1,11),,1max(,1,1,1（2）求解y Ux b Ly ==,（数组b 先是存放原方程组右端向量，后来存放中间向量y))1,...,2,1(/)(:/:),...,3,2(:,1),min(1.1.11),1max(,1--=-===-=+++-++-+--=++-∑∑n n i c x c b x c b x n i b c b b i s t n s i i t t s t i i i ns n n ti r i t t s t i i i使用反幂法，直接可以求得矩阵按模最小的特征值s λ。
北航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. 下列各数都是经过四舍五入得到的近似数，试指出它们有几位有效数字以及它们的绝对误差限、相对误差限。
（1）；（2）；（3）；（4）；（5）；（6）；（7）；1. （1）5，，；（2）2，，；（3）4，，；（4）5，，；（5）1，，；（6）2，，；（7）6，，2. 为使下列各数的近似值的相对误差限不超过，问各近似值分别应取几位有效数字？2. ；；3. 设均为第1题所给数据，估计下列各近似数的误差限。
（1）；（2）；（3）3. （1）；（2）；（3）4. 计算，取，利用下列等价表达式计算，哪一个的结果最好？为什么？（1）；（2）；（3）（4）4. 第（3）个结果最好。
5. 序列满足递推关系式若（三位有效数字），计算时误差有多大？这个计算过程稳定吗？5. 不稳定。
从计算到时，误差约为6. 求方程的两个根，使其至少具有四位有效数字（要求利用。
6. ，7. 某生产部门生产的一件产品需用七个零件，而这七个零件的质量取决于零件参数的标定值，它们的参数允许有一定的误差：若每一零件的标定值取做区间中点，在生产过程中每一零件的参数都有可能产生误差。
由此将零件分成不同的等级：A，B，C三等，等级由标定值的相对误差限表示，A等为1%，B等为5%，C等为10%。
试确定三个等级的零件分别满足的区间。
8. 将一个八位二进制数(10111101)2转换成十进制数时，可以用公式：（1）用多项式求值的秦九韶方法求C的值；（2）写出将任意一个八位二进制数(b1b2b3b4b5b6b7b8)2转化为十进制数的算法。
9. 利用等式变换使下列表达式的计算结果比较精确。
（1）；（2）（3）；（4）9. （1）；（2）；（3）；（4）10. 设，求证：（1）（2）利用（1）中的公式正向递推计算时误差增大；反向递推时误差函数减小。
习题二1. 判断下列方程有几个实根，并求出其隔根区间。
（1）；（2）（3）；（4）1. （1），，；（2）；（3），，；（4）为根。
数值分析第二次大作业史立峰SY1505327一、 方案（1）利用循环结构将sin(0.50.2)()1.5cos( 1.2)(){i j i j ij i j i j a +≠+==(i,j=1,2,……,10)进行赋值，得到需要变换的矩阵A ；（2）然后，对矩阵A 利用Householder 矩阵进行相似变换，把A 化为上三角矩阵A(n-1)。
对A 拟上三角化，得到拟上三角矩阵A(n-1)，具体算法如下：记A(1)=A ，并记A(r)的第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)(Λ++=全为零，则令A(r+1) =A(r),转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. 令()nTr 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 -=ωT rr T r r r r p u u A A --=+ω)()1(5. 继续。
（3）使用带双步位移的QR 方法计算矩阵A (n-1)的全部特征值，也是A 的全部特征值，具体算法如下：1. 给定精度水平0>ε和迭代最大次数L 。
2. 记n n ij n a A A ⨯-==][)1()1()1(，令n m k ==,1。
3. 如果ε≤-)(1,k m m a ，则得到A 的一个特征值)(k mm a ，置1:-=m m （降阶），转4；否则转5。
数值分析大作业一、算法设计方案1、矩阵初始化矩阵[]501501⨯=ij a A 的下半带宽r=2,上半带宽s=2，设置矩阵[][]5011++s r C ，在矩阵C 中检索矩阵A 中的带内元素ij a 的方法是：j s j i ij c a ,1++-=。
这样所需要的存储单元数大大减少，从而极大提高了运算效率。
2、利用幂法求出5011λλ，幂法迭代格式：011111111nT k k k k k k kk T k k k u R u u y u u Ay y u ηηβ--------⎧∈⎪⎪=⎪=⎨⎪=⎪⎪=⎩非零向量 当1210/-≤-k k βββ时，迭代终止。
首先对于矩阵A 利用幂法迭代求出一个λ，然后求出矩阵B ，其中I A B λ-=(I 为单位矩阵），对矩阵B 进行幂法迭代，求出λ'，之后令λλλ+'=''，比较的大小与λλ''，大者为501λ，小者为1λ。
3、利用反幂法求出iks λλ,反幂法迭代格式： 011111111nTk k k k k k kk T k k k u R u u y u Au y y u ηηβ--------⎧∈⎪⎪=⎪=⎨⎪=⎪⎪=⎩非零向量 当1210/-≤-k k βββ时，迭代终止，1s k λβ=。
每迭代一次都要求解一次线性方程组1-=k k y Au ，求解过程为：（1）作分解LU A =对于n k ,...,2,1=执行[][]s k n r k k k i c c c c c n s k k k j c cc c k s ks k t k s k r i t t s t i k s k i k s k i js j t k s j r k t t s t k j s j k j s j k <+++=-=++=-=+++----=++-++-++-++----=++-++-++-∑∑);,min(,...,2,1/)(:),min(,...,1,:,1,11),,1max(,1,1,1,11),,1max(,1,1,1（2）求解y Ux b Ly ==,（数组b 先是存放原方程组右端向量，后来存放中间向量y))1,...,2,1(/)(:/:),...,3,2(:,1),min(1.1.11),1max(,1--=-===-=+++-++-+--=++-∑∑n n i c x c b x c b x n i b c b b i s t n s i i t t s t i i i ns n n ti r i t t s t i i i使用反幂法，直接可以求得矩阵按模最小的特征值s λ。
《数值分析》第三次大作业一、算法的设计方案： （一）、总体方案设计：（1）解非线性方程组。
将给定的(,)i i x y 当作已知量代入题目给定的非线性方程组，求得与(,)i i x y 相对应的数组t[i][j],u[i][j]。
（2）分片二次代数插值。
通过分片二次代数插值运算，得到与数组t[11][21],u[11][21]]对应的数组z[11][21]，得到二元函数z=(,)i i f x y 。
（3）曲面拟合。
利用x[i],y[j],z[11][21]建立二维函数表，再根据精度的要求选择适当k 值，并得到曲面拟合的系数矩阵C[r][s]。
（4）观察和(,)i i p x y 的逼近效果。
观察逼近效果只需要重复上面（1）和（2）的过程，得到与新的插值节点(,)i i x y 对应的(,)i i f x y ，再与对应的(,)i i p x y 比较即可，这里求解(,)i i p x y 可以直接使用（3）中的C[r][s]和k 。
（二）具体算法设计：（1）解非线性方程组牛顿法解方程组()0F x =的解*x ，可采用如下算法： 1）在*x 附近选取(0)xD ∈，给定精度水平0ε>和最大迭代次数M 。
2）对于0,1,k M =L 执行① 计算()()k F x和()()k F x '。
② 求解关于()k x∆的线性方程组()()()()()k k k F x x F x '∆=-③ 若()()k k x x ε∞∞∆≤，则取*()k x x ≈，并停止计算；否则转④。
④ 计算(1)()()k k k xx x +=+∆。
⑤ 若k M <，则继续，否则，输出M 次迭代不成功的信息，并停止计算。
（2）分片双二次插值给定已知数表以及需要插值的节点，进行分片二次插值的算法：设已知数表中的点为： 00(0,1,,)(0,1,,)i j x x ih i n y y j j m τ=+=⎧⎪⎨=+=⎪⎩L L ，需要插值的节点为(,)x y 。
“数值分析“计算实习大作业第二题——SY1415215孔维鹏一、计算说明本程序采用带双步位移的QR方法求解矩阵A的所有特征值，然后采用反幂法求解矩阵A的实特征值对应的特征向量。
在采用带双步位移的QR方法求解特征值时，对教材上所提供的具体算法作稍微的改动，以简化程序，具体算法如下所示：1、计算出A拟上三角化后的矩阵，给定精度水平和最大迭代次数L；2、记，令k=1，m=n；3、如果，则可直接计算出最后1或2个特征值，转8，否则转4；4、如果，则可得一个特征值，置m=m-1；转3，否则转5；5、如果，则可得两个特征值，置m=m-2；转3，否则转6；6、记，计算7、k=k+1,转38、A的全部特征值已经求出，停止计算。
二、计算源程序（FORTRAN）PROGRAM SY1415215_2PARAMETER (N=10)DIMENSION A(N,N),A1(N,N),A2(N,N),C(2,N),Q(N,N),R(N,N),CR(N),CM(N)!C为存储特征值的数组,1为实部，为虚部REAL(8) A,A1,A2,C,Q,R,CME=1E-12 !精度水平L=1000 !迭代最大次数OPEN(1,FILE='数值分析大作业第二题计算结果.TXT')DO I=1,NDO J=1,NIF(I==J) THENA(I,J)=*COS(I+*J)ELSEA(I,J)=SIN*I+*J)ENDIFENDDOENDDOA1=AWRITE(*,"('矩阵A为：')")WRITE(1,"('矩阵A为：')")DO I=1,NDO J=1,NWRITE(*,"(2X,,2X,\)") A(I,J)WRITE(1,"(2X,,2X,\)") A(I,J)ENDDOWRITE(*,"(' ')")WRITE(1,"(' ')")ENDDO!使用矩阵的拟上三角化的算法将矩阵A化为拟上三角矩阵A（n-1）CALL HESSENBERG(A,N)WRITE(*,"('拟上三角化后矩阵A（n-1）为：')")WRITE(1,"('拟上三角化后矩阵A（n-1）为：')")DO I=1,NDO J=1,NWRITE(*,"(2X,,2X,\)") A(I,J)WRITE(1,"(2X,,2X,\)") A(I,J)ENDDOWRITE(*,"('')")WRITE(1,"('')")ENDDO!计算对矩阵A（n-1）实行QR方法迭代结束后所得矩阵A2=ACALL QRD(A2,N,Q,R)WRITE(*,"('对矩阵A（n-1）实行QR方法迭代结束后所得Q为：')") WRITE(1,"('对矩阵A（n-1）实行QR方法迭代结束后所得Q为：')") DO I=1,NDO J=1,NWRITE(*,"(2X,,2X,\)") Q(I,J)WRITE(1,"(2X,,2X,\)") Q(I,J)ENDDOWRITE(*,"('')")WRITE(1,"('')")ENDDOWRITE(*,"('对矩阵A（n-1）实行QR方法迭代结束后所得R为：')") WRITE(1,"('对矩阵A（n-1）实行QR方法迭代结束后所得R为：')") DO I=1,NDO J=1,NWRITE(*,"(2X,,2X,\)") R(I,J)WRITE(1,"(2X,,2X,\)") R(I,J)ENDDOWRITE(*,"('')")WRITE(1,"('')")ENDDO!使用带双步位移的QR方法求解矩阵A（n-1）的特征值K=1M=NDO WHILE(K<=L)IF(M<=2) THENIF(M==1) THENC(1,M)=A(M,M)ELSE IF(M==2) THENCALL CALCUS(A,N,M,C)ENDIFEXITELSE IF(ABS(A(M,M-1))<E) THENC(1,M)=A(M,M)M=M-1ELSE IF(ABS(A(M-1,M-2))<E) THENCALL CALCUS(A,N,M,C)M=M-2ELSECALL CALM(A,M,N)ENDIFK=K+1ENDDOWRITE(*,"('矩阵A的全部特征值为：')")WRITE(1,"('矩阵A的全部特征值为：')")DO J=1,NWRITE(*,",'+',,'i')") C(1,J),C(2,J)WRITE(1,",'+',,'i')") C(1,J),C(2,J)ENDDO!使用反幂法求解A的相应于实特征值的特征向量J=1DO I=1,NIF(C(2,I)==0)THENCR(J)=C(1,I)J=J+1ENDIFENDDOJC=J-1WRITE(*,"('矩阵A的实特征值为：')")WRITE(1,"('矩阵A的实特征值为：')")DO I=1,JCWRITE(*,"") CR(I)WRITE(1,"") CR(I)ENDDODO II=1,JCDO I=1,NDO J=1,NIF(I==J) THENA(I,J)=A1(I,J)-CR(II)ELSEA(I,J)=A1(I,J)ENDIFENDDOENDDOCALL INPOVERMETHOD(A,N,CM)WRITE(*,"('与实特征值',,'对应的特征向量为：')") CR(II) WRITE(1,"('与实特征值',,'对应的特征向量为：')") CR(II) DO I=1,NWRITE(*,"(2X,,2X,\)") CM(I)WRITE(1,"(2X,,2X,\)") CM(I)ENDDOWRITE(*,"('')")WRITE(1,"('')")ENDDOCLOSE(1)END!***************拟上三角化子函数*************************！SUBROUTINE HESSENBERG(A,N)DIMENSION A(N,N),P(N),Q(N),W(N),U(N),AT(N,N)REAL(8) A,P,Q,W,U,ATREAL(8) S0,S1,S2,S3,S4,TDO L=1,N-2JUDGE=0DO I=L+2,NIF(A(I,L)/=0) THENJUDGE=1EXITENDIFENDDOIF(JUDGE==0) THENA=ACYCLEELSE IF(JUDGE/=0) THEN!计算DRS0=0DO I=L+1,NS0=S0+A(I,L)**2ENDDODR=SQRT(S0)!计算CRIF(A(L+1,L)==0)THENCR=DRELSECR=-SGN(A(L+1,L))*DRENDIF!计算HRHR=CR**2-CR*A(L+1,L)!给u赋值IF(I<L+1) THENU(I)=0ELSE IF(I==L+1) THEN U(I)=A(I,L)-CRELSE IF(I>L+1) THEN U(I)=A(I,L)ENDIFENDDO!计算PDO I=1,NDO J=1,NAT(I,J)=A(J,I)ENDDOENDDODO I=1,NS1=0DO J=1,NS1=S1+AT(I,J)*U(J)ENDDOP(I)=S1/HRENDDO!计算QDO I=1,NS2=0DO J=1,NS2=S2+A(I,J)*U(J)ENDDOQ(I)=S2/HRENDDO!计算TS3=0DO I=1,NS3=S3+P(I)*U(I)ENDDOT=S3/HR!计算WDO I=1,NW(I)=Q(I)-T*U(I)!计算A（r+1）DO I=1,NDO J=1,NA(I,J)=A(I,J)-W(I)*U(J)-U(I)*P(J)ENDDOENDDOENDIFENDDORETURNEND!***************符号函数子程序*****************！FUNCTION SGN(X)REAL(8) XIF(X>0) THENSGN=1ELSE IF(X<0) THENSGN=-1ELSE IF(X==0) THENSGN=0ENDIFEND!*********计算二阶子阵特征值s1，s2子函数*******！SUBROUTINE CALCUS(X,N,M,Y)DIMENSION X(N,N),Y(2,N)REAL(8) A,B,C,D,X,YA=1C=X(M-1,M-1)*X(M,M)-X(M-1,M)*X(M,M-1)B=-(X(M-1,M-1)+X(M,M))D=B**2-4*CIF(D>=0) THENY(1,M)=(-B-SQRT(D))/2Y(1,M-1)=(-B+SQRT(D))/2ELSEIF(D<0) THENY(1,M)=-B/2Y(1,M-1)=-B/2Y(2,M)=-SQRT(-D)/2Y(2,M-1)=-SQRT(-D)/2ENDIFRETURNEND!*********计算Mk，Ak+1子函数************！SUBROUTINE CALM(A,M,N)DIMENSION A(N,N),MK(M,M),X(M,M),QK(M,M),RK(M,M),S1(M,M),S2(M,M),QKT(M,M) REAL(8) A,MK,X,QK,RK,QKTREAL(8) S0,S1,S2DO I=1,MDO J=1,MIF(I==J) THENX(I,J)=1ELSEX(I,J)=0ENDIFENDDOENDDOS=A(M-1,M-1)+A(M,M)T=A(M-1,M-1)*A(M,M)-A(M,M-1)*A(M-1,M)DO I=1,MDO J=1,MS0=0DO K=1,MS0=S0+A(I,K)*A(K,J)ENDDOMK(I,J)=S0-S*A(I,J)+T*X(I,J)ENDDOENDDO!对Mk做QR分解CALL QRD(MK,M,QK,RK)DO I=1,MDO J=1,MQKT(I,J)=QK(J,I)ENDDOENDDODO I=1,MDO J=1,MS1(I,J)=0DO K=1,MS1(I,J)=S1(I,J)+QKT(I,K)*A(K,J)ENDDOENDDOENDDOA=S1DO I=1,MDO J=1,MS2(I,J)=0DO K=1,MS2(I,J)=S2(I,J)+A(I,K)*QK(K,J)ENDDOENDDOENDDOA=S2RETURNEND!************QR分解子程序***************！SUBROUTINE QRD(A,N,Q,R)DIMENSION A(N,N),AT(N,N),Q(N,N),U(N),W(N),P(N),R(N,N) REAL(8) A,AT,Q,U,W,P,RREAL(8) DR,S0,CR,HR,S1,S2DO I=1,NDO J=1,NIF(I==J) THENQ(I,J)=1ELSEQ(I,J)=0ENDIFENDDOENDDODO L=1,N-1JUDGE=0DO I=L+1,NIF(A(I,L)/=0) THENJUDGE=1EXITENDIF!A(I,L)中有一个不为零,判断条件为真，跳出循环转ENDDOIF(JUDGE==0) THENQ=QA=ACYCLE!A(I,L)全为零，结束本循环，进入下一个ELSE IF(JUDGE/=0) THEN!计算DRS0=0DO I=L,NS0=S0+A(I,L)**2ENDDODR=SQRT(S0)!计算CRIF(A(L,L)==0)THENCR=DRELSECR=-SGN(A(L,L))*DRENDIF!计算HRHR=CR**2-CR*A(L,L)!给u赋值DO I=1,NIF(I<L) THENU(I)=0ELSE IF(I==L) THENU(I)=A(I,L)-CRELSE IF(I>L) THENU(I)=A(I,L)ENDIFENDDO!计算WDO I=1,NS1=0DO J=1,NS1=S1+Q(I,J)*U(J)ENDDOW(I)=S1ENDDO!计算Q（r+1）DO I=1,NDO J=1,NQ(I,J)=Q(I,J)-W(I)*U(J)/HRENDDOENDDO!计算PDO I=1,NDO J=1,NAT(I,J)=A(J,I)ENDDOENDDODO I=1,NS2=0DO J=1,NS2=S2+AT(I,J)*U(J)ENDDOP(I)=S2/HRENDDO!计算A（r+1）DO I=1,NDO J=1,NA(I,J)=A(I,J)-U(I)*P(J)ENDDOENDDOENDIFENDDOQ=QR=ARETURNEND!*************运用反幂法求解矩阵A实特征值的特征向量***********！SUBROUTINE INPOVERMETHOD(A,N,Y)DIMENSION A(N,N),U(N),Y(N),U1(N,N),L1(N,N)REAL(8) E,Z,Z1,Z2,S1,S2,BREAL(8) A,U,Y,U1,L1U(I)=1ENDDO!任取非零向量UCALL DETA(A,N,U1,L1)Z2=EIZ1=E=1K=1DO WHILE (E>1E-12)S1=0DO I=1,NS1=S1+U(I)**2ENDDOB=SQRT(S1) !1DO I=1,NY(I)=U(I)/BENDDO!2CALL DOOLITTLE(U1,L1,Y,N,U) !3利用DOOLITTLE分解法法求解Au=yS2=0DO I=1,NS2=S2+Y(I)*U(I)ENDDOZ1=Z2Z2=S2 !4E=ABS(1/Z2-1/Z1)/ABS(1/Z2) !判断是否满足精度K=K+1ENDDORETURNENDSUBROUTINE DOOLITTLE(U,L,B1,N,X)DIMENSION B(N),U(N,N),X(N),Y(N),B1(N)REAL(8) L(N,N)REAL(8) B,U,X,Y,B1REAL(8) S1,S2,S3,S4B=B1Y(1)=B(1)S3=0DO M=1,I-1S3=S3+L(I,M)*Y(M)ENDDOY(I)=B(I)-S3ENDDOX(N)=Y(N)/U(N,N)DO I=N-1,1,-1S4=0DO M=I+1,NS4=S4+U(I,M)*X(M)ENDDOX(I)=(Y(I)-S4)/U(I,I)ENDDORETURNENDSUBROUTINE DETA(A1,N,U,L) DIMENSION A(N,N),U(N,N),A1(N,N) REAL(8) L(N,N)REAL(8) X,S1,S2REAL(8) A,U,A1X=1A=A1!对矩阵A进行Doolittle分解DO K=1,NDO J=K,NS1=0DO M=1,K-1S1=S1+L(K,M)*U(M,J)ENDDOU(K,J)=A(K,J)-S1A(K,J)=U(K,J)ENDDOIF (K==N) THENEXITELSEDO I=K+1,NS2=0DO M=1,K-1S2=S2+L(I,M)*U(M,K)ENDDOL(I,K)=(A(I,K)-S2)/U(K,K)A(I,K)=L(I,K)ENDDOENDIFENDDORETURNEND三、计算结果矩阵A为：+00 +00 +00 +00 +00 +00 +00 +00 +00 +00+00 +00 +00 +00 +00 +00 +00 +00 +00 +00+00 +00 +01 +00 +00 +00 +00 +00 +00+00 +00 +00 +01 +00 +00 +00 +00 +00 +00 +00 +00 +00 +00 +00 +00 +00+00 +00 +00 +00 +01 +00 +00 +00 +00+00 +00 +00 +00 +00 +00 +01 +00 +00 +00 +00 +00 +00 +00 +00 +00 +00 +00 +00 +00 +00 +00 +00 +00 +00 +00 +00 +00 +00 +00+00 +00 +00 +00 +00 +00 +00 +00 +01拟上三角化后矩阵A（n-1）为：+00 +01 +00 +00 +01 +00 +00 +00+01 +01 +01 +00 +00 +01 +00 +01 +00 +00+01 +01 +01 +00 +01 +00 +00 +00 +00+01 +01 +01 +01 +00 +00 +00 +00+01 +00 +00 +00 +00 +00+00 +01 +00 +00 +00 +00+00 +00 +00+00 +00 +00 +00+00 +00 +00+00 +00对矩阵A（n-1）实行QR方法迭代结束后所得Q为：+00 +00 +00 +00 +00 +00+00 +00 +00 +00+00 +00 +00 +00+00 +00 +00 +00 +00 +00+00 +00+00 +00 +00 +00 +00+00 +00 +00 +00+00+00 +00 +00+00 +00对矩阵A（n-1）实行QR方法迭代结束后所得R为：+01 +01 +01 +00 +01 +00 +00 +00 +00 +01 +00 +00 +00 +00 +00 +00 +00 +00 +01 +01 +00 +01 +00 +01 +00 +00 +00 +00 +01 +00 +00 +00 +00+00 +00 +01 +01 +00 +00 +00+00 +00 +01 +00 +00 +00 +00+00 +00 +00 +00 +00 +00 +00+00 +00 +00 +00 +01 +00 +00+00 +00 +00 +00 +00 +00+00 +00 +00 +00矩阵A的全部特征值为：+01+ +00i+01++00i+01++00i+01+ +00i+01+ +00i+00++00i+00++00i+00+ +00i+00+ +00i+ +00i矩阵A的实特征值为：+01+01+01+00+00与实特征值 +01对应的特征向量为：+00 +00 +00 +00 +00 +00 +00 +00 +00与实特征值 +01对应的特征向量为：+00 +00 +00 +00 +00 +00 +00与实特征值+01对应的特征向量为：+00 +00 +00与实特征值 +00对应的特征向量为：+00 +00 +00 +00 +00与实特征值 +00对应的特征向量为：+00 +00 +00 +00 +00 +00 +00与实特征值对应的特征向量为：+00 +00 +00 +00 +00 +00 +00 +00。
北航研究生数值分析作业第二题北航研究生数值分析作业第二题：一、算法设计方案1.按照题目给出的矩阵定义对矩阵A赋初值：对应的函数为a_init()；2.对矩阵A进行householder变换，使其拟上三角化：对应的函数为householder()；3.输出拟上三角化后的A：对应的函数为aout(int)；4.对拟上三角化后的矩阵A使用带双步位移的QR分解法逐次迭代（最大迭代次数L=500），逐个求出其特征值，对应的函数为eigen_a()；中间包含两个子程序：calc_mk()和qr_analyze()，分别用来计算矩阵M k和对M k进行QR 分解并得到A k+1;5.输出QR分解过程完毕后的A及求得的特征向量：对应的函数为aout()和eigenvalout()；6.对于在第三步中求得的每个实特征值，使用带原点平移的反幂法求出其对应的特征向量，对应的函数为eigenvec()；其中包含一个解方程(A-μI)=y k-1的程序段。
这部分也用迭代完成，仍然将最大迭代次数L设置为500；7.输出矩阵A的特征向量，结束计算：对应的函数为eigenvecout()。
算法编译环境：vlsual c++6.0二、源程序如下：#include#include#define N 10 //矩阵阶数；#define EPSL 1.0e-12 //迭代的精度水平；#define L 500 //迭代最大次数；#define OUTPUTMODE 1 //输出格式：0--输出至屏幕，1--输出至文件double a[N][N], a2[N][N], eigen[N][N]; //声明矩阵A；double sa_re[N] = {0}, sa_im[N] = {0}; //声明矩阵的特征值数组；double u_init[N] = {2,1,2,1,2,1,2,1,2,1}; //定义反幂法中使用的初始向量u；//主程序开始；int main(){FILE *p;void a_init();void householder();void equal_zero(double matrix[N][N], int);void eigenvec();int eigen_a();void aout(int);void eigenvalout(int);void eigenvecout(int);if(OUTPUTMODE){p = fopen("Result.txt", "w+");fprintf(p, "计算结果：\n");fclose(p);}a_init(); //对矩阵A进行初始化；householder(); //对矩阵A进行拟上三角化；equal_zero(a, N); //对矩阵A的元素进行归零处理，消除误差；aout(OUTPUTMODE); //输出A；if(eigen_a()) printf("迭代超过最大次数，特征值求解结果可能不正确。
一、算法设计方案1、解非线性方程组将各拟合节点（x i ，y j ）分别带入非线性方程组，求出与(,)i i x y 相对应的数组te[i][j]，ue[i][j]，求解非线性方程组选择Newton 迭代法，迭代过程中需要求解线性方程组，选择选主元的Doolittle 分解法。
2、二元二次分偏插值对数表z(t,u)进行分片二次代数插值，求得对应（t ij ，u ij ）处的值，即为),(j i y x f 的值。
根据给定的数表，可将整个插值区域分成 16 个小 的区域，故先判断(t i j , u ij ) 所在，的区域，再作此区域的插值，计算 z ij ，相应的Lagrange 形式的插值多项式为：°112211(,)()()(,)m n krkrk m r n p t u l t l u f t u ++=-=-=∑∑其中11()m wk w m k ww kt t l t t t +=-≠-=-∏ (k=m-1, m, m+1) °11()n wr w n r ww ry y l u y y +=-≠-=-∏ (r=n-1, n, n+1)3、曲面拟合从k=1开始逐渐增大k 的值，使用最小二乘法曲面拟合法对z=f(x,y)进行拟合，当710-<σ时结束计算。
拟合基函数φr (x)ψs (y)选择为φr (x)=x r ，ψs (y)=y s 。
拟合系数矩阵c 通过连续两次解线性方程组求得。
[]rsc *=C ，11()()T T T --=C B B B UG G G其中011101011[()]1kk r i k x x x x x x x ϕ⎡⎤⎢⎥⎢⎥==⎢⎥⎢⎥⎢⎥⎣⎦B L LM M M M L ，0011101011[()]1k k s j k y y y y G y y y ψ⎡⎤⎢⎥⎢⎥==⎢⎥⎢⎥⎢⎥⎣⎦L LM M M M L [(,)]i j f x y =U4、观察比较计算)5,,2,1,8,,2,1)(,(),,(****⋅⋅⋅=⋅⋅⋅=j i y x p y x f j i j i 的值并输出结果，以观察),(y x p 逼近),(y x f 的效果。
《数值分析（A）》计算实习题目三一、题目关于x，y，t，u，v，w的下列方程组0.5cost+u+v+w-x=2.67t+0.5sinu+v+w-y=1.070.5t+u+cosv+w-x=3.74t+0.5u+v+sinw-y=0.79以及关于z，t，u的下列二维数表确定了一个二元函数z=f(x，y)。
1．试用数值方法求出f(x，y)在区域D={(x,y)︱0≤x≤0.8,0.5≤y≤1.5}上的一个近似表达式,0(,)kr srs r s p x y c x y==∑要求p(x,y)一最小的k 值达到以下的精度102027((,)(,))10iji j i j f x yp x y σ-===-≤∑∑其中x i =0.08i ，y j =0.5+0.05j 。
2．计算f(x i *,y j *),p(x i *,y j *)(i=1,2,…,8;j=1,2,…,5)的值，以观察p(x,y)逼近f(x,y)的效果，其中x i *=0.1i,y j *=0.5+0.2j 。
二、算法设计方案1.将0.08(0,1,,10)i x i i *== 和0.50.05(0,1,,20)j y j j *=+= 代入非线性方程组中，用牛顿法解出i t 和j u ；2.以采取分片二次插值，选择（m ，n ）满足,2322m i m h h t t t m -<≤+≤≤ ,2322n j n u u u n ττ-<≤+≤≤如果12ih t t ≤+或42ih t t >-，则m=1或4；如果12ju u τ≤+或42ju u τ>-，则n=1或n=4。
选择(,)(1,,1;1,,1)k r t u km m m r n n n =-+=-+为插值节点，相应的Lagrange形式的插值多项式为),()(~)(),(111122r k r m m k n n r k u t f u l t l u t p ∑∑+-=+-==其中11()m w k w m k ww kt t l t t t +=-≠-=-∏(k=m-1, m, m+1)∏+≠-=--=11)(~n rw n w wr w r y y y y u l (r=n-1, n, n+1)并将i t 和j u 代入22(,)p t u ，便得到了数表,,(,)i j i j x y f x y 。
3.进行曲面拟和系数矩阵[]rs c *=C ，11()()TTT--=C B B B U G G G其中0011101011[()]1kk r i k x x x x x x x ϕ⎡⎤⎢⎥⎢⎥==⎢⎥⎢⎥⎢⎥⎣⎦B，0011101011[()]1kk s j k y y y y G y y y ψ⎡⎤⎢⎥⎢⎥==⎢⎥⎢⎥⎢⎥⎣⎦[(,)]i j f x y =Uk 从0逐渐增大，直到710σ-≤，便得到了要求精度的系数rs c 。
4.由前面得到的函数关系，根据重新取值的x ，y 可以分别得到新的数表，比较两组数据观察逼近效果三、全部源程#include "stdio.h"#include "stdafx.h"#include "math.h"//高斯列主元法解线性方程组void gauss(double a[4][4],double b[4],double x[4]) {intk,i,j,ik;double mik,temp,maxa,sum,tempb;for (k=0;k<3;k++){temp=fabs(a[k][k]);for(i=k;i<4;i++)if (fabs(a[i][k])>temp){ik=i;maxa=fabs(a[i][k]);}if(ik!=k){for (j=k;j<4;j++){temp=a[k][j];a[k][j]=a[ik][j];a[ik][j]=temp;}tempb=b[k];b[k]=b[ik];b[ik]=tempb;}for (i=k+1;i<4;i++){mik=a[i][k]/a[k][k];for (j=k+1;j<4;j++)a[i][j]=a[i][j]-mik*a[k][j];b[i]=b[i]-mik*b[k];}}x[3]=b[3]/a[3][3];for (k=2;k>=0;k--){sum=0;for (j=k+1;j<4;j++){sum+=a[k][j]*x[j];x[k]=(b[k]-sum)/a[k][k];}}}//求无穷范数double FindMax(double Arr[4]){inti;double MaxVal;MaxVal=fabs(Arr[0]);for (i=0;i<4;i++)if(fabs(Arr[i])>MaxVal)MaxVal=fabs(Arr[i]);return (MaxVal);}void newton(double x, double y,doubletuvw[4])//牛顿法解非线性方程组{double A[4][4],b[4],var[4]={1,2,1,2},D_var[4],MaxVar,MaxD_Var;inti;while(1){A[0][0]=-0.5*sin(var[0]);A[0][1]=1;A[0][2]=1;A[0][3]=1;A[1][0]=1;A[1][1]=0.5*cos(var[1]);A[1][2]=1;A[1][3]=1;A[2][0]=0.5;A[2][1]=1;A[2][2]=-sin(var[2]);A[2][3]=1;A[3][0]=1;A[3][1]=0.5;A[3][2]=1;A[3][3]=cos(var[3]);//系数矩阵（求导数）b[0]=-(0.5*cos(var[0])+var[1]+var[2]+var[3]-x-2.67);b[1]=-(var[0]+0.5*sin(var[1])+var[2]+var[3]-y-1.07);b[2]=-(0.5*var[0]+var[1]+cos(var[2])+var[3]-x-3.74);b[3]=-(var[0]+0.5*var[1]+var[2]+sin(var[3])-y-0.79);gauss(A,b,D_var);//高斯列主元法解线性方程组MaxVar=FindMax(var);MaxD_Var=FindMax(D_var);if ((MaxD_Var/MaxV ar)<(1e-12)) //精度判断{for (i=0;i<4;i++)tuvw[i]=var[i];break;}for (i=0;i<4;i++)var[i]=var[i]+D_var[i];}}//拉格朗日插值系数Lkdouble lagrange_t(double x,inti,int k){int m;double value=1,t[6]={0,0.2,0.4,0.6,0.8,1.0};for (m=i-1;m<=i+1;m++)if(m!=k)value=value*(x-t[m])/(t[k]-t[m]);return value;}//拉格朗日插值系数Lrdouble lagrange_u(double y,intj,int r){int m;double value=1;double u[6]={0,0.4,0.8,1.2,1.6,2};for (m=j-1;m<=j+1;m++)if(m!=r)value=value*(y-u[m])/(u[r]-u[m]);return value;}//二元二次分片代数插值double interpolation(double tt,doubleuu){inti,j,k,r;double zvalue;double z[6][6]={ {-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}};switch (int(10*tt)){case 0:case 1:case 2:i=1;break;case 3:case 4:i=2;break;case 5:case 6:i=3;break;case 7:case 8:case 9:i=4;break;default: printf("error t\n");}switch (int(5*uu)){case 0:case 1:case 2:j=1;break;case 3:case 4:j=2;break;case 5:case 6:j=3;break;case 7:case 8:case 9:j=4;break;default:printf("error u\n");}zvalue=0;for (k=i-1;k<=i+1;k++)for (r=j-1;r<=j+1;r++)zvalue+=z[k][r]*lagrange_t(tt,i,k)*lagrange_u(uu,j,r);//拉格朗日插值系数Lk,Lr return zvalue;}//初始化矩阵Bvoid init_B(double B[][10],int k){double x[11],sum;inti,j,m;for (i=0;i<11;i++)x[i]=0.08*i;for (i=0;i<11;i++){B[i][0]=1;for (j=1;j<=k;j++){sum=1;for (m=1;m<=j;m++)sum=sum*x[i];B[i][j]=sum;}}}//初始化矩阵Gvoid init_G(double G[][10],int k) {inti,j,m;double y[21],sum;for (i=0;i<21;i++)y[i]=0.5+0.05*i;for (i=0;i<21;i++){G[i][0]=1;for (j=1;j<=k;j++){sum=1;for (m=1;m<=j;m++)sum=sum*y[i];G[i][j]=sum;}}}//求逆矩阵void inverse(double A[][10],int n) {double b[10][20],c[10][10],temp;inti,j,k;for (i=0;i<n;i++)for (j=0;j<2*n;j++)if(j<n)b[i][j]=A[i][j];else b[i][j]=0;for (i=0;i<n;i++)b[i][n+i]=1;for (k=0;k<n;k++){temp=b[k][k];i=k;while (b[k][k]==0){b[k][k]=b[i+1][k];i++;}if (i>k){b[i][k]=temp;for (j=0;j<k;j++){temp=b[k][j];b[k][j]=b[i][j];b[i][j]=temp;}for (j=k+1;j<2*n;j++){temp=b[k][j];b[k][j]=b[i][j];b[i][j]=temp;}}for (i=k+1;i<n;i++)for (j=2*n-1;j>=k;j--)b[i][j]=b[i][j]-b[i][k]*b[k][j]/b[k][k];for (j=2*n-1;j>=k;j--)b[k][j]=b[k][j]/b[k][k];}for (k=n-1;k>0;k--)for (i=0;i<k;i++)for(j=2*n-1;j>=k;j--)b[i][j]=b[i][j]-b[i][k]*b[k][j];for (i=0;i<n;i++)for (j=0;j<n;j++)c[i][j]=b[i][n+j];for (i=0;i<n;i++)for (j=0;j<n;j++)A[i][j]=c[i][j];}//最小二乘拟合曲面double surface_fitting(double x,doubley,double C[][10],intk,double Z[][21]) {double B[11][10],G[21][10],BTB[10][10],GTG[10][10],U[11][21],UU[21][21]; inti,j,m;double sum;init_B(B,k);init_G(G,k);//初始化矩阵B和Gfor (i=0;i<11;i++)for (j=0;j<21;j++)U[i][j]=Z[i][j];for (i=0;i<=k;i++)for (j=0;j<=k;j++){sum=0;for(m=0;m<11;m++)sum=sum+B[m][i]*B[m][j];BTB[i][j]=sum;}for (i=0;i<=k;i++)for (j=0;j<=k;j++){sum=0;for(m=0;m<21;m++)sum=sum+G[m][i]*G[m][j];GTG[i][j]=sum;}inverse(BTB,k+1);//求逆inverse(GTG,k+1);//求逆//求矩阵BT*Ufor (i=0;i<=k;i++)for (j=0;j<21;j++){sum=0;for(m=0;m<11;m++)sum+=B[m][i]*U[m][j];UU[i][j]=sum;}//求矩阵BT*U*Gfor (i=0;i<=k;i++)for (j=0;j<=k;j++){sum=0;for(m=0;m<21;m++)sum+=UU[i][m]*G[m][j];C[i][j]=sum;}//求矩阵BIB*BT*U*Gfor (i=0;i<=k;i++)for (j=0;j<=k;j++){sum=0;for(m=0;m<=k;m++)sum+=BTB[i][m]*C[m][j];UU[i][j]=sum;}//求矩阵BIB*BT*U*G*GTGfor (i=0;i<=k;i++)for (j=0;j<=k;j++){sum=0;for(m=0;m<=k;m++)sum+=UU[i][m]*GTG[m][j];C[i][j]=sum;}//求和，得出系数Cdouble result=0,r,s;for (i=0;i<=k;i++){r=1;for (m=1;m<=i;m++)r=r*x;for (j=0;j<=k;j++){s=1;for (m=1;m<=j;m++)s=s*y;result+=C[i][j]*r*s;}}return result;}//主函数void main(){inti,j;double x[11],y[21],z[11][21],C[10][10],tuvw[4];for (i=0;i<11;i++)x[i]=0.08*i;for (j=0;j<21;j++)y[j]=0.5+0.05*j; //初始化x、yprintf("数表xi,yi,f(xi,yi)的值：\n");for (i=0;i<11;i++)for (j=0;j<21;j++){newton(x[i],y[j],tuvw); //牛顿解非线性方程组z[i][j]=interpolation(tuvw[0],tuvw[1]);//二元二次插值printf("x[%d]=%2.2f,y[%d]=%2.2f,f(x%d,y%d)=%2.12e\n",i,x[i],j,y[j],i,j,z[i][j]);}//去曲面拟合double delta=0;int k=0;printf("\n");printf("选择过程中的k、σ：\n");while(1){delta=0;for (i=0;i<11;i++)for (j=0;j<21;j++)delta+=(z[i][j]-surface_fitting(x[i],y[j],C,k,z))*(z[i][j]-surface_fitting(x[i],y[j],C,k,z));printf("k=%d,σ=%2.12e\n",k,delta);if (delta<=1e-7)break;k++;}printf("\n达到精度要求时\nk=%d,σ=%2.12e\n",k,delta);printf("Crs：\n");for (i=0;i<=k;i++){for (j=0;j<=k;j++)printf("C[%d,%d]=%2.12e\n",i,j,C[i][j]);}//输出数表x*,y*,f*,p*double xx[8],yy[5];for (i=0;i<8;i++)xx[i]=0.1*(i+1);for(j=0;j<5;j++)yy[j]=0.5+0.2*(j+1);printf("\n数表x*i,y*j,f(x*i,y*j),p(x*i,y*j)：\n");for (i=0;i<8;i++)for (j=0;j<5;j++){newton(xx[i],yy[j],tuvw);printf("x*[%d]=%2.2f,y*[%d]=%2.2f,f*(x%d,y%d)=%2.12e,p*(x%d,y%d)=%2.12e\n", i+1,xx[i],j+1,yy[j],i+1,j+1,interpolation(tuvw[0],tuvw[1]),i+1,j+1,surface_fitting(xx[i], yy[j],C,k,z));}}四、程序运行结果：数表xi,yi,f(xi,yi)的值：x[0]=0.00, y[0]=0.50, f(x0,y0)=4.465040184799e-001x[0]=0.00, y[1]=0.55, f(x0,y1)=3.246832629274e-001x[0]=0.00, y[2]=0.60, f(x0,y2)=2.101596866825e-001x[0]=0.00, y[3]=0.65, f(x0,y3)=1.030436083159e-001x[0]=0.00, y[4]=0.70, f(x0,y4)=3.401895562659e-003x[0]=0.00, y[5]=0.75, f(x0,y5)=-8.873581363801e-002x[0]=0.00, y[6]=0.80, f(x0,y6)=-1.733716327497e-001x[0]=0.00, y[7]=0.85, f(x0,y7)=-2.505346114666e-001x[0]=0.00, y[8]=0.90, f(x0,y8)=-3.202765063876e-001x[0]=0.00, y[9]=0.95, f(x0,y9)=-3.826680661097e-001x[0]=0.00, y[10]=1.00, f(x0,y10)=-4.377957667384e-001x[0]=0.00, y[11]=1.05, f(x0,y11)=-4.857589414438e-001x[0]=0.00, y[12]=1.10, f(x0,y12)=-5.266672548835e-001x[0]=0.00, y[13]=1.15, f(x0,y13)=-5.606384797965e-001x[0]=0.00, y[14]=1.20, f(x0,y14)=-5.877965387677e-001x[0]=0.00, y[15]=1.25, f(x0,y15)=-6.0826********e-001x[0]=0.00, y[16]=1.30, f(x0,y16)=-6.221894528764e-001x[0]=0.00, y[17]=1.35, f(x0,y17)=-6.296883781856e-001x[0]=0.00, y[18]=1.40, f(x0,y18)=-6.308997600028e-001x[0]=0.00, y[19]=1.45, f(x0,y19)=-6.259561525454e-001x[0]=0.00, y[20]=1.50, f(x0,y20)=-6.149885466094e-001x[1]=0.08, y[0]=0.50, f(x1,y0)=6.380152265102e-001x[1]=0.08, y[1]=0.55, f(x1,y1)=5.0661********e-001x[1]=0.08, y[2]=0.60, f(x1,y2)=3.821763692772e-001x[1]=0.08, y[3]=0.65, f(x1,y3)=2.648634911536e-001x[1]=0.08, y[4]=0.70, f(x1,y4)=1.547802002848e-001x[1]=0.08, y[5]=0.75, f(x1,y5)=5.199268349093e-002x[1]=0.08, y[6]=0.80, f(x1,y6)=-4.346804020491e-002x[1]=0.08, y[7]=0.85, f(x1,y7)=-1.316010567885e-001x[1]=0.08, y[8]=0.90, f(x1,y8)=-2.124310883088e-001x[1]=0.08, y[9]=0.95, f(x1,y9)=-2.860045510580e-001x[1]=0.08, y[10]=1.00, f(x1,y10)=-3.523860789794e-001x[1]=0.08, y[11]=1.05, f(x1,y11)=-4.116554565222e-001x[1]=0.08, y[12]=1.10, f(x1,y12)=-4.639049115188e-001x[1]=0.08, y[13]=1.15, f(x1,y13)=-5.092367247005e-001x[1]=0.08, y[14]=1.20, f(x1,y14)=-5.477611179623e-001x[1]=0.08, y[16]=1.30, f(x1,y16)=-6.048572588895e-001 x[1]=0.08, y[17]=1.35, f(x1,y17)=-6.236734213318e-001 x[1]=0.08, y[18]=1.40, f(x1,y18)=-6.361682484133e-001 x[1]=0.08, y[19]=1.45, f(x1,y19)=-6.424676566901e-001 x[1]=0.08, y[20]=1.50, f(x1,y20)=-6.426971026996e-001 x[2]=0.16, y[0]=0.50, f(x2,y0)=8.400813957652e-001x[2]=0.16, y[1]=0.55, f(x2,y1)=6.997641656726e-001x[2]=0.16, y[2]=0.60, f(x2,y2)=5.660614423514e-001x[2]=0.16, y[3]=0.65, f(x2,y3)=4.391716081175e-001x[2]=0.16, y[4]=0.70, f(x2,y4)=3.192421380408e-001x[2]=0.16, y[5]=0.75, f(x2,y5)=2.063761923874e-001x[2]=0.16, y[6]=0.80, f(x2,y6)=1.006385238914e-001x[2]=0.16, y[7]=0.85, f(x2,y7)=2.060740067835e-003x[2]=0.16, y[8]=0.90, f(x2,y8)=-8.935402476698e-002 x[2]=0.16, y[9]=0.95, f(x2,y9)=-1.736269688648e-001 x[2]=0.16, y[10]=1.00, f(x2,y10)=-2.507999561599e-001 x[2]=0.16, y[11]=1.05, f(x2,y11)=-3.209322694446e-001 x[2]=0.16, y[12]=1.10, f(x2,y12)=-3.840977350046e-001 x[2]=0.16, y[13]=1.15, f(x2,y13)=-4.403821754175e-001 x[2]=0.16, y[14]=1.20, f(x2,y14)=-4.898811523126e-001 x[2]=0.16, y[15]=1.25, f(x2,y15)=-5.326979655338e-001 x[2]=0.16, y[16]=1.30, f(x2,y16)=-5.689418792921e-001 x[2]=0.16, y[17]=1.35, f(x2,y17)=-5.987265495151e-001 x[2]=0.16, y[18]=1.40, f(x2,y18)=-6.221686297503e-001 x[2]=0.16, y[19]=1.45, f(x2,y19)=-6.393865356972e-001 x[2]=0.16, y[20]=1.50, f(x2,y20)=-6.504993507878e-001 x[3]=0.24, y[0]=0.50, f(x3,y0)=1.0515********e+000x[3]=0.24, y[1]=0.55, f(x3,y1)=9.029*********e-001x[3]=0.24, y[2]=0.60, f(x3,y2)=7.605802668593e-001x[3]=0.24, y[3]=0.65, f(x3,y3)=6.247151981455e-001x[3]=0.24, y[4]=0.70, f(x3,y4)=4.955197560009e-001x[3]=0.24, y[5]=0.75, f(x3,y5)=3.731340427746e-001x[3]=0.24, y[6]=0.80, f(x3,y6)=2.576567488723e-001x[3]=0.24, y[7]=0.85, f(x3,y7)=1.491505594102e-001x[3]=0.24, y[8]=0.90, f(x3,y8)=4.764698677337e-002x[3]=0.24, y[9]=0.95, f(x3,y9)=-4.684932320146e-002 x[3]=0.24, y[10]=1.00, f(x3,y10)=-1.343567603849e-001 x[3]=0.24, y[11]=1.05, f(x3,y11)=-2.149133449274e-001 x[3]=0.24, y[12]=1.10, f(x3,y12)=-2.885737006348e-001 x[3]=0.24, y[13]=1.15, f(x3,y13)=-3.554063647857e-001 x[3]=0.24, y[14]=1.20, f(x3,y14)=-4.154913964886e-001 x[3]=0.24, y[15]=1.25, f(x3,y15)=-4.689182499695e-001 x[3]=0.24, y[16]=1.30, f(x3,y16)=-5.157838831247e-001x[3]=0.24, y[18]=1.40, f(x3,y18)=-5.902469305629e-001 x[3]=0.24, y[19]=1.45, f(x3,y19)=-6.180615482412e-001 x[3]=0.24, y[20]=1.50, f(x3,y20)=-6.397468392579e-001 x[4]=0.32, y[0]=0.50, f(x4,y0)=1.271246751481e+000x[4]=0.32, y[1]=0.55, f(x4,y1)=1.115002018146e+000x[4]=0.32, y[2]=0.60, f(x4,y2)=9.646077272154e-001x[4]=0.32, y[3]=0.65, f(x4,y3)=8.203473694749e-001x[4]=0.32, y[4]=0.70, f(x4,y4)=6.824476781794e-001x[4]=0.32, y[5]=0.75, f(x4,y5)=5.510852085975e-001x[4]=0.32, y[6]=0.80, f(x4,y6)=4.263923859018e-001x[4]=0.32, y[7]=0.85, f(x4,y7)=3.084629956332e-001x[4]=0.32, y[8]=0.90, f(x4,y8)=1.973571296919e-001x[4]=0.32, y[9]=0.95, f(x4,y9)=9.310562085941e-002x[4]=0.32, y[10]=1.00, f(x4,y10)=-4.285992234035e-003 x[4]=0.32, y[11]=1.05, f(x4,y11)=-9.483392529689e-002 x[4]=0.32, y[12]=1.10, f(x4,y12)=-1.785729903640e-001 x[4]=0.32, y[13]=1.15, f(x4,y13)=-2.555537790546e-001 x[4]=0.32, y[14]=1.20, f(x4,y14)=-3.258401501575e-001 x[4]=0.32, y[15]=1.25, f(x4,y15)=-3.895069883634e-001 x[4]=0.32, y[16]=1.30, f(x4,y16)=-4.466382045995e-001 x[4]=0.32, y[17]=1.35, f(x4,y17)=-4.973249517677e-001 x[4]=0.32, y[18]=1.40, f(x4,y18)=-5.416640326994e-001 x[4]=0.32, y[19]=1.45, f(x4,y19)=-5.797564797951e-001 x[4]=0.32, y[20]=1.50, f(x4,y20)=-6.117062881476e-001 x[5]=0.40, y[0]=0.50, f(x5,y0)=1.498321052481e+000x[5]=0.40, y[1]=0.55, f(x5,y1)=1.334998632066e+000x[5]=0.40, y[2]=0.60, f(x5,y2)=1.177125123739e+000x[5]=0.40, y[3]=0.65, f(x5,y3)=1.025*********e+000x[5]=0.40, y[4]=0.70, f(x5,y4)=8.789600231743e-001x[5]=0.40, y[5]=0.75, f(x5,y5)=7.391451087035e-001x[5]=0.40, y[6]=0.80, f(x5,y6)=6.0574********e-001x[5]=0.40, y[7]=0.85, f(x5,y7)=4.788838610666e-001x[5]=0.40, y[8]=0.90, f(x5,y8)=3.586506258818e-001x[5]=0.40, y[9]=0.95, f(x5,y9)=2.451022361964e-001x[5]=0.40, y[10]=1.00, f(x5,y10)=1.382683509285e-001 x[5]=0.40, y[11]=1.05, f(x5,y11)=3.815486540699e-002 x[5]=0.40, y[12]=1.10, f(x5,y12)=-5.525282116814e-002 x[5]=0.40, y[13]=1.15, f(x5,y13)=-1.419868808137e-001 x[5]=0.40, y[14]=1.20, f(x5,y14)=-2.220944390959e-001 x[5]=0.40, y[15]=1.25, f(x5,y15)=-2.956352324598e-001 x[5]=0.40, y[16]=1.30, f(x5,y16)=-3.626795115028e-001 x[5]=0.40, y[17]=1.35, f(x5,y17)=-4.233061642240e-001 x[5]=0.40, y[18]=1.40, f(x5,y18)=-4.776010361325e-001x[5]=0.40, y[20]=1.50, f(x5,y20)=-5.675647436551e-001 x[6]=0.48, y[0]=0.50, f(x6,y0)=1.731892740382e+000x[6]=0.48, y[1]=0.55, f(x6,y1)=1.562034577208e+000x[6]=0.48, y[2]=0.60, f(x6,y2)=1.397216918208e+000x[6]=0.48, y[3]=0.65, f(x6,y3)=1.237801006739e+000x[6]=0.48, y[4]=0.70, f(x6,y4)=1.0840********e+000x[6]=0.48, y[5]=0.75, f(x6,y5)=9.363227723149e-001x[6]=0.48, y[6]=0.80, f(x6,y6)=7.947044490537e-001x[6]=0.48, y[7]=0.85, f(x6,y7)=6.593871980282e-001x[6]=0.48, y[8]=0.90, f(x6,y8)=5.304875868399e-001x[6]=0.48, y[9]=0.95, f(x6,y9)=4.080886854542e-001x[6]=0.48, y[10]=1.00, f(x6,y10)=2.922442012295e-001 x[6]=0.48, y[11]=1.05, f(x6,y11)=1.829822068536e-001 x[6]=0.48, y[12]=1.10, f(x6,y12)=8.030849403542e-002 x[6]=0.48, y[13]=1.15, f(x6,y13)=-1.579041305164e-002 x[6]=0.48, y[14]=1.20, f(x6,y14)=-1.0534********e-001 x[6]=0.48, y[15]=1.25, f(x6,y15)=-1.883980906096e-001 x[6]=0.48, y[16]=1.30, f(x6,y16)=-2.650071493189e-001 x[6]=0.48, y[17]=1.35, f(x6,y17)=-3.352378389040e-001 x[6]=0.48, y[18]=1.40, f(x6,y18)=-3.991645038868e-001 x[6]=0.48, y[19]=1.45, f(x6,y19)=-4.568681433016e-001 x[6]=0.48, y[20]=1.50, f(x6,y20)=-5.084349932782e-001 x[7]=0.56, y[0]=0.50, f(x7,y0)=1.971221786400e+000x[7]=0.56, y[1]=0.55, f(x7,y1)=1.795329599501e+000x[7]=0.56, y[2]=0.60, f(x7,y2)=1.624067113228e+000x[7]=0.56, y[3]=0.65, f(x7,y3)=1.457830582708e+000x[7]=0.56, y[4]=0.70, f(x7,y4)=1.296954649752e+000x[7]=0.56, y[5]=0.75, f(x7,y5)=1.141718105447e+000x[7]=0.56, y[6]=0.80, f(x7,y6)=9.923495333243e-001x[7]=0.56, y[7]=0.85, f(x7,y7)=8.490326633294e-001x[7]=0.56, y[8]=0.90, f(x7,y8)=7.119113522641e-001x[7]=0.56, y[9]=0.95, f(x7,y9)=5.810941589219e-001x[7]=0.56, y[10]=1.00, f(x7,y10)=4.566585132334e-001 x[7]=0.56, y[11]=1.05, f(x7,y11)=3.386544961394e-001 x[7]=0.56, y[12]=1.10, f(x7,y12)=2.271082557696e-001 x[7]=0.56, y[13]=1.15, f(x7,y13)=1.220250891932e-001 x[7]=0.56, y[14]=1.20, f(x7,y14)=2.339221963760e-002 x[7]=0.56, y[15]=1.25, f(x7,y15)=-6.881870197104e-002 x[7]=0.56, y[16]=1.30, f(x7,y16)=-1.546493442129e-001 x[7]=0.56, y[17]=1.35, f(x7,y17)=-2.341526664587e-001 x[7]=0.56, y[18]=1.40, f(x7,y18)=-3.0739********e-001 x[7]=0.56, y[19]=1.45, f(x7,y19)=-3.744348623481e-001 x[7]=0.56, y[20]=1.50, f(x7,y20)=-4.353605565359e-001x[8]=0.64, y[1]=0.55, f(x8,y1)=2.034201133607e+000x[8]=0.64, y[2]=0.60, f(x8,y2)=1.856955143619e+000x[8]=0.64, y[3]=0.65, f(x8,y3)=1.684358164161e+000x[8]=0.64, y[4]=0.70, f(x8,y4)=1.516776352400e+000x[8]=0.64, y[5]=0.75, f(x8,y5)=1.354519041151e+000x[8]=0.64, y[6]=0.80, f(x8,y6)=1.197844086673e+000x[8]=0.64, y[7]=0.85, f(x8,y7)=1.046963049419e+000x[8]=0.64, y[8]=0.90, f(x8,y8)=9.020*********e-001x[8]=0.64, y[9]=0.95, f(x8,y9)=7.632264776629e-001x[8]=0.64, y[10]=1.00, f(x8,y10)=6.306048219543e-001 x[8]=0.64, y[11]=1.05, f(x8,y11)=5.042528145972e-001 x[8]=0.64, y[12]=1.10, f(x8,y12)=3.842167155457e-001 x[8]=0.64, y[13]=1.15, f(x8,y13)=2.705204766410e-001 x[8]=0.64, y[14]=1.20, f(x8,y14)=1.631685723996e-001 x[8]=0.64, y[15]=1.25, f(x8,y15)=6.214855811676e-002 x[8]=0.64, y[16]=1.30, f(x8,y16)=-3.256661939682e-002 x[8]=0.64, y[17]=1.35, f(x8,y17)=-1.210165348444e-001 x[8]=0.64, y[18]=1.40, f(x8,y18)=-2.032513996228e-001 x[8]=0.64, y[19]=1.45, f(x8,y19)=-2.793303595584e-001 x[8]=0.64, y[20]=1.50, f(x8,y20)=-3.493199575400e-001 x[9]=0.72, y[0]=0.50, f(x9,y0)=2.464684222660e+000x[9]=0.72, y[1]=0.55, f(x9,y1)=2.278058979399e+000x[9]=0.72, y[2]=0.60, f(x9,y2)=2.0952********e+000x[9]=0.72, y[3]=0.65, f(x9,y3)=1.916718127997e+000x[9]=0.72, y[4]=0.70, f(x9,y4)=1.742854628776e+000x[9]=0.72, y[5]=0.75, f(x9,y5)=1.573998427334e+000x[9]=0.72, y[6]=0.80, f(x9,y6)=1.410434835231e+000x[9]=0.72, y[7]=0.85, f(x9,y7)=1.252401750608e+000x[9]=0.72, y[8]=0.90, f(x9,y8)=1.100094409628e+000x[9]=0.72, y[9]=0.95, f(x9,y9)=9.536698512613e-001x[9]=0.72, y[10]=1.00, f(x9,y10)=8.132510552489e-001 x[9]=0.72, y[11]=1.05, f(x9,y11)=6.789307429659e-001 x[9]=0.72, y[12]=1.10, f(x9,y12)=5.507748485043e-001 x[9]=0.72, y[13]=1.15, f(x9,y13)=4.288256769731e-001 x[9]=0.72, y[14]=1.20, f(x9,y14)=3.131047717398e-001 x[9]=0.72, y[15]=1.25, f(x9,y15)=2.036155140327e-001 x[9]=0.72, y[16]=1.30, f(x9,y16)=1.003454782409e-001 x[9]=0.72, y[17]=1.35, f(x9,y17)=3.268565186572e-003 x[9]=0.72, y[18]=1.40, f(x9,y18)=-8.765306591329e-002 x[9]=0.72, y[19]=1.45, f(x9,y19)=-1.724672478188e-001 x[9]=0.72, y[20]=1.50, f(x9,y20)=-2.512302207523e-001 x[10]=0.80, y[0]=0.50, f(x10,y0)=2.717811109469e+000 x[10]=0.80, y[1]=0.55, f(x10,y1)=2.526399501256e+000x[10]=0.80, y[3]=0.65, f(x10,y3)=2.154329377280e+000x[10]=0.80, y[4]=0.70, f(x10,y4)=1.974574556652e+000x[10]=0.80, y[5]=0.75, f(x10,y5)=1.799510579099e+000x[10]=0.80, y[6]=0.80, f(x10,y6)=1.629448220554e+000x[10]=0.80, y[7]=0.85, f(x10,y7)=1.464650043751e+000x[10]=0.80, y[8]=0.90, f(x10,y8)=1.305334967651e+000x[10]=0.80, y[9]=0.95, f(x10,y9)=1.151682621307e+000x[10]=0.80, y[10]=1.00, f(x10,y10)=1.003837419906e+000 x[10]=0.80, y[11]=1.05, f(x10,y11)=8.619123372279e-001 x[10]=0.80, y[12]=1.10, f(x10,y12)=7.259923711112e-001 x[10]=0.80, y[13]=1.15, f(x10,y13)=5.961377115201e-001 x[10]=0.80, y[14]=1.20, f(x10,y14)=4.723866279136e-001 x[10]=0.80, y[15]=1.25, f(x10,y15)=3.547580958979e-001 x[10]=0.80, y[16]=1.30, f(x10,y16)=2.432541841813e-001 x[10]=0.80, y[17]=1.35, f(x10,y17)=1.378622225247e-001 x[10]=0.80, y[18]=1.40, f(x10,y18)=3.855677032640e-002 x[10]=0.80, y[19]=1.45, f(x10,y19)=-5.469859593446e-002 x[10]=0.80, y[20]=1.50, f(x10,y20)=-1.419496597088e-001选择过程中的k、σ：k=0,σ=1.442880771836e+002k=1,σ=3.220908973634e+000k=2,σ=4.659960033243e-003k=3,σ=1.721175379293e-004k=4,σ=3.309534300922e-006k=5,σ=2.541512318326e-008达到精度要求时k=5,σ=2.541512318326e-008Crs：C[0,0]=2.021*********e+000C[0,1]=-3.668423846364e+000C[0,2]=7.092421418056e-001C[0,3]=8.486122535542e-001C[0,4]=-4.159009451978e-001C[0,5]=6.743269192521e-002C[1,0]=3.191954374357e+000C[1,1]=-7.413689485984e-001C[1,2]=-2.696557302843e+000C[1,3]=1.630583174760e+000C[1,4]=-4.844127955148e-001C[1,5]=6.0553********e-002C[2,0]=2.564098287839e-001C[2,1]=1.582655804232e+000C[2,2]=-4.694162700325e-001C[2,3]=-7.498418726027e-002C[2,4]=9.883782174438e-002C[2,5]=-2.036812738515e-002C[3,0]=-2.676476564375e-001C[3,1]=-7.394417682663e-001C[3,2]=1.096321268007e+000C[3,3]=-8.283654926345e-001C[3,4]=3.138025505468e-001C[3,5]=-4.814400518080e-002C[4,0]=2.151515221121e-001C[4,1]=-1.652099197963e-001C[4,2]=-1.012971356977e-001C[4,3]=2.739136226010e-001C[4,4]=-1.569928317913e-001C[4,5]=2.962164662313e-002C[5,0]=-5.501011264460e-002C[5,1]=1.381011343328e-001C[5,2]=-1.250693502516e-001C[5,3]=2.885186714411e-002C[5,4]=9.856276497885e-003C[5,5]=-3.877033908793e-003数表x*i,y*j,f(x*i,y*j),p(x*i,y*j)：x*[1]=0.10, y*[1]=0.70, f*(x1,y1)=1.947204079177e-001, p*(x1,y1)=1.947303552407e-001 x*[1]=0.10, y*[2]=0.90, f*(x1,y2)=-1.830370791887e-001, p*(x1,y2)=-1.830418407586e-001 x*[1]=0.10, y*[3]=1.10, f*(x1,y3)=-4.454976469148e-001, p*(x1,y3)=-4.455000487616e-001 x*[1]=0.10, y*[4]=1.30, f*(x1,y4)=-5.975667076413e-001, p*(x1,y4)=-5.975588666853e-001 x*[1]=0.10, y*[5]=1.50, f*(x1,y5)=-6.464595939011e-001, p*(x1,y5)=-6.464461277567e-001 x*[2]=0.20, y*[1]=0.70, f*(x2,y1)=4.0597********e-001, p*(x2,y1)=4.0598********e-001 x*[2]=0.20, y*[2]=0.90, f*(x2,y2)=-2.251595837462e-002, p*(x2,y2)=-2.252111976633e-002 x*[2]=0.20, y*[3]=1.10, f*(x2,y3)=-3.382208160396e-001, p*(x2,y3)=-3.382240227895e-001 x*[2]=0.20, y*[4]=1.30, f*(x2,y4)=-5.444378315219e-001, p*(x2,y4)=-5.444304561199e-001 x*[2]=0.20, y*[5]=1.50, f*(x2,y5)=-6.473613385679e-001, p*(x2,y5)=-6.473480147123e-001 x*[3]=0.30, y*[1]=0.70, f*(x3,y1)=6.347771951509e-001, p*(x3,y1)=6.347874560151e-001 x*[3]=0.30, y*[2]=0.90, f*(x3,y2)=1.588011688394e-001, p*(x3,y2)=1.587962910079e-001 x*[3]=0.30, y*[3]=1.10, f*(x3,y3)=-2.0736********e-001, p*(x3,y3)=-2.0736********e-001 x*[3]=0.30, y*[4]=1.30, f*(x3,y4)=-4.653579068978e-001, p*(x3,y4)=-4.653499224584e-001 x*[3]=0.30, y*[5]=1.50, f*(x3,y5)=-6.202709530749e-001, p*(x3,y5)=-6.202571262457e-001 x*[4]=0.40, y*[1]=0.70, f*(x4,y1)=8.789600231743e-001, p*(x4,y1)=8.789698723457e-001 x*[4]=0.40, y*[2]=0.90, f*(x4,y2)=3.586506258818e-001, p*(x4,y2)=3.586460358660e-001 x*[4]=0.40, y*[3]=1.10, f*(x4,y3)=-5.525282116814e-002, p*(x4,y3)=-5.525541795963e-002 x*[4]=0.40, y*[4]=1.30, f*(x4,y4)=-3.626795115028e-001, p*(x4,y4)=-3.626710548561e-001x*[4]=0.40, y*[5]=1.50, f*(x4,y5)=-5.675647436551e-001, p*(x4,y5)=-5.675505482659e-001 x*[5]=0.50, y*[1]=0.70, f*(x5,y1)=1.136610910158e+000, p*(x5,y1)=1.136620368290e+000 x*[5]=0.50, y*[2]=0.90, f*(x5,y2)=5.749803409475e-001, p*(x5,y2)=5.749758252283e-001 x*[5]=0.50, y*[3]=1.10, f*(x5,y3)=1.159923767920e-001, p*(x5,y3)=1.159893591603e-001 x*[5]=0.50, y*[4]=1.30, f*(x5,y4)=-2.385683040123e-001, p*(x5,y4)=-2.385604036565e-001 x*[5]=0.50, y*[5]=1.50, f*(x5,y5)=-4.914343936557e-001, p*(x5,y5)=-4.914208309330e-001 x*[6]=0.60, y*[1]=0.70, f*(x6,y1)=1.406041798905e+000, p*(x6,y1)=1.406050718527e+000 x*[6]=0.60, y*[2]=0.90, f*(x6,y2)=8.0594********e-001, p*(x6,y2)=8.0593********e-001 x*[6]=0.60, y*[3]=1.10, f*(x6,y3)=3.044292210453e-001, p*(x6,y3)=3.044259035644e-001 x*[6]=0.60, y*[4]=1.30, f*(x6,y4)=-9.501613009962e-002, p*(x6,y4)=-9.500892391062e-002 x*[6]=0.60, y*[5]=1.50, f*(x6,y5)=-3.939023077456e-001, p*(x6,y5)=-3.938897085270e-001 x*[7]=0.70, y*[1]=0.70, f*(x7,y1)=1.685783515309e+000, p*(x7,y1)=1.685791278477e+000 x*[7]=0.70, y*[2]=0.90, f*(x7,y2)=1.049881153064e+000, p*(x7,y2)=1.049877642456e+000 x*[7]=0.70, y*[3]=1.10, f*(x7,y3)=5.082937839397e-001, p*(x7,y3)=5.082911724749e-001 x*[7]=0.70, y*[4]=1.30, f*(x7,y4)=6.614879670648e-002, p*(x7,y4)=6.615638110524e-002 x*[7]=0.70, y*[5]=1.50, f*(x7,y5)=-2.768343417776e-001, p*(x7,y5)=-2.768218193171e-001 x*[8]=0.80, y*[1]=0.70, f*(x8,y1)=1.974574556652e+000, p*(x8,y1)=1.974581369620e+000 x*[8]=0.80, y*[2]=0.90, f*(x8,y2)=1.305334967651e+000, p*(x8,y2)=1.305331820825e+000 x*[8]=0.80, y*[3]=1.10, f*(x8,y3)=7.259923711112e-001, p*(x8,y3)=7.259895225092e-001 x*[8]=0.80, y*[4]=1.30, f*(x8,y4)=2.432541841813e-001, p*(x8,y4)=2.432608163707e-001 x*[8]=0.80, y*[5]=1.50, f*(x8,y5)=-1.419496597088e-001, p*(x8,y5)=-1.419384254301e-001 Press any key to continue。
