数值分析大作业
北航数值分析大作业 第一题 幂法与反幂法
数 值 分 析（B ） 大 作 业（一）姓名： 学号： 电话：1、算法设计：①求1λ、501λ和s λ的值：s λ：s λ表示矩阵的按模最小特征值，为求得s λ直接对待求矩阵A 应用反幂法即可。
1λ、501λ：若矩阵A 的特征值满足关系 1n λλ<<且1n λλ≠，要求1λ、及501λ时，可按如下方法求解： a ． 对矩阵A 用幂法，求得按模最大的特征值1m λ。
b ． 按平移量1m λ对矩阵A 进行原点平移得矩阵1m BA I λ=+，对矩阵B 用反幂法求得B 的按模最小特征值2m λ。
c ． 321m m m λλλ=-则：113min(,)m m λλλ=，13max(,)n m m λλλ=即为所求。
②求和A 的与数5011140k k λλμλ-=+最接近的特征值ik λ（k=0，1，…39）：求矩阵A 的特征值中与P 最接近的特征值的大小，采用原点平移的方法：先求矩阵 B=A-PI 对应的按模最小特征值k β，则k β+P 即为矩阵A 与P 最接近的特征值。
在本次计算实习中则是先求平移矩阵k B A I μ=-，对该矩阵应用反幂法求得s λ，则与k μ最接近的A 的特征值为：s P λ+重复以上过程39次即可求得ik λ（k=0，1，…39）的值。
③求A 的（谱范数）条件数2cond()A 和行列式det A ：在（1）中用反幂法求矩阵A 的按模最小特征值时，要用到Doolittle 分解方法，在Doolittle 分解完成后得到的两个矩阵分别为L 和U ，则A 的行列式可由U 阵求出，即：det(A)=det(U)。
求得det(A)不为0，因此A 为非奇异的实对称矩阵，则： max 2()scond A λλ=，max λ和s λ分别为模最大特征值与模最小特征值。
2、程序源代码：#include "Stdio.h"#include "Conio.h"#include "math.h"//****************************************************************************// // 在存储带状矩阵时，下面的几个量在程序中反复用到，为方便编程故把它们定义成宏.// // M ：转换后的矩阵的行数，M=R+S+1。
数值分析大作业一
数值分析大作业一一、算法设计方案1、求λ1和λ501的值：思路：采用幂法求出按模最大特征值λmax，该值必为λ1或λ501，若λmax小于0，则λmax=λ1；否则λmax=λ501。
再经过原点平移，使用幂法迭代出矩阵A-λmax I的特征值，此时求出的按模最大特征值即为λ1和λ501的另一个值。
2、求λs的值：采用反幂法求出按模最小的特征值λmin即为λs，其中的方程组采用LU分解法进行求解。
3、求与μk最接近的特征值：对矩阵A采用带原点平移的反幂法求解最小特征值，其中平移量为：μk。
4、A的条件数cond(A)=| λmax/λmin|；5、A的行列式的值：先将A进行LU分解，再求U矩阵对角元素的乘积即为A 行列式的值。
二、源程序#include<iostream>#include<iomanip>#include<math.h>#define N 501#define E 1.0e-12 //定义精度常量#define r 2#define s 2using namespace std;double a[N];double cc[5][N];void init();double mifa();double fmifa();int max(int aa,int bb);int min(int aa,int bb);int max_3(int aa,int bb,int cc);void LU();void main(){double a1,a2,d1,d501=0,ds,det=1,miu[39],lamta,cond;int i,k;init();/*************求λ1和λ501********************/a1=mifa();if(a1<0)d1=a1; //若小于0则表示λ1的值elsed501=a1; //若大于0则表示λ501的值for(i=0;i<N;i++)a[i]=a[i]-a1;a2=mifa()+a1;if(a2<0)d1=a2; //若小于0则表示λ1的值elsed501=a2; //若大于0则表示λ501的值cout<<"λ1="<<setiosflags(ios::scientific)<<setprecision(12)<<d1<<"\t";cout<<"λ501="<<setiosflags(ios::scientific)<<setprecision(12)<<d501<<endl;/**************求λs*****************/init();ds=fmifa();cout<<"λs="<<setiosflags(ios::scientific)<<setprecision(12)<<ds<<endl;/**************求与μk最接近的特征值λik**************/cout<<"与μk最接近的特征值λik:"<<endl;for(k=0;k<39;k++){miu[k]=d1+(k+1)*(d501-d1)/40;init();for(i=0;i<N;i++)a[i]=a[i]-miu[k];lamta=fmifa()+miu[k];cout<<"λi"<<k+1<<"\t\t"<<setiosflags(ios::scientific)<<setprecision(12)<<lamta<<en dl;}/**************求A的条件数**************/cout<<"矩阵A的条件式";cond=abs(max(abs(d1),abs(d501))/ds);cout<<"cond="<<setiosflags(ios::scientific)<<setprecision(12)<<cond<<endl;/**************求A的行列式**************/cout<<"矩阵A的行列式";init();LU();for(i=0;i<N;i++){det*=cc[2][i];}cout<<"det="<<setiosflags(ios::scientific)<<setprecision(12)<<det<<endl;system("pause");}/**************初始化函数，给a[N]赋值*************/void init(){int i;for(i=1;i<=501;i++)a[i-1]=(1.64-0.024*i)*sin((double)(0.2*i))-0.64*exp((double)(0.1/i)); }/**************幂法求最大绝对特征值**************/double mifa(){int i,k=0;double u[N],y[N]={0},b=0.16,c=-0.064,Beta_=0,error;for(i=0;i<501;i++)u[i]=1; //令u[N]=1for(k=1;k<2000;k++) //控制最大迭代次数为2000{/***求y(k-1)***/double sum_u=0,gh_sum_u;for(i=0;i<N;i++){sum_u+=u[i]*u[i]; }gh_sum_u=sqrt(sum_u);for(i=0;i<N;i++){y[i]=u[i]/gh_sum_u;}/****求新的uk****/u[0]=a[0]*y[0]+b*y[1]+c*y[2];u[1]=b*y[0]+a[1]*y[1]+b*y[2]+c*y[3]; //前两列和最后两列单独拿出来求中D间的循环求for(i=2;i<N-2;i++){u[i]=c*y[i-2]+b*y[i-1]+a[i]*y[i]+b*y[i+1]+c*y[i+2];}u[N-2]=c*y[N-4]+b*y[N-3]+a[N-2]*y[N-2]+b*y[N-1];u[N-1]=c*y[N-3]+b*y[N-2]+a[N-1]*y[N-1];/***求beta***/double Beta=0;for(i=0;i<N;i++){Beta+=y[i]*u[i];}//cout<<"Beta"<<k<<"="<<Beta<<"\t"; 输出每次迭代的beta /***求误差***/error=abs(Beta-Beta_)/abs(Beta);if(error<=E) //若迭代误差在精度水平内则可以停止迭代{return Beta;} //控制显示位数Beta_=Beta; //第个eta的值都要保存下来，为了与后个值进行误差计算 }if(k==2000){cout<<"error"<<endl;return 0;} //若在最大迭代次数范围内都不能满足精度要求说明不收敛}/**************反幂法求最小绝对特¬征值**************/double fmifa(){int i,k,t;double u[N],y[N]={0},yy[N]={0},b=0.16,c=-0.064,Beta_=0,error;for(i=0;i<501;i++)u[i]=1; //令u[N]=1for(k=1;k<2000;k++){double sum_u=0,gh_sum_u;for(i=0;i<N;i++){sum_u+=u[i]*u[i]; }gh_sum_u=sqrt(sum_u);for(i=0;i<N;i++){y[i]=u[i]/gh_sum_u;yy[i]=y[i]; //用重新赋值，避免求解方程组的时候改变y的值}/****LU分解法解方程组Au=y，求新的***/LU();for(i=2;i<=N;i++){double temp_b=0;for(t=max(1,i-r);t<=i-1;t++)temp_b+=cc[i-t+s][t-1]*yy[t-1];yy[i-1]=yy[i-1]-temp_b;}u[N-1]=yy[N-1]/cc[s][N-1];for(i=N-1;i>=1;i--){double temp_u=0;for(t=i+1;t<=min(i+s,N);t++)temp_u+=cc[i-t+s][t-1]*u[t-1];u[i-1]=(yy[i-1]-temp_u)/cc[s][i-1];}double Beta=0;for(i=0;i<N;i++){Beta+=y[i]*u[i];}error=abs(Beta-Beta_)/abs(Beta);if(error<=E){return (1/Beta);}Beta_=Beta;}if(k==2000){cout<<"error"<<endl;return 0;} }/**************求两数最大值的子程序**************/int max(int aa,int bb){return(aa>bb?aa:bb);}/**************求两数最小值的子程序**************/int min(int aa,int bb){return(aa<bb?aa:bb);}/**************求三数最大值的子程序**************/int max_3(int aa,int bb,int cc){ int tt;if(aa>bb)tt=aa;else tt=bb;if(tt<cc) tt=cc;return(tt);}/**************LU分解**************/void LU(){int i,j,k,t;double b=0.16,c=-0.064;/**赋值压缩后矩阵cc[5][501]**/for(i=2;i<N;i++)cc[0][i]=c;for(i=1;i<N;i++)cc[1][i]=b;for(i=0;i<N;i++)cc[2][i]=a[i];for(i=0;i<N-1;i++)cc[3][i]=b;for(i=0;i<N-2;i++)cc[4][i]=c;for(k=1;k<=N;k++){for(j=k;j<=min(k+s,N);j++){double temp=0;for(t=max_3(1,k-r,j-s);t<=k-1;t++)temp+=cc[k-t+s][t-1]*cc[t-j+s][j-1];cc[k-j+s][j-1]=cc[k-j+s][j-1]-temp;}//if(k<500){for(i=k+1;i<=min(k+r,N);i++){double temp2=0;for(t=max_3(1,i-r,k-s);t<=k-1;t++)temp2+=cc[i-t+s][t-1]*cc[t-k+s][k-1];cc[i-k+s][k-1]=(cc[i-k+s][k-1]-temp2)/cc[s][k-1];}}}}三、程序结果。
数值分析大作业
数值分析上机作业（一）一、算法的设计方案1、幂法求解λ1、λ501幂法主要用于计算矩阵的按模最大的特征值和相应的特征向量，即对于|λ1|≥|λ2|≥.....≥|λn|可以采用幂法直接求出λ1，但在本题中λ1≤λ2≤……≤λ501，我们无法判断按模最大的特征值。
但是由矩阵A的特征值条件可知|λ1|和|λ501|之间必然有一个是最大的，通过对矩阵A使用幂法迭代一定次数后得到满足精度ε=10−12的特征值λ0，然后在对矩阵A做如下的平移：B=A-λ0I由线性代数(A-PI)x=(λ-p)x可得矩阵B的特征值为：λ1-λ0、λ2-λ0…….λ501-λ0。
对B矩阵采用幂法求出B矩阵按模最大的特征值为λ∗=λ501-λ0，所以λ501=λ∗+λ0，比较λ0与λ501的大小，若λ0>λ501则λ1=λ501，λ501=λ0；若λ0<λ501，则令t=λ501，λ1=λ0，λ501=t。
求矩阵M按模最大的特征值λ的具体算法如下：任取非零向量u0∈R nηk−1=u T(k−1)∗u k−1y k−1=u k−1ηk−1u k=Ay k−1βk=y Tk−1u k(k=1,2,3……)当|βk−βk−1||βk|≤ε=10−12时，迭终终止，并且令λ1=βk2、反幂法计算λs和λik由已知条件可知λs是矩阵A 按模最小的特征值，可以应用反幂法直接求解出λs。
使用带偏移量的反幂法求解λik，其中偏移量为μk=λ1+kλ501−λ140(k=1,2,3…39)，构造矩阵C=A-μk I，矩阵C的特征值为λik−μk，对矩阵C使用反幂法求得按模最小特征值λ0，则有λik=1λ0+μk。
求解矩阵M按模最小特征值的具体算法如下：任取非零向量u 0∈R n ηk−1= u T (k−1)∗u k−1y k−1=u k−1ηk−1 Au k =y k−1βk =y T k−1u k (k=1,2,3……)在反幂法中每一次迭代都要求解线性方程组Au k =y k−1，当K 足够大时，取λn =1βk 。
北航数值分析全部三次大作业
北航数值分析全部三次大作业第一次大作业是关于解线性方程组的数值方法。
我们被要求实现各种常用的线性方程组求解算法，例如高斯消元法、LU分解法和迭代法等。
我首先学习了这些算法的原理和实现方法，并借助Python编程语言编写了这些算法的代码。
在实验中，我们使用了不同规模和条件的线性方程组进行测试，并比较了不同算法的性能和精度。
通过这个作业，我深入了解了线性方程组求解的原理和方法，提高了我的编程和数值计算能力。
第二次大作业是关于数值积分的方法。
数值积分是数值分析中的重要内容，它可以用于计算曲线的长度、函数的面积以及求解微分方程等问题。
在这个作业中，我们需要实现不同的数值积分算法，例如矩形法、梯形法和辛普森法等。
我学习了这些算法的原理和实现方法，并使用Python编写了它们的代码。
在实验中，我们计算了不同函数的积分值，并对比了不同算法的精度和效率。
通过这个作业，我深入了解了数值积分的原理和方法，提高了我的编程和数学建模能力。
第三次大作业是关于常微分方程的数值解法。
常微分方程是数值分析中的核心内容之一，它可以用于描述众多物理、化学和生物现象。
在这个作业中，我们需要实现不同的常微分方程求解算法，例如欧拉法、龙格-库塔法和Adams法等。
我学习了这些算法的原理和实现方法，并使用Python编写了它们的代码。
在实验中，我们解决了一些具体的常微分方程问题，并比较了不同算法的精度和效率。
通过这个作业，我深入了解了常微分方程的原理和方法，提高了我的编程和问题求解能力。
总的来说，北航数值分析课程的三次大作业非常有挑战性，但也非常有意义。
通过这些作业，我在数值计算和编程方面得到了很大的提升，也更加深入地了解了数值分析的理论和方法。
虽然这些作业需要大量的时间和精力，但我相信这些努力将会对我未来的学习和工作产生积极的影响。
北航数值分析报告大作业第三题(fortran)
北航数值分析报告大作业第三题(fortran)“数值分析“计算实习大作业第三题——SY1415215孔维鹏一、计算说明1、将x i=0.08i，y j=0.5+0.05j分别代入方程组（A.3）得到关于t，u，v，w的的方程组，调用离散牛顿迭代子函数求出与x i，y j对应的t i，u j。
2、调用分片二次代数插值子函数在点(t i,u j)处插值得到z(x i,y j)=f(x i,y j)，得到数表(x i,y j,f(x i,y j))。
3、对于k=1,2,3,4?，分别调用最小二乘拟合子函数计算系数矩阵c rs 及误差σ，直到满足精度，即求得最小的k值及系数矩阵c rs。
4、将x i?=0.1i，y j?=0.5+0.2j分别代入方程组（A.3）得到关于t?，u?，v?，w?的的方程组，调用离散牛顿迭代子函数求出与x i?，y j?对应的t i?，u j?，调用分片二次代数插值子函数在点(t i?,u j?)处插值得到z?(x i?,y j?)=f(x i?,y j?)；调用步骤3中求得的系数矩阵c rs求得p(x i?,y j?)，打印数表(x i?,y j?,f(x i?,y j?),p(x i?,y j?))。
二、源程序（FORTRAN）PROGRAM SY1415215DIMENSIONX(11),Y(21),T(6),U(6),Z(6,6),UX(11,21),TY(11,21),FXY(11,21), C(6,6) DIMENSIONX1(8),Y1(5),FXY1(8,5),PXY1(8,5),UX1(8,5),TY1(8,5)REAL(8) X,Y,T,U,Z,FXY,UX,TY,C,E,X1,Y1,FXY1,PXY1,UX1,TY1OPEN (1,FILE='第三题计算结果.TXT')DO I=1,11X(I)=0.08*(I-1)ENDDODO I=1,21Y(I)=0.5+0.05*(I-1)ENDDO!*****求解非线性方程组，得到z=f（t，u）的函数*******DO I=1,11DO J=1,21CALL DISNEWTON_NONLINEAR(X(I),Y(J),UX(I,J),TY(I,J)) ENDDO ENDDO!*************分片二次插值得到z=f（x，y）***********DO I=1,11DO J=1,21CALL INTERPOLATION(UX(I,J),TY(I,J),FXY(I,J))ENDDO ENDDOWRITE (1,"('数表(x,y,f(x,y)):')")WRITE (1,"(3X,'X',7X,'Y',10X,'F(X,Y)')")DO I=1,11DO J=1,21WRITE(1,'(1X,F5.2,2X,F5.3,2X,E20.13)') X(I),Y(J),FXY(I,J) ENDDOWRITE (1,"('')")ENDDO!***********最小二乘拟合得到P（x，y）**************N=11M=21WRITE (1,'(" ","K和σ分别为：")')DO K=1,20CALL LSFITTING(X,Y,FXY,C,N,M,K,K,E) WRITE (1,'(I3,2X,E20.13)') K-1,EIF(ETA).OR.(A(L,K)==TA)) THENTA=A(L,K)TL=LDO J=K,NT(K,J)=A(K,J)A(K,J)=A(TL,J)A(TL,J)=T(K,J)ENDDOTB(K)=B(K)B(K)=B(TL)B(TL)=TB(K)ENDIF ENDDODO I=K+1,NM(I,K)=A(I,K)/A(K,K)A(I,K)=0DO J=K+1,NA(I,J)=A(I,J)-M(I,K)*A(K,J) ENDDOB(I)=B(I)-M(I,K)*B(K)ENDDOENDDO!回代过程X(N)=B(N)/A(N,N)DO K=N-1,1,-1S=0.0DO J=K+1,NS=S+A(K,J)*X(J)ENDDOX(K)=(B(K)-S)/A(K,K)ENDDORETURNEND!***********求向量的无穷数************ SUBROUTINE NORM(X,N,A) DIMENSION X(N)REAL(8) X,AA=ABS(X(1))DO I=2,NIF(ABS(X(I))>ABS(X(I-1))) THENA=ABS(X(I)) ENDIFENDDORETURNEND!**************分片二次代数插值************** SUBROUTINE INTERPOLATION(U,V,W) PARAMETER (N=6,M=6)DIMENSION X(N),Y(M),Z(M,N),LK(3),LR(3)REAL(8) X,Y,Z,H,TREAL(8) U,V,W,LK,LR !U,V分别为插值点处的坐标，W为插值结果INTEGER R!**********************数据赋值********************** DATA Y/0.0,0.2,0.4,0.6,0.8,1.0/DATA X/0.0,0.4,0.8,1.2,1.6,2.0/DATA Z/-0.5,-0.42,-0.18,0.22,0.78,1.5,&&-0.34,-0.5,-0.5,-0.34,-0.02,0.46,&&0.14,-0.26,-0.5,-0.58,-0.5,-0.26,&&0.94,0.3,-0.18,-0.5,-0.66,-0.66,&&2.06,1.18,0.46,-0.1,-0.5,-0.74,&&3.5,2.38,1.42,0.62,-0.02,-0.5/H=0.4T=0.2!******************计算K,R************************* IF(UX(N-1)-H/2) THENK=N-1ELSEDO I=3,N-2IF((U>X(I)-H/2).AND.(UY(M-1)-T/2) THENR=M-1 ELSEDO J=3,M-2IF((V>Y(J)-T/2).AND.(VN) P=N IF(P>20) P=20IF(Q>M) Q=MIF(Q>20) Q=20XX=0YY=0D1=NAPX(1)=0.0DO I=1,NAPX(1)=APX(1)+X(I)ENDDOAPX(1)=APX(1)/D1DO J=1,MV(1,J)=0.0DO I=1,NV(1,J)=V(1,J)+Z(I,J)ENDDOV(1,J)=V(1,J)/D1ENDDOIF(P>1) THEND2=0.0APX(2)=0.0DO I=1,NG=X(I)-APX(1)D2=D2+G*GAPX(2)=APX(2)+(X(I)-XX)*G*G ENDDO APX(2)=APX(2)/D2BX(2)=D2/D1DO J=1,MV(2,J)=0.0DO I=1,NG=X(I)-APX(1)V(2,J)=V(2,J)+Z(I,J)*G ENDDOV(2,J)=V(2,J)/D2ENDDOD1=D2ENDIFDO K=3,PD2=0.0APX(K)=0.0DO J=1,MV(K,J)=0.0ENDDODO I=1,NG1=1.0G2=X(I)-APX(1)DO J=3,KG=(X(I)-APX(J-1))*G2-BX(J-1)*G1 G1=G2 G2=GENDDOD2=D2+G*GAPX(K)=APX(K)+X(I)*G*GDO J=1,M V(K,J)=V(K,J)+Z(I,J)*G ENDDOENDDODO J=1,MV(K,J)=V(K,J)/D2ENDDOAPX(K)=APX(K)/D2BX(K)=D2/D1D1=D2ENDDOD1=MAPY(1)=0.0DO I=1,MAPY(1)=APY(1)+Y(I)ENDDOAPY(1)=APY(1)/D1DO J=1,PU(J,1)=0.0DO I=1,MU(J,1)=U(J,1)+V(J,I) ENDDO U(J,1)=U(J,1)/D1ENDDOIF(Q>1)THEND2=0.0APY(2)=0.0DO I=1,MG=Y(I)-APY(1)D2=D2+G*G APY(2)=APY(2)+(Y(I))*G*G ENDDO APY(2)=APY(2)/D2BY(2)=D2/D1DO J=1,PU(J,2)=0.0DO I=1,MG=Y(I)-APY(1)U(J,2)=U(J,2)+V(J,I)*GENDDOU(J,2)=U(J,2)/D2ENDDOD1=D2ENDIFDO K=3,QD2=0.0APY(K)=0.0DO J=1,PU(J,K)=0.0ENDDODO I=1,MG1=1.0G2=Y(I)-APY(1)DO J=3,KG=(Y(I)-APY(J-1))*G2-BY(J-1)*G1 G1=G2 G2=GENDDOD2=D2+G*GAPY(K)=APY(K)+Y(I)*G*G DO J=1,PU(J,K)=U(J,K)+V(J,I)*G ENDDOENDDODO J=1,PU(J,K)=U(J,K)/D2ENDDOAPY(K)=APY(K)/D2BY(K)=D2/D1D1=D2ENDDOV(1,1)=1.0V(2,1)=-APY(1)V(2,2)=1.0DO I=1,PDO J=1,QA(I,J)=0.0ENDDOENDDODO I=3,QV(I,I)=V(I-1,I-1)V(I,I-1)=-APY(I-1)*V(I-1,I-1)+V(I-1,I-2)IF(I>=4) THENDO K=I-2,2,-1V(I,K)=-APY(I-1)*V(I-1,K)+V(I-1,K-1)-BY(I-1)*V(I-2,K) ENDDO ENDIFV(I,1)=-APY(I-1)*V(I-1,1)-BY(I-1)*V(I-2,1)ENDDO DO I=1,PIF(I==1) THENT(1)=1.0T1(1)=1.0ELSEIF(I==2) THENT(1)=-APX(1)T(2)=1.0T2(1)=T(1)T2(2)=T(2)ELSET(I)=T2(I-1)T(I-1)=-APX(I-1)*T2(I-1)+T2(I-2) IF(I>=4) THENDO K=I-2,2,-1T(K)=-APX(I-1)*T2(K)+T2(K-1)-BX(I-1)*T1(K) ENDDOENDIFT(1)=-APX(I-1)*T2(1)-BX(I-1)*T1(1)T2(I)=T(I)DO K=I-1,1,-1T1(K)=T2(K)T2(K)=T(K)ENDDOENDIFDO J=1,QDO K=I,1,-1DO L=J,1,-1A(K,L)=A(K,L)+U(I,J)*T(K)*V(J,L) ENDDOENDDOENDDOENDDODT1=0.0DO I=1,NX1=X(I)DO J=1,MY1=Y(J)X2=1.0DD=0.0DO K=1,PG=A(K,Q)DO KK=Q-1,1,-1G=G*Y1+A(K,KK)ENDDOG=G*X2DD=DD+GX2=X2*X1ENDDODT=DD-Z(I,J)DT1=DT1+DT*DTENDDOENDDORETURNEND三、计算结果数表(x,y,f(x,y)): X Y UX TY F(X,Y) 0.00 0.500 1.345 0.243 0.17E+000.00 0.550 1.322 0.269 0.66E+000.00 0.600 1.299 0.295 0.35E+000.00 0.650 1.277 0.322 0.94E+000.00 0.700 1.255 0.350 0.30E-020.00 0.750 1.235 0.377 -0.87E-010.00 0.800 1.215 0.406 -0.58E+000.00 0.850 1.196 0.434 -0.72E+000.00 0.900 1.177 0.463 -0.54E+000.00 0.950 1.159 0.492 -0.86E+000.00 1.050 1.125 0.550 -0.74E+00 0.00 1.100 1.109 0.580 -0.06E+00 0.00 1.150 1.093 0.609 -0.00E+00 0.00 1.200 1.0790.639 -0.18E+00 0.00 1.250 1.064 0.669 -0.52E+00 0.00 1.3001.050 0.699 -0.19E+00 0.00 1.350 1.037 0.729 -0.48E+00 0.001.400 1.024 0.759 -0.68E+00 0.00 1.450 1.011 0.790 -0.52E+00 0.00 1.500 1.000 0.820 -0.29E+000.08 0.500 1.415 0.228 0.67E+00 0.08 0.550 1.391 0.253 0.08E+00 0.08 0.600 1.368 0.279 0.02E+00 0.08 0.650 1.346 0.306 0.47E+00 0.08 0.700 1.325 0.333 0.57E+00 0.08 0.750 1.304 0.360 0.48E-01 0.08 0.800 1.284 0.388 -0.73E-01 0.08 0.850 1.265 0.416 -0.16E+00 0.08 0.900 1.246 0.444 -0.29E+00 0.08 0.950 1.229 0.473 -0.36E+00 0.08 1.000 1.211 0.502 -0.08E+00 0.08 1.050 1.194 0.531 -0.29E+00 0.08 1.100 1.178 0.560 -0.78E+00 0.08 1.150 1.163 0.589 -0.93E+00 0.08 1.200 1.148 0.619 -0.44E+00 0.08 1.250 1.133 0.649 -0.92E+00 0.08 1.300 1.119 0.679 -0.71E+000.08 1.400 1.093 0.739 -0.37E+00 0.08 1.450 1.080 0.769-0.83E+00 0.08 1.500 1.068 0.799 -0.92E+000.16 0.500 1.483 0.214 0.31E+00 0.16 0.550 1.460 0.239 0.64E+00 0.16 0.600 1.437 0.264 0.91E+00 0.16 0.650 1.414 0.290 0.06E+00 0.16 0.700 1.393 0.316 0.70E+00 0.16 0.750 1.372 0.343 0.59E+00 0.16 0.800 1.352 0.370 0.12E+00 0.16 0.850 1.333 0.398 0.77E-02 0.16 0.900 1.315 0.426 -0.83E-01 0.16 0.950 1.297 0.454-0.58E+00 0.16 1.000 1.279 0.483 -0.20E+00 0.16 1.050 1.2620.512 -0.11E+00 0.16 1.100 1.246 0.541 -0.74E+00 0.16 1.1501.231 0.570 -0.09E+00 0.16 1.200 1.216 0.600 -0.59E+00 0.16 1.250 1.201 0.629 -0.66E+00 0.16 1.300 1.187 0.659 -0.71E+00 0.16 1.350 1.174 0.689 -0.32E+00 0.16 1.400 1.161 0.718-0.56E+00 0.16 1.450 1.148 0.748 -0.31E+00 0.16 1.500 1.136 0.778 -0.75E+000.24 0.500 1.551 0.201 0.66E+01 0.24 0.550 1.527 0.2250.03E+000.24 0.650 1.482 0.275 0.64E+00 0.24 0.700 1.460 0.3010.47E+00 0.24 0.750 1.439 0.327 0.34E+00 0.24 0.800 1.419 0.354 0.24E+00 0.24 0.850 1.400 0.381 0.69E+00 0.24 0.900 1.381 0.409 0.04E-01 0.24 0.950 1.363 0.437 -0.42E-01 0.24 1.000 1.346 0.465 -0.06E+00 0.24 1.050 1.329 0.494 -0.59E+00 0.24 1.100 1.313 0.523 -0.83E+00 0.24 1.150 1.297 0.552 -0.15E+00 0.24 1.200 1.282 0.581 -0.19E+00 0.24 1.250 1.267 0.610 -0.84E+00 0.24 1.300 1.253 0.640 -0.66E+00 0.24 1.350 1.240 0.669 -0.30E+00 0.24 1.400 1.227 0.699 -0.86E+00 0.24 1.450 1.214 0.729 -0.84E+00 0.24 1.500 1.202 0.759 -0.77E+000.32 0.500 1.617 0.188 0.28E+01 0.32 0.550 1.593 0.212 0.49E+01 0.32 0.600 1.570 0.236 0.68E+00 0.32 0.650 1.547 0.261 0.75E+00 0.32 0.700 1.526 0.286 0.60E+00 0.32 0.750 1.505 0.312 0.77E+00 0.32 0.800 1.485 0.339 0.05E+00 0.32 0.850 1.466 0.365 0.99E+00 0.32 0.900 1.447 0.393 0.27E+00 0.32 1.000 1.411 0.448 -0.01E-02 0.32 1.050 1.395 0.477-0.41E-01 0.32 1.100 1.378 0.505 -0.18E+00 0.32 1.150 1.3630.534 -0.25E+00 0.32 1.200 1.347 0.563 -0.29E+00 0.32 1.2501.333 0.592 -0.90E+00 0.32 1.300 1.319 0.621 -0.00E+00 0.32 1.350 1.305 0.650 -0.40E+00 0.32 1.400 1.292 0.680 -0.54E+00 0.32 1.450 1.279 0.710 -0.79E+00 0.32 1.500 1.267 0.739-0.91E+000.40 0.500 1.681 0.177 0.91E+01 0.40 0.550 1.658 0.1990.00E+01 0.40 0.600 1.634 0.223 0.83E+01 0.40 0.650 1.612 0.247 0.02E+01 0.40 0.700 1.591 0.272 0.94E+00 0.40 0.750 1.570 0.298 0.49E+00 0.40 0.800 1.550 0.324 0.94E+00 0.40 0.850 1.530 0.350 0.40E+00 0.40 0.900 1.512 0.377 0.33E+00 0.40 0.950 1.493 0.405 0.99E+00 0.40 1.000 1.476 0.432 0.68E+00 0.40 1.050 1.459 0.460 0.08E-01 0.40 1.100 1.443 0.488 -0.84E-01 0.40 1.150 1.427 0.517-0.98E+00 0.40 1.200 1.412 0.545 -0.27E+00 0.40 1.250 1.397 0.574 -0.06E+000.40 1.350 1.369 0.632 -0.66E+00 0.40 1.400 1.356 0.662-0.37E+00 0.40 1.450 1.343 0.691 -0.43E+00 0.40 1.500 1.331 0.721 -0.12E+000.48 0.500 1.745 0.166 0.69E+01 0.48 0.550 1.721 0.188 0.02E+01 0.48 0.600 1.698 0.211 0.74E+01 0.48 0.650 1.676 0.235 0.40E+01 0.48 0.700 1.654 0.259 0.23E+01 0.48 0.750 1.634 0.284 0.56E+00 0.48 0.800 1.613 0.310 0.28E+00 0.48 0.850 1.594 0.336 0.49E+00 0.48 0.900 1.575 0.363 0.31E+00 0.48 0.950 1.557 0.390 0.66E+00 0.48 1.000 1.539 0.417 0.30E+00 0.48 1.050 1.522 0.444 0.34E+00 0.48 1.100 1.506 0.472 0.07E-01 0.48 1.150 1.490 0.500 -0.62E-01 0.48 1.200 1.475 0.529 -0.45E+00 0.48 1.250 1.460 0.557 -0.86E+00 0.48 1.300 1.446 0.586 -0.39E+00 0.48 1.350 1.432 0.615 -0.22E+00 0.48 1.400 1.419 0.644 -0.67E+00 0.48 1.450 1.406 0.674-0.55E+00 0.48 1.500 1.394 0.703 -0.14E+000.56 0.500 1.808 0.156 0.48E+010.56 0.600 1.761 0.200 0.10E+01 0.56 0.650 1.739 0.2230.68E+01 0.56 0.700 1.717 0.247 0.94E+01 0.56 0.750 1.696 0.272 0.33E+01 0.56 0.800 1.676 0.297 0.11E+00 0.56 0.850 1.657 0.323 0.63E+00 0.56 0.900 1.638 0.349 0.97E+00 0.56 0.950 1.620 0.375 0.52E+00 0.56 1.000 1.602 0.402 0.56E+00 0.56 1.050 1.585 0.429 0.47E+00 0.56 1.100 1.568 0.457 0.20E+00 0.56 1.150 1.552 0.485 0.13E+00 0.56 1.200 1.537 0.513 0.09E-01 0.56 1.250 1.522 0.541 -0.47E-01 0.56 1.300 1.508 0.570 -0.99E+00 0.56 1.350 1.4940.599 -0.82E+00 0.56 1.400 1.481 0.627 -0.26E+00 0.56 1.4501.468 0.657 -0.71E+00 0.56 1.500 1.455 0.686 -0.98E+000.64 0.500 1.870 0.147 0.74E+01 0.64 0.550 1.846 0.1680.10E+01 0.64 0.600 1.823 0.190 0.54E+01 0.64 0.650 1.801 0.213 0.42E+01 0.64 0.700 1.779 0.236 0.56E+01 0.64 0.750 1.758 0.260 0.03E+01 0.64 0.800 1.738 0.285 0.42E+01 0.64 0.850 1.718 0.310 0.41E+010.64 0.950 1.681 0.362 0.36E+00 0.64 1.000 1.664 0.388 0.18E+00 0.64 1.050 1.646 0.415 0.28E+00 0.64 1.100 1.630 0.443 0.07E+00 0.64 1.150 1.614 0.470 0.66E+00 0.64 1.200 1.598 0.498 0.09E+00 0.64 1.250 1.584 0.526 0.50E-01 0.64 1.300 1.569 0.554 -0.88E-01 0.64 1.350 1.555 0.583 -0.76E+00 0.64 1.400 1.542 0.611 -0.66E+00 0.64 1.450 1.529 0.640 -0.33E+00 0.64 1.500 1.516 0.669 -0.56E+00 0.72 0.500 1.931 0.139 0.94E+01 0.72 0.550 1.907 0.159 0.84E+01 0.72 0.600 1.884 0.181 0.36E+01 0.72 0.650 1.862 0.203 0.40E+01 0.72 0.700 1.840 0.226 0.47E+01 0.72 0.750 1.819 0.249 0.56E+01 0.72 0.800 1.799 0.273 0.19E+01 0.72 0.850 1.779 0.298 0.37E+01 0.72 0.900 1.760 0.323 0.86E+01 0.72 0.950 1.742 0.349 0.76E+00 0.72 1.000 1.724 0.375 0.24E+00 0.72 1.050 1.707 0.402 0.55E+00 0.72 1.100 1.691 0.429 0.97E+00 0.72 1.150 1.675 0.456 0.27E+00 0.72 1.200 1.659 0.484 0.31E+000.72 1.300 1.630 0.539 0.49E+00 0.72 1.350 1.616 0.5680.72E-02 0.72 1.400 1.602 0.596 -0.69E-01 0.72 1.450 1.589 0.625 -0.67E+00 0.72 1.500 1.576 0.653 -0.20E+000.80 0.500 1.992 0.131 0.31E+01 0.80 0.550 1.968 0.1510.44E+01 0.80 0.600 1.945 0.172 0.41E+01 0.80 0.650 1.922 0.193 0.45E+01 0.80 0.700 1.900 0.216 0.00E+01 0.80 0.750 1.879 0.239 0.10E+01 0.80 0.800 1.859 0.263 0.16E+01 0.80 0.850 1.840 0.287 0.52E+01 0.80 0.900 1.821 0.312 0.02E+01 0.80 0.950 1.802 0.337 0.38E+01 0.80 1.000 1.784 0.363 0.89E+01 0.80 1.050 1.767 0.389 0.28E+00 0.80 1.100 1.751 0.416 0.09E+00 0.80 1.150 1.734 0.4430.23E+00 0.80 1.200 1.719 0.470 0.93E+00 0.80 1.250 1.704 0.498 0.15E+00 0.80 1.300 1.689 0.525 0.86E+00 0.80 1.350 1.675 0.553 0.64E+00 0.80 1.400 1.662 0.582 0.74E-01 0.80 1.450 1.649 0.610 -0.37E-01 0.80 1.500 1.636 0.638 -0.81E+00K和σ分别为：0 0.93E+031 0.61E+012 0.92E-023 0.53E-034 0.16E-055 0.77E-07系数矩阵Crs（按行）为：0.00E+01 -0.83E+01 0.56E+00 0.97E+00 -0.03E+00 0.70E-010.91E+01 -0.99E+00 -0.96E+01 0.17E+01 -0.66E+00 0.10E-01 0.77E+00 0.42E+01 -0.10E+00 -0.81E+00 0.81E+00 -0.62E-01-0.25E+00 -0.21E+00 0.97E+00 -0.18E+00 0.49E+00 -0.63E-010.34E+00 -0.56E+00 0.69E-01 0.51E+00 -0.77E-01 0.27E-01-0.94E-01 0.94E+00 -0.58E+00 0.69E-01 -0.50E-01 0.53E-02 数表（x,y,f(x,y),p(x,y)):X Y F(X,Y) P(X,Y)0.100 0.700 0.58E+00 0.05E+000.100 1.100 -0.66E+00 -0.26E+00 0.100 1.300 -0.68E+00-0.31E+00 0.100 1.500 -0.52E+00 -0.49E+000.200 0.700 0.54E+00 0.19E+00 0.200 0.900 -0.63E-01 -0.65E-01 0.200 1.100 -0.90E+00 -0.90E+00 0.200 1.300 -0.84E+00 -0.90E+00 0.200 1.500 -0.03E+00 -0.04E+000.300 0.700 0.82E+00 0.09E+00 0.300 0.900 0.48E+00 0.11E+00 0.300 1.100 -0.63E+00 -0.88E+00 0.300 1.300 -0.72E+00 -0.96E+00 0.300 1.500 -0.34E+00 -0.84E+000.400 0.700 0.79E+00 0.89E+00 0.400 0.900 0.56E+00 0.63E+00 0.400 1.100 -0.83E-01 -0.04E-01 0.400 1.300 -0.72E+00 -0.71E+00 0.400 1.500 -0.85E+00 -0.07E+000.500 0.700 0.56E+01 0.92E+01 0.500 0.900 0.51E+00 0.23E+00 0.500 1.100 0.59E+00 0.27E+00 0.500 1.300 -0.53E+00 -0.11E+00 0.500 1.500 -0.67E+00 -0.33E+000.600 0.900 0.14E+00 0.75E+00 0.600 1.100 0.19E+00 0.32E+00 0.600 1.300 -0.70E-01 -0.82E-01 0.600 1.500 -0.08E+00 -0.75E+00 0.700 0.700 0.89E+01 0.29E+01 0.700 0.900 0.91E+01 0.11E+010.700 1.100 0.60E+00 0.97E+00 0.700 1.300 0.22E-01 0.06E-01 0.7001.500 -0.53E+00 -0.80E+00 0.800 0.700 0.09E+01 0.06E+01 0.800 0.900 0.32E+01 0.50E+01 0.800 1.100 0.03E+00 0.79E+00 0.800 1.300 0.25E+00 0.50E+00 0.800 1.500 -0.14E+00 -0.28E+00。
北航数值分析大作业一
北京航空航天大学数值分析大作业一学院名称自动化专业方向控制工程学号ZY*******学生姓名许阳教师孙玉泉日期2021 年11月26 日设有501501⨯的实对称矩阵A ，⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=5011A a b c b c c b c b a其中，064.0,16.0),501,,2,1(64.0)2.0sin()024.064.1(1.0-==⋅⋅⋅=--=c b i e i i a ii 。
矩阵A 的特征值为)501,,2,1(⋅⋅⋅=i i λ，并且有||min ||,501150121i i s λλλλλ≤≤=≤⋅⋅⋅≤≤1λ，501λ和s λ的值。
A 的与数4015011λλλμ-+=kk 最接近的特征值)39,,2,1(⋅⋅⋅=k k i λ。
A 的(谱范数)条件数2)A (cond 和行列式detA 。
一 方案设计1 求1λ，501λ和s λ的值。
s λ为按模最小特征值，||min ||5011i i s λλ≤≤=。
可使用反幂法求得。
1λ，501λ分别为最大特征值及最小特征值。
可使用幂法求出按模最大特征值，如结果为正，即为501λ，结果为负，那么为1λ。
使用位移的方式求得另一特征值即可。
2 求A 的与数4015011λλλμ-+=kk 最接近的特征值)39,...,2,1(=k k i λ。
题目可看成求以k μ为偏移量后，按模最小的特征值。
即以k μ为偏移量做位移，使用反幂法求出按模最小特征值后，加上k μ，即为所求。
3 求A 的(谱范数)条件数2)(A cond 和行列式detA 。
矩阵A 为非奇异对称矩阵，可知，||)(min max2λλ=A cond(1-1)其中m ax λ为按模最大特征值，min λ为按模最小特征值。
detA 可由LU 分解得到。
因LU 均为三角阵，那么其主对角线乘积即为A 的行列式。
二 算法实现1 幂法使用如下迭代格式：⎪⎪⎩⎪⎪⎨⎧⋅===⋅⋅⋅=------||max |)|sgn(max ||max /),,(111111)0()0(10k k k k k k k k Tn u u Ay u u u y u u u β任取非零向量 (2-1)终止迭代的控制理论使用εβββ≤--||/||1k k k ， 实际使用εβββ≤--||/||||||1k k k(2-2)由于不保存A 矩阵中的零元素，只保存主对角元素a[501]及b,c 值。
数值分析大作业四
《数值分析》大作业四一、算法设计方案：复化梯形积分法，选取步长为1/500=0.002，迭代误差控制在E ≤1.0e-10①复化梯形积分法：11()[()()2()]2n bak hf x dx f a f b f a kh -=⎰≈+++∑，截断误差为：322()''()''(),[,]1212T b a b a R f h f a b n ηηη--=-=-∈其中。
复化Simpson 积分法，选取步长为1/50=0.02，迭代误差控制在E ≤1.0e-10②Simpson 积分法：121211()[()()4()2()]3m m bi i a i i hf x dx f a f b f x f x --==≈+++∑∑⎰， 截断误差为：4(4)(),[,]180s b a R h f a b ηη-=-∈。
③Guass积分法选用Gauss-Legendre 求积公式：111()()ni i i f x dx A f x -=≈∑⎰截断误差为：R= ()()n 2n 422n!2×(2[2!]2n 1f n n ⨯（2）η（））＋ η∈(1,1)。
选择9个节点：-0.9681602395,-0.8360311073,-0.6133714327,-0.3242534234,0,0.3242534234,0.6133714327,0.8360311073,0.9681602395， 对应的求积系数依次为：0.0812743884,0.1806481607,0.2606106964,0.3123470770,0.3302393550,0.3123470770,0.2606106964,0.1806481607,0.0812743884。
二、程序源代码：#include<stdio.h>#include<math.h>#include<stdlib.h>#define E 1.0e-10/****定义函数g和K*****/double g(double a){double b;b=exp(4*a)+(exp(a+4)-exp(-a-4))/(a+4);return b;}double K(double a,double b){double c;c=exp(a*b);return c;}/******复化梯形法******/void Tixing( ){double u[1001],x[1001],h,c[1001],e;int i,j,k;FILE *fp;fp=fopen("f:/result0. xls ","w");h=1.0/1500;for(i=0;i<3001;i++){x[i]=i*h-1;u[i]=g(x[i]);}for(k=0;k<100;k++){e=0;for(i=0;i<1001;i++){for(j=1,c[i]=0;j<N-1;j++)c[i]+=K(x[i],x[j])*u[j];u[i]=g(x[i])-h*c[i]-h/2*(K(x[i],x[0])*u[0]+K(x[i],x[N-1])*u[N-1]);e+=h*(exp(4*x[i])-u[i])*(exp(4*x[i])-u[i]);}if(e<=E) break;}for(i=0;i<1001;i++)fprintf(fp,"%.12lf,%.12lf\n",x[i],u[i]);fclose(fp);}/******复化Simpson法******/void simpson( ){double u[101],x[101],h,c[101],d[101],e;int i,j,k;FILE *fp;fp=fopen("f:/result1.xls","w");h=1.0/50;for(i=0;i<101;i++){x[i]=i*h-1;u[i]=g(x[i]);}for(k=0;k<50;k++){e=0;for(i=0;i<101;i++){for(j=1,c[i]=0,d[i]=0;j<51;j++){c[i]+=K(x[i],x[2*j-1])*u[2*j-1];if(j<50)d[i]+=K(x[i],x[2*j])*u[2*j];}u[i]=g(x[i])-4*h/3*c[i]-2*h/3*d[i]-h/3*(K(x[i],x[0])*u[0]+K(x[i],x[M-1])*u[M-1]);e+=h*(exp(4*x[i])-u[i])*(exp(4*x[i])-u[i]);}if(e<=E) break;}for(i=0;i<101;i++)fprintf(fp,"%.12lf,%.12lf\n",x[i],u[i]);fclose(fp);}/******Gauss积分法******/void gauss( ){double x[9]={-0.9681602395,-0.8360311073,-0.6133714327,-0.3242534234,0,\0.3242534234,0.6133714327,0.8360311073,0.9681602395},A[9]={0.0812743884,0.1806481607,0.2606106964,0.3123470770,0.3302393550,\0.3123470770,0.2606106964,0.1806481607,0.0812743884},u[9],c[9],e;int i,j,k;FILE *fp;fp=fopen("f:/result2. xls ","w");for(i=0;i<9;i++)u[i]=g(x[i]);for(k=0;k<50;k++){e=0;for(i=0;i<9;i++){for(j=0,c[i]=0;j<9;j++)c[i]+=A[j]*K(x[i],x[j])*u[j];u[i]=g(x[i])-c[i];e+=A[i]*(exp(4*x[i])-u[i])*(exp(4*x[i])-u[i]);}if(e<=E) break;}for(i=0;i<9;i++)fprintf(fp,"%.12lf,%.12lf\n",x[i],u[i]);fclose(fp);}/******主函数******/main(){Tixing ( );Simpson( );Gauss( );return 0;}三、运算结果复化梯形数据-10.018323-0.920.02523-0.9980.018471-0.9180.025433-0.9960.018619-0.9160.025637-0.9940.018768-0.9140.025843-0.9920.018919-0.9120.026051-0.990.019071-0.910.02626-0.9880.019224-0.9080.026471-0.9860.019378-0.9060.026683-0.9840.019534-0.9040.026897-0.9820.019691-0.9020.027113-0.980.019849-0.90.027331-0.9780.020008-0.8980.02755-0.9760.020169-0.8960.027772-0.9740.020331-0.8940.027995-0.9720.020494-0.8920.028219-0.970.020658-0.890.028446-0.9680.020824-0.8880.028674-0.9660.020992-0.8860.028905-0.9640.02116-0.8840.029137-0.9620.02133-0.8820.029371-0.960.021501-0.880.029607-0.9580.021674-0.8780.029844-0.9560.021848-0.8760.030084-0.9540.022023-0.8740.030326-0.9520.0222-0.8720.030569-0.950.022378-0.870.030815-0.9480.022558-0.8680.031062-0.9460.022739-0.8660.031311-0.9440.022922-0.8640.031563-0.9420.023106-0.8620.031816-0.940.023291-0.860.032072-0.9380.023478-0.8580.032329-0.9360.023667-0.8560.032589-0.9340.023857-0.8540.032851-0.9320.024048-0.8520.033114-0.930.024241-0.850.03338-0.9280.024436-0.8480.033648-0.9260.024632-0.8460.033918-0.9240.02483-0.8440.034191-0.9220.025029-0.8420.034465-0.840.034742-0.760.047841-0.8380.035021-0.7580.048225-0.8360.035302-0.7560.048613 -0.8340.035586-0.7540.049003 -0.8320.035872-0.7520.049396 -0.830.03616-0.750.049793 -0.8280.03645-0.7480.050193 -0.8260.036743-0.7460.050596 -0.8240.037038-0.7440.051002 -0.8220.037335-0.7420.051412 -0.820.037635-0.740.051825 -0.8180.037937-0.7380.052241 -0.8160.038242-0.7360.052661 -0.8140.038549-0.7340.053084 -0.8120.038858-0.7320.05351 -0.810.039171-0.730.05394 -0.8080.039485-0.7280.054373 -0.8060.039802-0.7260.054809 -0.8040.040122-0.7240.05525 -0.8020.040444-0.7220.055693 -0.80.040769-0.720.056141 -0.7980.041096-0.7180.056591 -0.7960.041426-0.7160.057046 -0.7940.041759-0.7140.057504 -0.7920.042094-0.7120.057966 -0.790.042432-0.710.058431 -0.7880.042773-0.7080.058901 -0.7860.043116-0.7060.059374 -0.7840.043463-0.7040.05985 -0.7820.043812-0.7020.060331 -0.780.044164-0.70.060816 -0.7780.044518-0.6980.061304 -0.7760.044876-0.6960.061796 -0.7740.045236-0.6940.062293 -0.7720.045599-0.6920.062793 -0.770.045966-0.690.063297 -0.7680.046335-0.6880.063805 -0.7660.046707-0.6860.064318 -0.7640.047082-0.6840.064834 -0.7620.04746-0.6820.065355-0.680.06588-0.60.090722 -0.6780.066409-0.5980.091451-0.6760.066942-0.5960.092185 -0.6740.06748-0.5940.092926 -0.6720.068022-0.5920.093672 -0.670.068568-0.590.094424 -0.6680.069119-0.5880.095183 -0.6660.069674-0.5860.095947 -0.6640.070234-0.5840.096718 -0.6620.070798-0.5820.097494 -0.660.071366-0.580.098277 -0.6580.071939-0.5780.099067 -0.6560.072517-0.5760.099862 -0.6540.0731-0.5740.100664 -0.6520.073687-0.5720.101473 -0.650.074278-0.570.102288 -0.6480.074875-0.5680.103109 -0.6460.075476-0.5660.103937 -0.6440.076082-0.5640.104772 -0.6420.076694-0.5620.105614 -0.640.077309-0.560.106462 -0.6380.07793-0.5580.107317 -0.6360.078556-0.5560.108179 -0.6340.079187-0.5540.109048 -0.6320.079823-0.5520.109924 -0.630.080464-0.550.110806 -0.6280.08111-0.5480.111696 -0.6260.081762-0.5460.112593 -0.6240.082418-0.5440.113498 -0.6220.08308-0.5420.114409 -0.620.083748-0.540.115328 -0.6180.08442-0.5380.116254 -0.6160.085098-0.5360.117188 -0.6140.085782-0.5340.118129 -0.6120.086471-0.5320.119078 -0.610.087165-0.530.120035 -0.6080.087865-0.5280.120999 -0.6060.088571-0.5260.12197 -0.6040.089282-0.5240.12295 -0.6020.089999-0.5220.123938-0.550.110806-0.470.152592 -0.5480.111696-0.4680.153817-0.5460.112593-0.4660.155053-0.5440.113498-0.4640.156298-0.5420.114409-0.4620.157553-0.540.115328-0.460.158819-0.5380.116254-0.4580.160095-0.5360.117188-0.4560.16138-0.5340.118129-0.4540.162677-0.5320.119078-0.4520.163983-0.530.120035-0.450.1653-0.5280.120999-0.4480.166628-0.5260.12197-0.4460.167966-0.5240.12295-0.4440.169315-0.5220.123938-0.4420.170675-0.520.124933-0.440.172046-0.5180.125936-0.4380.173428-0.5160.126948-0.4360.174821-0.5140.127967-0.4340.176225-0.5120.128995-0.4320.17764-0.510.130031-0.430.179067-0.5080.131076-0.4280.180505-0.5060.132128-0.4260.181955-0.5040.13319-0.4240.183416-0.5020.134259-0.4220.18489-0.50.135338-0.420.186375-0.4980.136425-0.4180.187871-0.4960.13752-0.4160.18938-0.4940.138625-0.4140.190901-0.4920.139738-0.4120.192435-0.490.140861-0.410.19398-0.4880.141992-0.4080.195538-0.4860.143132-0.4060.197109-0.4840.144282-0.4040.198692-0.4820.145441-0.4020.200288-0.480.146609-0.40.201897-0.4780.147786-0.3980.203518-0.4760.148973-0.3960.205153-0.4740.15017-0.3940.206801-0.4720.151376-0.3920.208462-0.390.210136-0.310.289382 -0.3880.211824-0.3080.291706-0.3860.213525-0.3060.294049 -0.3840.21524-0.3040.296411 -0.3820.216969-0.3020.298792 -0.380.218711-0.30.301192 -0.3780.220468-0.2980.303611 -0.3760.222239-0.2960.306049 -0.3740.224024-0.2940.308508 -0.3720.225823-0.2920.310985 -0.370.227637-0.290.313483 -0.3680.229465-0.2880.316001 -0.3660.231308-0.2860.318539 -0.3640.233166-0.2840.321098 -0.3620.235039-0.2820.323677 -0.360.236927-0.280.326277 -0.3580.23883-0.2780.328897 -0.3560.240748-0.2760.331539 -0.3540.242682-0.2740.334202 -0.3520.244631-0.2720.336886 -0.350.246596-0.270.339592 -0.3480.248576-0.2680.34232 -0.3460.250573-0.2660.345069 -0.3440.252586-0.2640.347841 -0.3420.254614-0.2620.350635 -0.340.256659-0.260.353451 -0.3380.258721-0.2580.35629 -0.3360.260799-0.2560.359151 -0.3340.262894-0.2540.362036 -0.3320.265005-0.2520.364944 -0.330.267134-0.250.367875 -0.3280.269279-0.2480.37083 -0.3260.271442-0.2460.373809 -0.3240.273622-0.2440.376811 -0.3220.27582-0.2420.379838 -0.320.278035-0.240.382888 -0.3180.280268-0.2380.385964 -0.3160.28252-0.2360.389064 -0.3140.284789-0.2340.392189 -0.3120.287076-0.2320.395339-0.230.398514-0.150.548804-0.2280.401715-0.1480.553212-0.2260.404942-0.1460.557655 -0.2240.408194-0.1440.562134 -0.2220.411473-0.1420.56665 -0.220.414778-0.140.571201 -0.2180.418109-0.1380.575789 -0.2160.421467-0.1360.580414 -0.2140.424853-0.1340.585076 -0.2120.428265-0.1320.589775 -0.210.431705-0.130.594512 -0.2080.435172-0.1280.599287 -0.2060.438668-0.1260.604101 -0.2040.442191-0.1240.608953 -0.2020.445743-0.1220.613844 -0.20.449323-0.120.618774 -0.1980.452932-0.1180.623744 -0.1960.45657-0.1160.628754 -0.1940.460237-0.1140.633805 -0.1920.463934-0.1120.638895 -0.190.46766-0.110.644027 -0.1880.471416-0.1080.6492 -0.1860.475203-0.1060.654414 -0.1840.47902-0.1040.659671 -0.1820.482867-0.1020.664969 -0.180.486746-0.10.67031 -0.1780.490655-0.0980.675694 -0.1760.494596-0.0960.681121 -0.1740.498569-0.0940.686592 -0.1720.502573-0.0920.692107 -0.170.50661-0.090.697666 -0.1680.510679-0.0880.70327 -0.1660.514781-0.0860.708919 -0.1640.518916-0.0840.714613 -0.1620.523084-0.0820.720352 -0.160.527285-0.080.726138 -0.1580.53152-0.0780.731971 -0.1560.535789-0.0760.73785 -0.1540.540093-0.0740.743776 -0.1520.544431-0.0720.749751-0.070.7557730.01 1.040796 -0.0680.7618430.012 1.049156-0.0660.7679620.014 1.057583 -0.0640.7741310.016 1.066077 -0.0620.7803480.018 1.07464 -0.060.7866160.02 1.083272 -0.0580.7929340.022 1.091973 -0.0560.7993030.024 1.100743 -0.0540.8057230.026 1.109585 -0.0520.8121950.028 1.118497 -0.050.8187190.03 1.127481 -0.0480.8252950.032 1.136537 -0.0460.8319240.034 1.145666 -0.0440.8386060.036 1.154868 -0.0420.8453410.038 1.164144 -0.040.8521310.04 1.173494 -0.0380.8589760.042 1.18292 -0.0360.8658750.044 1.192421 -0.0340.872830.046 1.201999 -0.0320.879840.048 1.211654 -0.030.8869070.05 1.221386 -0.0280.8940310.052 1.231196 -0.0260.9012120.054 1.241085 -0.0240.9084510.056 1.251054 -0.0220.9157480.058 1.261102 -0.020.9231030.06 1.271232 -0.0180.9305170.062 1.281442 -0.0160.9379910.064 1.291735 -0.0140.9455250.066 1.30211 -0.0120.953120.068 1.312569 -0.010.9607750.07 1.323112 -0.0080.9684930.072 1.333739 -0.0060.9762720.074 1.344452 -0.0040.9841130.076 1.355251 -0.0020.9920180.078 1.366136 00.9999860.08 1.377109 0.002 1.0080180.082 1.38817 0.004 1.0161140.084 1.39932 0.006 1.0242760.086 1.41056 0.008 1.0325030.088 1.4218890.09 1.433310.17 1.973853 0.092 1.4448230.172 1.9897080.094 1.4564280.174 2.005689 0.096 1.4681260.176 2.021799 0.098 1.4799180.178 2.038039 0.1 1.4918050.18 2.054408 0.102 1.5037870.182 2.07091 0.104 1.5158660.184 2.087543 0.106 1.5280410.186 2.104311 0.108 1.5403150.188 2.121213 0.11 1.5526870.19 2.138251 0.112 1.5651580.192 2.155425 0.114 1.577730.194 2.172738 0.116 1.5904020.196 2.19019 0.118 1.6031760.198 2.207781 0.12 1.6160530.2 2.225515 0.122 1.6290340.202 2.24339 0.124 1.6421180.204 2.261409 0.126 1.6553080.206 2.279573 0.128 1.6686040.208 2.297883 0.13 1.6820060.21 2.31634 0.132 1.6955160.212 2.334945 0.134 1.7091350.214 2.3537 0.136 1.7228630.216 2.372605 0.138 1.7367010.218 2.391662 0.14 1.750650.22 2.410872 0.142 1.7647120.222 2.430236 0.144 1.7788860.224 2.449756 0.146 1.7931740.226 2.469433 0.148 1.8075770.228 2.489268 0.15 1.8220960.23 2.509262 0.152 1.8367310.232 2.529417 0.154 1.8514840.234 2.549733 0.156 1.8663550.236 2.570213 0.158 1.8813460.238 2.590857 0.16 1.8964570.24 2.611667 0.162 1.911690.242 2.632645 0.164 1.9270450.244 2.65379 0.166 1.9425230.246 2.675106 0.168 1.9581260.248 2.6965930.25 2.7182520.33 3.743385 0.252 2.7400850.332 3.7734530.254 2.7620940.334 3.803761 0.256 2.7842790.336 3.834314 0.258 2.8066430.338 3.865111 0.26 2.8291860.34 3.896156 0.262 2.8519110.342 3.927451 0.264 2.8748180.344 3.958996 0.266 2.8979090.346 3.990796 0.268 2.9211850.348 4.02285 0.27 2.9446480.35 4.055162 0.272 2.96830.352 4.087734 0.274 2.9921420.354 4.120567 0.276 3.0161750.356 4.153664 0.278 3.0404010.358 4.187026 0.28 3.0648220.36 4.220657 0.282 3.0894390.362 4.254558 0.284 3.1142540.364 4.288731 0.286 3.1392680.366 4.323179 0.288 3.1644830.368 4.357903 0.29 3.18990.37 4.392906 0.292 3.2155220.372 4.42819 0.294 3.2413490.374 4.463758 0.296 3.2673840.376 4.499612 0.298 3.2936280.378 4.535753 0.3 3.3200830.38 4.572185 0.302 3.346750.382 4.608909 0.304 3.3736320.384 4.645928 0.306 3.4007290.386 4.683245 0.308 3.4280440.388 4.720861 0.31 3.4555790.39 4.75878 0.312 3.4833350.392 4.797003 0.314 3.5113130.394 4.835533 0.316 3.5395160.396 4.874373 0.318 3.5679460.398 4.913524 0.32 3.5966040.4 4.95299 0.322 3.6254930.402 4.992773 0.324 3.6546130.404 5.032876 0.326 3.6839670.406 5.0733 0.328 3.7135570.408 5.114050.41 5.1551260.497.099276 0.412 5.1965330.4927.1562980.414 5.2382720.4947.213778 0.416 5.2803460.4967.27172 0.418 5.3227590.4987.330127 0.42 5.3655120.57.389004 0.422 5.4086080.5027.448353 0.424 5.4520510.5047.508179 0.426 5.4958420.5067.568486 0.428 5.5399850.5087.629277 0.43 5.5844830.517.690556 0.432 5.6293380.5127.752327 0.434 5.6745540.5147.814595 0.436 5.7201330.5167.877362 0.438 5.7660770.5187.940634 0.44 5.8123910.528.004414 0.442 5.8590770.5228.068707 0.444 5.9061380.5248.133516 0.446 5.9535770.5268.198845 0.448 6.0013960.5288.264699 0.45 6.04960.538.331082 0.452 6.0981910.5328.397998 0.454 6.1471730.5348.465452 0.456 6.1965480.5368.533447 0.458 6.2463190.5388.601989 0.46 6.296490.548.671081 0.462 6.3470640.5428.740728 0.464 6.3980450.5448.810935 0.466 6.4494340.5468.881705 0.468 6.5012370.5488.953044 0.47 6.5534560.559.024956 0.472 6.6060940.5529.097445 0.474 6.6591550.5549.170517 0.476 6.7126420.5569.244175 0.478 6.7665580.5589.318426 0.48 6.8209080.569.393272 0.482 6.8756950.5629.46872 0.484 6.9309210.5649.544774 0.486 6.9865910.5669.621439 0.4887.0427080.5689.6987190.579.776620.6513.46367 0.5729.8551470.65213.571810.5749.9343050.65413.68082 0.57610.01410.65613.79071 0.57810.094530.65813.90147 0.5810.175610.6614.01313 0.58210.257340.66214.12569 0.58410.339730.66414.23915 0.58610.422780.66614.35352 0.58810.50650.66814.46881 0.5910.590890.6714.58502 0.59210.675960.67214.70217 0.59410.761710.67414.82026 0.59610.848150.67614.9393 0.59810.935280.67815.05929 0.611.023110.6815.18025 0.60211.111650.68215.30218 0.60411.20090.68415.42509 0.60611.290870.68615.54898 0.60811.381560.68815.67387 0.6111.472980.6915.79977 0.61211.565130.69215.92667 0.61411.658020.69416.0546 0.61611.751660.69616.18355 0.61811.846050.69816.31354 0.6211.94120.716.44457 0.62212.037110.70216.57665 0.62412.133790.70416.7098 0.62612.231250.70616.84401 0.62812.32950.70816.97931 0.6312.428530.7117.11569 0.63212.528360.71217.25316 0.63412.628990.71417.39174 0.63612.730420.71617.53143 0.63812.832680.71817.67225 0.6412.935750.7217.81419 0.64213.039650.72217.95728 0.64413.144390.72418.10151 0.64613.249960.72618.24691 0.64813.356390.72818.393470.7318.541210.8125.53363 0.73218.690130.81225.738720.73418.840250.81425.94545 0.73618.991580.81626.15385 0.73819.144120.81826.36392 0.7419.297890.8226.57568 0.74219.452890.82226.78914 0.74419.609140.82427.00431 0.74619.766640.82627.22121 0.74819.925410.82827.43985 0.7520.085450.8327.66025 0.75220.246780.83227.88242 0.75420.409410.83428.10638 0.75620.573340.83628.33213 0.75820.738580.83828.5597 0.7620.905160.8428.78909 0.76221.073070.84229.02033 0.76421.242330.84429.25342 0.76621.412950.84629.48839 0.76821.584940.84829.72524 0.7721.758310.8529.964 0.77221.933080.85230.20467 0.77422.109250.85430.44728 0.77622.286830.85630.69184 0.77822.465840.85830.93836 0.7822.646290.8631.18686 0.78222.828190.86231.43735 0.78423.011550.86431.68986 0.78623.196380.86631.9444 0.78823.382690.86832.20098 0.7923.570510.8732.45962 0.79223.759830.87232.72034 0.79423.950670.87432.98315 0.79624.143040.87633.24807 0.79824.336960.87833.51513 0.824.532440.8833.78432 0.80224.729490.88234.05568 0.80424.928110.88434.32922 0.80625.128340.88634.60496 0.80825.330170.88834.882910.8935.163090.94643.99154 0.89235.445520.94844.344880.89435.730220.9544.701070.89636.017210.95245.060110.89836.306510.95445.422040.936.598120.95645.786870.90236.892080.95846.154630.90437.188410.9646.525350.90637.487110.96246.899050.90837.788210.96447.275750.9138.091730.96647.655470.91238.397680.96848.038240.91438.70610.9748.424090.91639.016990.97248.813040.91839.330380.97449.205110.9239.646280.97649.600330.92239.964720.97849.998720.92440.285720.9850.400320.92640.60930.98250.805140.92840.935480.98451.213210.9341.264280.98651.624560.93241.595720.98852.039210.93441.929820.9952.45720.93642.26660.99252.878540.93842.606090.99453.303270.9442.948310.99653.73140.94243.293270.99854.162980.94443.64101154.59802复化Simpson数据：-1 0.018319929 -0.34 0.256658088 0.32 3.596641805 -0.98 0.0198445 -0.32 0.278035042 0.34 3.896195298-0.96 0.021494322 -0.3 0.301192133 0.36 4.220697765-0.94 0.023283225 -0.28 0.326278124 0.38 4.572227037-0.92 0.025220379 -0.26 0.353453177 0.4 4.95303418-0.9 0.027320224 -0.24 0.382891765 0.42 5.365557596-0.88 0.029594431 -0.22 0.41478194 0.44 5.812438891-0.86 0.032059069 -0.16 0.527292277 0.54 8.671138204-0.84 0.034728638 -0.14 0.571209036 0.56 9.39333156-0.82 0.037621263 -0.12 0.61878367 0.58 10.17567433-0.8 0.040754615 -0.1 0.670320427 0.6 11.02317608-0.78 0.044149394 -0.08 0.726149698 0.62 11.94126383-0.76 0.047826844 -0.06 0.78662861 0.64 12.93581634-0.74 0.051810827 -0.04 0.85214479 0.66 14.01320231-0.72 0.056126648 -0.02 0.92311742 0.68 15.1803205-0.7 0.060802006 0 1.0000013 0.7 16.44464467 -0.68 0.065866854 0.02 1.083288424 0.72 17.81427057 -0.66 0.071353499 0.04 1.173512427 0.74 19.29796874 -0.64 0.077297255 0.06 1.271250748 0.76 20.90523965 -0.62 0.083735917 0.08 1.377129533 0.78 22.64637562 -0.6 0.090711017 0.1 1.491826493 0.8 24.53252554 -0.58 0.098266855 0.12 1.616076341 0.82 26.57576756 -0.56 0.106452202 0.14 1.750674449 0.84 28.78918506 -0.54 0.11531904 0.16 1.896482943 0.86 31.18695183 -0.52 0.12492459 0.18 2.054435268 0.88 33.78442141 -0.5 0.135329888 0.2 2.225543071 0.9 36.59822683 -0.48 0.14660204 0.22 2.410901825 0.92 39.64638571 -0.46 0.158812728 0.24 2.611698647 0.94 42.94841704 -0.44 0.17204064 0.26 2.829219145 0.96 46.52546475 -0.42 0.18636997 0.28 3.064856356 0.98 50.40043451 -0.4 0.201892977 0.3 3.320119013 1 54.59813904 -0.38 0.218708553 0.46 6.296539601-0.36 0.236924875 0.48 6.820959636-0.2 0.449328351 0.5 7.389057081-0.18 0.486751777 0.52 8.0044696750102030405060四、讨论①在满足相同精度要求的情况下复化梯形积分法比复化Simpson 积分法计算所需节点数多，计算量大。
数值分析大作业课题一
迭代格式的比较1、自定义函数：function y=fun1(x)y=(3*x+1)/x^2;主程序：x0=input('Please input the initial value£ºx0=');k=input('Please input the number of iterations£ºk='); x(1)=x0;i=1;while i<=kx(i+1)=fun1(x(i));i=i+1;xk=x(i);endfprintf(' x=%f\n',x);运行结果：Please input the initial value：x0=1Please input the number of iterations：k=5x=1.000000x=4.000000x=0.812500x=5.207101x=0.613018x=7.554877结论：发散2、自定义函数：function y=fun2(x)y=(x^3-1)/3;主程序：x0=input('Please input the initial value£ºx0=');k=input('Please input the number of iterations£ºk='); x(1)=x0;i=1;while i<=kx(i+1)=fun2(x(i));i=i+1;xk=x(i);endfprintf(' x=%f\n',x);运行结果：Please input the initial value：x0=1Please input the number of iterations：k=8x=1.000000x=0.000000x=-0.333333x=-0.345679x=-0.347102x=-0.347273x=-0.347294x=-0.347296x=-0.347296结论：收敛3、自定义函数：function y=fun3(x)y=(3*x+1)^(1/3);主程序：x0=input('Please input the initial value£ºx0=');k=input('Please input the number of iterations£ºk='); x(1)=x0;i=1;while i<=kx(i+1)=fun3(x(i));i=i+1;xk=x(i);endfprintf(' x=%f\n',x);运行结果：Please input the initial value：x0=1Please input the number of iterations：k=11x=1.000000x=1.587401x=1.792790x=1.854542x=1.872325x=1.877384x=1.878819x=1.879225x=1.879340x=1.879372x=1.879382x=1.879384结论：收敛4、自定义函数：function y=fun4(x)y=1/(x^2-3);主程序：x0=input('Please input the initial value£ºx0=');k=input('Please input the number of iterations£ºk='); x(1)=x0;i=1;while i<=kx(i+1)=fun4(x(i));i=i+1;xk=x(i);endfprintf(' x=%f\n',x);运行结果：Please input the initial value：x0=1Please input the number of iterations：k=9x=1.000000x=-0.500000x=-0.363636x=-0.348703x=-0.347414x=-0.347306x=-0.347297x=-0.347296x=-0.347296x=-0.347296结论：收敛5、自定义函数：function y=fun5(x)y=sqrt(3+1/x);主程序：x0=input('Please input the initial value£ºx0=');k=input('Please input the number of iterations£ºk='); x(1)=x0;i=1;while i<=kx(i+1)=fun5(x(i));i=i+1;xk=x(i);endfprintf(' x=%f\n',x);运行结果：Please input the initial value：x0=1Please input the number of iterations：k=6x=1.000000x=2.000000x=1.870829x=1.880033x=1.879336x=1.879389x=1.879385结论：收敛6、自定义函数：function y=fun6(x)y=x-1/3*((x^3-3*x-1)/(x^2-1));主程序：x0=input('Please input the initial value£ºx0=');k=input('Please input the number of iterations£ºk='); x(1)=x0;i=1;while i<=kx(i+1)=fun6(x(i));i=i+1;xk=x(i);endfprintf(' x=%f\n',x);运行结果：Please input the initial value：x0=3Please input the number of iterations：k=6x=3.000000x=2.291667x=1.965507x=1.884402x=1.879404x=1.879385x=1.879385结论：收敛要求2，输出迭代次数：自定义函数：interation.m%x0为初始值，k为最大迭代次数，e为控制精度function interation(x0,k,e)x(1)=x0;i=1;while i<=kx(i+1)=fun6(x(i));if abs(x(i+1)-x(i))<efprintf(' the number of iterations: n=%f\n',i);breakendi=i+1;endx=x(i);fprintf(' a root of the original equation: x=%8.6f\n',x);取形式6为例，运行结果：>>interation(5,24,0.00001)the number of iterations: n=7.000000a root of the original equation: x=1.879385。
BUAA数值分析大作业三
北京航空航天大学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 相同。
北航数值分析大作业3(学硕)
《数值分析》作业三院系：机械学院学号：SY1307145姓名：龙安林2013年11 月24 日1. 算法设计1) 开始；2) 计算数组[][]0.08,0.050.5,0,1,2,,10;0,1,2,,20x i i y j j i j ==+=⋯=⋯（）； 3) 将点[][],0,1,2,,10;0,1,2,,20x i y j i j =⋯=⋯（）,（）带入非线性方程组： 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 +++-=⎧⎪+++-=⎪⎨+++-=⎪⎪+++-=⎩ 得出相应的点,t u （）； 4) 选择拉格朗日插值法，将,t u （）作为中间变量，在题目所给出的二维数表中进行二次代数插值，得到[][],)(z f x i y j =；5) 输出数表：[][][][]()()0,1,2,,10;0,1,2,,20,,,x i y j f x i y j i j =⋯=⋯； 6) 令k=0；7) 以()()(),,,0,1,r r r s x x y y r s ϕψ===…,k 为拟合基函数，将上述数表作为拟合条件，对于给定的k 值，得到矩阵B 、G 、U ；8) 令-1-1(),()T T T A B B B U C AG G G ==，用选主元的LU 分解法分别计算矩阵A 和C 的各列，最后得到系数矩阵C ；9) 以公式：()()()00,k ki j rs r i s j s r p x y C x y ϕψ===∑∑计算每个点的拟合值；10) 利用公式：()()()2102000,,i j i j i j f x y p x y σ===-∑∑计算拟合误差，当σ≤10-7时，循环结束，否则k=k+1，转（6）；11) 令[][]()**0.10.50.2 1,2,81,2,5x i i y j j i j ==+=⋯=⋯；，；，；12) 计算()()()******,,,,,i j i j i jf x y p x y delta x y ，输出数表，观察逼近效果； 13) 结束。
数值分析大作业
数值分析大作业数值分析大作业姓名：黄晨晨学号：S1*******学院：储运与建筑工程学院学院班级：储建研17-2实验3.1 Gauss消去法的数值稳定性实验实验目的：理解高斯消元过程中出现小主元即很小时引起方程组解数值不定性实验内容：求解方程组Ax=b，其中（1）A1=0.3×10?1559.14315.291?6.130?1211.29521211，b1=59.1746.7812;（2）A2=10?7013 2.099999999999625?15?10102，b2=85.90000000000151;实验要求：(1)计算矩阵的条件数，判断系数矩阵是良态的还是病态的(2)用Gauss列主元消去法求得L和U及解向量x1,x2∈R4(3)用不选主元的高斯消去法求得L和U及解向量x1,x2∈R4(4)观察小主元并分析对计算结果的影响（1）计算矩阵的条件数，判断系数矩阵是良态的还是病态的代码：format longeformat compactA1=[0.3*10^-15,59.14,3,1;5.291,-6.130,-1,2;11.2,9,5,2;1,2,1,1] b1=[59.17;46.78;1;2]n=4C1=cond(A1,1) %C1为A1矩阵1范数下的条件数C2=cond(A1,2) %C2为A1矩阵2范数下的条件数C3=cond(A1,inf) %C3为1矩阵谱范数下的条件数结果：C1 =1.362944708720448e+02C2 =6.842955771253409e+01C3 =8.431146*********e+01显然A1矩阵为病态矩阵将矩阵A2，b2输入上述代码中求得A2矩阵的条件数为：C1 =1.928316831682894e+01C2 =8.993938090170119e+00C3 =1.835643564356072e+01显然A2矩阵也为病态矩阵(2)用Gauss列主元消去法求得L和U及解向量x1,x2∈R4代码：for k=1:n-1a=max(abs(A1(k:n,k)))if a==0returnendfor i=k:nif abs(A1(i,k))==ay=A1(i,:)A1(i,:)=A1(k,:)A1(k,:)=yx=b1(i,:)b1(i,:)=b1(k,:)b1(k,:)=xbreakendendif A1(k,k)~=0A1(k+1:n,k)=A1(k+1:n,k)/A1(k,k)A1(k+1:n,k+1:n)=A1(k+1:n,k+1:n)-A1(k+1:n,k)*A1(k,k+1:n) elsebreakendendL=tril(A1,0);for i=1:nL(i,i)=1;endLU=triu(A1,0)y1=L\b1x1=U\y1得到如下结果：x1 =3.845714853511634e+001.609517394778522e+00-1.547605454206655e+011.041130489899787e+01将A2=[10,-7,0,1;-3,2.0999********,6,2;5,-1,5,-1;0,1,0,2]b2=[8;5.900000000001;5;1]代入上述代码求得结果如下：X2 =4.440892098500626e-16-9.999999999999993e-019.999999999999997e-011.000000000000000e+00(3)用不选主元的高斯消去法求得L和U及解向量x1,x2∈R4代码：format longeformat compactA1=[0.3*10^-15,59.14,3,1;5.291,-6.130,-1,2;11.2,9,5,2;1,2,1,1] b1=[59.17;46.78;1;2][L,U]=lu(A1)y1=L\b1x1=U\y1求得如下结果：x1=3.845714853511634e+001.609517394778522e+00-1.547605454206655e+011.041130489899787e+01将A2=[10,-7,0,1;-3,2.0999********,6,2;5,-1,5,-1;0,1,0,2] b2=[8;5.900000000001;5;1]代入上述代码，求得结果如下：x 2 =4.440892098500626e-16 -9.999999999999993e-01 9.999999999999997e-01 9.999999999999999e-01（2）（3）求得结果相同，可验证结果正确。
数值分析第一次大作业
《数值分析》计算实习报告第一题院系：机械工程及自动化学院_学号： _____姓名： _ ______2017年11月7日一、算法设计方案1、求λ1，λ501和λs 的值1)利用幂法计算出矩阵A 按模最大的特征值，设其为λm 。
2)令矩阵B =A −λm I （I 为单位矩阵），同样利用幂法计算出矩阵B 按模最大的特征值λm ′。
3)令λm ′′=λm ′+λm 。
由计算过程可知λm 和λm ′′分别为矩阵A 所有特征值按大小排序后，序列两端的值。
即，λ1=min{λm ,λm ′′}，λ501=max{λm ,λm ′′}。
4) 利用反幂法计算λs 。
其中，反幂法每迭代一次都要求解线性方程组1k k Au y -=，由于矩阵A 为带状矩阵，故可用三角分解法解带状线性方程组的方法求解得到k u 。
2、求A 的与数μk =λ1+k λ501−λ140最接近的特征值λi k (k =1,2, (39)1) 令矩阵D k =A −μk I ，利用反幂法计算出矩阵D k 按模最小的特征值λi k ′,则λi k =λi k ′+μk 。
3、求A 的（谱范数）条件数cond(A )2和行列式det A1) cond(A)2=|λm λs |，前文已算出m λ和s λ，直接带入即可。
2) 反幂法计算λs 时，已经对矩阵A 进行过Doolittle 分解，得到A=LU 。
而L 为对角线上元素全为1的下三角矩阵，U 为上三角矩阵，可知det 1L =，5011det ii i U u ==∏，即有5011det det det ii i A L U u ====∏。
最后，为节省存储量，需对矩阵A 进行压缩，将A 中带内元素存储为数组C [5][501]。
二、源程序代码#include<windows.h>#include<iostream>#include<iomanip>#include<math.h>using namespace std;#define N 501#define K 39#define r 2#define s 2#define EPSI 1.0e-12//求两个整数中的最大值int int_max2(int a, int b){return(a>b ? a : b);}//求两个整数中的最小值int int_min2(int a, int b){return(a<b ? a : b);}//求三个整数中的最大值int int_max3(int a, int b, int c){int t;if (a>b)t = a;else t = b;if (t<c) t = c;return(t);}//定义向量内积double dianji(double x[], double y[]) {double sum = 0;for (int i = 0; i<N; i++)sum = sum + x[i] * y[i];return(sum);}//计算两个数之间的相对误差double erro(double lamd0, double lamd1){double e, d, l;e = fabs(lamd1 - lamd0);d = fabs(lamd1);l = e / d;return(l);}//矩阵A的压缩存储初始化成Cvoid init_c(double c[][N]){int i, j;for (i = 0; i<r + s + 1; i++)for (j = 0; j<N; j++)if (i == 0 || i == 4)c[i][j] = -0.064;else if (i == 1 || i == 3)c[i][j] = 0.16;elsec[i][j] = (1.64 - 0.024*(j + 1))*sin(0.2*(j + 1)) - 0.64*exp(0.1 / (j + 1)); }//矩阵复制void fuzhi_c(double c_const[][N], double c[][N]){int i, j;for (i = 0; i<r + s + 1; i++)for (j = 0; j<N; j++)c[i][j] = c_const[i][j];}//LU三角分解void LUDet_c(double c_const[][N], double c_LU[][N]){double sum;int k, i, j;fuzhi_c(c_const, c_LU);for (k = 1; k <= N; k++){for (j = k; j <= int_min2(k + s, N); j++){sum = 0;for (i = int_max3(1, k - r, j - s); i <= k - 1; i++)sum += c_LU[k - i + s][i - 1] * c_LU[i - j + s][j - 1];c_LU[k - j + s][j - 1] -= sum;}for (j = k + 1; j <= int_min2(k + r, N); j++){sum = 0;for (i = int_max3(1, j - r, k - s); i <= k - 1; i++)sum += c_LU[j - i + s][i - 1] * c_LU[i - k + s][k - 1];c_LU[j - k + s][k - 1] = (c_LU[j - k + s][k - 1] - sum) / c_LU[s][k - 1];}}}//三角分解法解带状线性方程组void jiefc(double c_const[][N], double b_const[], double x[]){int i, j;double b[N], c_LU[r + s + 1][N], sum;for (i = 0; i<N; i++)b[i] = b_const[i];LUDet_c(c_const, c_LU);for (i = 2; i <= N; i++){sum = 0;for (j = int_max2(i - 2, 1); j <= i - 1; j++)sum += c_LU[i - j + 2][j - 1] * b[j - 1];b[i - 1] -= sum;}x[N - 1] = b[N - 1] / c_LU[2][N - 1];for (i = N - 1; i >= 1; i--){sum = 0;for (j = i + 1; j <= int_min2(i + 2, N); j++)sum += c_LU[i - j + 2][j - 1] * x[j - 1];x[i - 1] = (b[i - 1] - sum) / c_LU[2][i - 1];}}//幂法求按模最大特征值double mifa_c(double c_const[][N]){double u[N], y[N];double sum, length_u, beta0, beta1;int i, j;for (i = 0; i<N; i++)//迭代初始向量u[i] = 0.5;length_u = sqrt(dianji(u, u));for (i = 0; i<N; i++)y[i] = u[i] / length_u;for (i = 1; i <= N; i++){sum = 0;for (j = int_max2(i - 2, 1); j <= int_min2(i + 2, N); j++)sum = sum + c_const[i - j + 2][j - 1] * y[j - 1];u[i - 1] = sum;}beta1 = dianji(u, y);do{beta0 = beta1;length_u = sqrt(dianji(u, u));for (i = 0; i<N; i++)y[i] = u[i] / length_u;for (i = 1; i <= N; i++){sum = 0;for (j = int_max2(i - 2, 1); j <= int_min2(i + 2, N); j++)sum = sum + c_const[i - j + 2][j - 1] * y[j - 1];u[i - 1] = sum;}beta1 = dianji(u, y);} while (erro(beta0, beta1) >= EPSI);return(beta1);}//反幂法求按模最小特征值double fmifa_c(double c_const[][N]){double u[N], y[N];double length_u, beta0, beta1;int i;for (i = 0; i<N; i++)//迭代初始向量u[i] = 0.5;length_u = sqrt(dianji(u, u));for (i = 0; i<N; i++)y[i] = u[i] / length_u;jiefc(c_const, y, u);beta1 = dianji(y, u);do{beta0 = beta1;length_u = sqrt(dianji(u, u));for (i = 0; i<N; i++)y[i] = u[i] / length_u;jiefc(c_const, y, u);beta1 = dianji(y, u);} while (erro(beta0, beta1) >= EPSI);beta1 = 1 / beta1;return(beta1);}//计算lamd_1、lamd_501、lamd_svoid calculate1(double c_const[][N], double &lamd_1, double &lamd_501, double &lamd_s) {int i;double lamd_mifa0, lamd_mifa1, c[r + s + 1][N];lamd_mifa0 = mifa_c(c_const);fuzhi_c(c_const, c);for (i = 0; i<N; i++)c[2][i] = c[2][i] - lamd_mifa0;lamd_mifa1 = mifa_c(c) + lamd_mifa0;if (lamd_mifa0<lamd_mifa1){lamd_1 = lamd_mifa0;lamd_501 = lamd_mifa1;}else{lamd_501 = lamd_mifa0;lamd_1 = lamd_mifa1;}lamd_s = fmifa_c(c_const);}//平移+反幂法求最接近u_k的特征值void calculate2(double c_const[][N], double lamd_1, double lamd_501, double lamd_k[]){int i, k;double c[r + s + 1][N], h, temp;temp = (lamd_501 - lamd_1) / 40;for (k = 1; k <= K; k++){h = lamd_1 + k*temp;fuzhi_c(c_const, c);for (i = 0; i<N; i++)c[2][i] = c[2][i] - h;lamd_k[k - 1] = fmifa_c(c) + h;}}//计算cond(A)和det(A)void calculate3(double c_const[][N], double lamd_1, double lamd_501, double lamd_s, double &cond_A, double &det_A){int i;double c_LU[r + s + 1][N];if (fabs(lamd_1)>fabs(lamd_501))cond_A = fabs(lamd_1 / lamd_s);elsecond_A = fabs(lamd_501 / lamd_s);LUDet_c(c_const, c_LU);det_A = 1;for (i = 0; i<N; i++)det_A *= c_LU[2][i];}//*主程序*//int main(){int i, count = 0;double c_const[5][N], lamd_k[K];double lamd_1, lamd_501, lamd_s;double cond_A, det_A;//设置白背景黑字system("Color f0");//矩阵A压缩存储到c[5][501]init_c(c_const);cout << setiosflags(ios::scientific) << setiosflags(ios::right) << setprecision(12) << endl;//计算lamd_1、lamd_501、lamd_scalculate1(c_const, lamd_1, lamd_501, lamd_s);cout << " 矩阵A的最小特征值：λ1 = " << setw(20) << lamd_1 << endl;cout << " 矩阵A的最大特征值：λ501 = " << setw(20) << lamd_501 << endl;cout << " 矩阵A的按模最小的特征值：λs = " << setw(20) << lamd_s << endl;//求最接近u_k的特征值calculate2(c_const, lamd_1, lamd_501, lamd_k);cout << endl << " 与数u_k最接近的特征值：" << endl;for (i = 0; i<K; i++){cout << " λ_ik_" << setw(2) << i + 1 << " = " << setw(20) << lamd_k[i] << " ";count++;if (count == 2){cout << endl;count = 0;}}//计算cond_A和det_Acalculate3(c_const, lamd_1, lamd_501, lamd_s, cond_A, det_A);cout << endl << endl;cout << " 矩阵A的条件数：cond(A) = " << setw(20) << cond_A << endl;cout << " 矩阵A的行列式的值：det(A) = " << setw(20) << det_A << endl << endl;return 0;}三，计算结果四，分析初始向量选择对计算结果的影响当选取初始向量0(1,1,,1)Tu=时，计算的结果如下：此结果即为上文中的正确计算结果。
数值分析报告第二题 北航 大作业
实用文档数值分析第二题目录数值分析第二题 (1)1 引言： (1)1.1 矩阵的拟上三角化 (1)1.2 矩阵的特征值求解 (2)1.3 矩阵的特征向量求解 (4)2 算法的程序实现 (5)2.1 主程序 (5)2.2 子程序的实现 (7)3计算结果 (8)3.2矩阵Q、R 以及乘积RQ (9)3.3 各实特征值及其相对应的特征向量 (10)4 实验结论 (12)附录源代码 (12)1 引言1.1 矩阵的拟上三角化为了减少计算量，对矩阵A 利用Householder 矩阵进行相似变换，把A 化为上三角矩阵 。
对A 拟上三角化，得到拟上三角矩阵)1(-n 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 -=ωT rr T r r r r p u u A A --=+ω)()1(5. 继续。
1.2 矩阵的特征值求解使用带双步位移的QR 方法计算矩阵)1(-n A 的全部特征值，也是A 的全部特征值，具体算法如下：1. 给定精度水平0>ε和迭代最大次数L 。
北航研究生数值分析编程大作业1
数值分析大作业一、算法设计方案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、利用反幂法求出ik s λλ,反幂法迭代格式：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 λ。
北理工数值分析大作业
数值分析上机作业第 1 章1.1计算积分，n=9。
（要求计算结果具有6位有效数字）程序：n=1:19;I=zeros(1,19);I(19)=1/2*((exp(-1)/20)+(1/20));I(18)=1/2*((exp(-1)/19)+(1/19));for i=2:10I(19-i)=1/(20-i)*(1-I(20-i));endformat longdisp(I(1:19))结果截图及分析：在MATLAB中运行以上代码，得到结果如下图所示：当计算到数列的第10项时，所得的结果即为n=9时的准确积分值。
取6位有效数字可得.1.2分别将区间[-10.10]分为100,200,400等份，利用mesh或surf命令画出二元函数z=的三维图形。
程序：>> x = -10:0.1:10;y = -10:0.1:10;[X,Y] = meshgrid(x,y);Z = exp(-abs(X))+cos(X+Y)+1./(X.^2+Y.^2+1);subplot(2,2,1);mesh(X,Y,Z);title('步长0.1')>> x = -10:0.2:10;y = -10:0.2:10;[X,Y] = meshgrid(x,y);Z = exp(-abs(X))+cos(X+Y)+1./(X.^2+Y.^2+1);subplot(2,2,1);mesh(X,Y,Z);title('步长 0.2')>>x = -10:0.05:10;y = -10:0.05:10;[X,Y] = meshgrid(x,y);Z = exp(-abs(X))+cos(X+Y)+1./(X.^2+Y.^2+1);subplot(2,2,1);mesh(X,Y,Z);title('步长0.05')结果截图及分析：由图可知，步长越小时，绘得的图形越精确。
数值分析期末大作业
一、问题提出设方程f(x)=x 3-3x-1=0有三个实根 x *1=1.8793 , x *2=-0.34727 ,x *3=-1.53209现采用下面六种不同计算格式，求 f(x)=0的根 x *1 或x *2 。
1、 x = 213xx + 2、x = 313-x3、 x = 313+x4、 x = 312-x 5、 x = x13+6、 x = x - ()1133123---x x x二、目的和意义1、通过实验进一步了解方程求根的算法；2、认识选择计算格式的重要性；3、掌握迭代算法和精度控制；4、明确迭代收敛性与初值选取的关系。
三、结构程序设计本程序实在matlab 软件上进行操作的。
首先建立一个空白的M-文件。
在编辑器中输入以下内容，并保存。
function [X1,m,n,q]=shizi1(p) x=zeros(100,1); x=double(x);x(1,1)=p;i=1;deltax=100;while (i<100 & deltax > 0.000001)x(i+1,1)=(3*x(i,1)+1)/x(i,1)^2deltax=abs(x(i+1,1)-x(i,1));i=i+1;endX1=x(1,1);m=i;n=x(i,1);q=deltax;以上是运行函数，下一步在建立一个执行M-文件，输入以下内容，并保存。
其中X1为初始值，m为迭代次数，n为最后得到的值，q为|x k+1-x k|。
clear all;clc;p=1.8;[X1,m,n,q]=shizi1(p)1、对第一个迭代公式，在执行文件中输入p=1.8；[X1,m,n,q]=shizi1(p)。
得到如下结果如下：初值为1.8，迭代100次，精度为10-6。
可见该迭代公式是发散的，将初值改为-1.5，其他均条件不变。
p=-1.5;[X1,m,n,q]=shizi1(p)改变初值后可以得到一个接近真值的结果x*3的结果ans=-1.5321。
北航数值分析大作业一
北航数值分析大作业一————————————————————————————————作者: ————————————————————————————————日期：ﻩ《数值分析B》大作业一ＳY110３1２0 朱舜杰一．算法设计方案：1.矩阵A的存储与检索将带状线性矩阵A［501］[5０1]转存为一个矩阵MatｒixC[5][501] . 由于C语言中数组角标都是从0开始的,所以在数组MatｒixＣ[５][501]中检索A的带内元素a ij的方法是：Ａ的带内元素aｉｊ＝C中的元素ｃｉ-j＋２，ｊ2.求解λ１,λ5０1,λｓ①首先分别使用幂法和反幂法迭代求出矩阵按摸最大和最小的特征值λmaｘ和λmｉｎ。
λmｉn即为λs;如果λmａx>0，则λ５01=λmax;如果λmax<０，则λ１=λmaｘ。
②使用带原点平移的幂法（mｉfa(）函数)，令平移量p=λmａx,求出对应的按摸最大的特征值λ，max,如果λmaｘ>０,则λ1=λ，mａｘ+p；如果λｍax<０,则λ５01=λ，mａx+p。
3.求解A的与数μk＝λ1+k（λ50１-λ1）/4０的最接近的特征值λiｋ（ｋ=1，2,…,3９)。
使用带原点平移的反幂法,令平移量ｐ＝μｋ,即可求出与μk最接近的特征值λｉk。
4.求解A的(谱范数)条件数cｏnd（A）２和行列式detA。
①ｃond（A）2=｜λ1/λn|，其中λ１和λn分别是矩阵Ａ的模最大和最小特征值。
②矩阵A的行列式可先对矩阵A进行LU分解后,ｄetＡ等于U所有对角线上元素的乘积。
二．源程序#includｅ＜ｓtdio.ｈ>＃incｌｕde<iostrｅaｍ.h>#incluｄｅ<ｓｔdｌib.h>#inｃlude<ｍath．h>#iｎclude<flｏat.ｈ＞#include<iomanip.h>＃inｃlude<tiｍｅ.h>#define E １.0e-１2/*定义全局变量相对误差限*/int maｘ２(intａ，ｉｎｔb）ﻩ／＊求两个整型数最大值的子程序*／｛ﻩif(a>b)retｕrn a；elｓeretｕrn b;｝int mｉn2(iｎt a,ｉnｔｂ) /*求两个整型数最小值的子程序＊／{iｆ(a>ｂ)rｅtuｒn b；eｌsereｔurn ａ;｝int mａx3(iｎt ａ，ｉnt ｂ,ｉｎｔc)／＊求三整型数最大值的子程序＊/｛int t;ﻩif(a>b）ﻩt=a;ﻩelsｅt=ｂ;ﻩif(t＜c) ｔ＝c;reｔｕｒｎ（t)；}voiｄasｓignment（dｏuble aｒray[5][5０1]) /*将矩阵A转存为数组Ｃ[５][501]*／{iｎt ｉ，j,ｋ;/／所有元素归零fｏｒ(i＝0;i<=4;){for(j=0;j<=5００;){arｒay[i][j]＝0;ﻩj++;}ﻩ i+＋;}//第0,4行赋值fｏr(ｊ＝２;j＜=５00;){ﻩk＝500-ｊ;arraｙ[0][ｊ]=－０.06４;arraｙ[4][k］=-0．064；ﻩｊ++；｝//第１,３行赋值ｆor(j＝1;ｊ<=５0０；)｛ﻩﻩｋ＝500-j;array[1][j]＝0.１6;arraｙ［3][ｋ]=0.1６;ﻩ j++；}//第2行赋值fｏr(j＝0;j<＝5０0;){k=j;ﻩﻩj++;arraｙ［2][k]＝(１．６4－０．０2４*j)*sin((doublｅ)(0．2*ｊ))-０．6４*exp((do ｕble)(0.1／ｊ）)；}}doublｅｍifa(douｂle u[５01］,ｄouble arraｙ[5][5０1］,dｏuｂle p) ／*带原点平移的幂法*/{int ｉ,j; /* u[５01]为初始迭代向量*/double a，ｂ，c=０; ／* arｒaｙ[5][５01]为矩阵A的转存矩阵*／ﻩdｏｕｂle y[501]; /＊ｐ为平移量＊／ﻩfor(；;)｛a=0；ﻩｂ=０；/*选用第一种迭代格式*／/／求ηk-1fｏｒ(ｉ=０；ｉ<=50０；ｉ＋+){ﻩa=a+u[i]*u[i］;ﻩ}ａ=sqrｔ(a);//求y k-1ｆoｒ(i=0;i<=5０0;i++){ﻩﻩｙ[i］=u[i］/a；}/／求u kｆor(i＝0;i<=500;i++)ﻩ{ﻩu［ｉ］=０；foｒ(j＝max2(i-2,0)；j<=min2(i+２,50０）;ｊ+＋）ﻩ{ｕ[i］＋=arrａy[ｉ-j+2][j]*y[j];ﻩ}ﻩｕ[i]=u[i]－p*ｙ[i］; /*引入平移量＊/ﻩ}/／求βkﻩfor（i=0;i<=50０；i++)ﻩ{ﻩﻩb+＝ｙ［ｉ]*u[i］;ﻩ｝ﻩif(faｂs(（b-c)/b)<=Ｅ) /＊达到精度水平,迭代终止＊/ ﻩｂreak;c＝b；ﻩ}returｎ（b+p)；／*直接返回A的特征值*/｝ｖoiｄchuｚhi(doｕble a[]) /＊用随机数为初始迭代向量赋值＊／{ｉｎｔi；srａnd((ｉnt)time（0));ｆｏr（i=0;i<=500;i+＋)ﻩ{ﻩﻩa[i]=(１0.0*rand（）/ＲAND＿MAＸ)；/*生成0~1０的随机数*/ }}void chｕzｈi2(double a[],ｉnt j）/*令初始迭代向量为e i*/{int i；ﻩfｏr（ｉ=0;i<=500；i++）ﻩ{a[i]＝0;ﻩ｝a[j]=1;}voｉd LU(dｏｕble ａrrａｙ[5][501]) /＊对矩阵A进行Doolittle分解＊／｛／*矩阵A转存在C［5][５0１]中*／ﻩint j,k,t;/*分解结果L，Ｕ分别存在C[５］[501]的上半部与下半部*/ﻩfｏｒ(k=0；k<=5０0;k+＋)ﻩ{ﻩﻩfor(j=k;j<＝min２(（k+2），500）;ｊ+＋)ﻩﻩ{ﻩﻩfor(t=maｘ3(0,ｋ-2，j－2)；ｔ＜=(k-1)；ｔ++)ﻩ｛ﻩﻩaｒraｙ[ｋ－ｊ+2]［j]－=ａrraｙ［ｋ－t+2］[ｔ]*aｒraｙ[ｔ-j+2]［j]; ﻩﻩ}ﻩ}ﻩiｆ(ｋ<５00)ﻩfor(j=ｋ+1;j＜=min２(（k+2),5００）;j+＋){ﻩﻩfｏr（ｔ＝maｘ3(0，k－2，ｊ-2);t<=(k－1)；t++){ﻩﻩarray[j-ｋ+２][k]-＝aｒrａｙ[j-t＋２][t]*arｒay［t－ｋ+2][k]；ﻩﻩ}array[j-ｋ+2]［ｋ]=array[j-k+2][k]／arｒay[2][ｋ];}ﻩ}}dｏuｂｌe fmｉfa(double u[5０１],double array［５][501],dｏｕｂle p）{ /*带原点平移的反幂法*/ﻩｉｎｔi，j;ﻩdｏｕbｌｅa,b，c=0;ﻩdｏublｅｙ［501］;//引入平移量fｏr(ｉ=0；i<=50０;i＋＋)ﻩ{ａｒrａy［2]［i]－＝p；ﻩ}//先将矩阵Dｏｏlittlｅ分解LＵ(array);ﻩｆｏr(;;){a=0;b=0;／/求ηk-1foｒ（ｉ＝０;i<＝５00;i++){ﻩﻩa＝ａ+u[i］*ｕ［i］;}a=sqrt(a）;/／求ｙk－1ﻩfor(i＝0;ｉ<=50０;ｉ＋+)ﻩ{ｙ[i]=u[i]/a;ﻩ}//回带过程,求解u kfor(i=0;i＜＝500;i++)ﻩ{ﻩu[i］=y[i];ﻩ｝for（i=1;i＜＝５00;i++）ﻩ｛ﻩfor(j=ｍaｘ2（0,(i－２）);j＜＝(ｉ-1）;j+＋)ﻩ{ﻩﻩu[i]－=arrａy［ｉ-j+2][ｊ]*u[j];ﻩ}ﻩﻩ}ﻩu[50０]=u[5０0]/ａｒray[2][500];for（i=４９9；i>=0;i－－)ﻩ｛foｒ(ｊ=i+1;ｊ＜=min２((i＋2）,５00）;j+＋)ﻩﻩ｛ﻩﻩu［ｉ]－＝aｒraｙ[i-j+2][j]＊u[j］;ﻩ}ﻩﻩu[ｉ]=u[ｉ］／ａrrａy［２][ｉ];ﻩ}／／求βｋfor(ｉ=0;i＜=500；i＋+){ﻩﻩb+=y[i]*u［ｉ];}if（faｂｓ((b-c)/b)<=E) /*达到精度要求,迭代终止*/bｒeａk;ﻩc＝b;}reｔurｎ(ｐ＋(1／ｂ)）; ／＊直接返回距离原点Ｐ最接近的A的特征值*／｝//主函数ｍain(){ int i;douｂｌe ｄ1，d501,dｓ,d,a；ﻩdouｂle u[５0１];ﻩdouｂlｅMatrｉxC[５][501];prinｔｆ(" 《数值分析》计算实习题目第一题\n＂);ﻩprintｆ（＂SＹ１1０312０朱舜杰\n＂);／/将矩阵A转存为MatrixＣassiｇｎment(MatrixＣ);/／用带原点平移的幂法求解λ1,λ５01chｕzhi(u）;ﻩd=mifａ（u,MaｔrixＣ,0）;chuzhｉ(u）；ﻩa=ｍifa(u，MatｒｉxC,d）；ﻩiｆ(d<0)ﻩ｛ﻩﻩｄ1=d;ﻩｄ５０1=a；ﻩ}elseﻩ｛ﻩﻩﻩｄ50１=d;ﻩﻩd1=a;ﻩ｝prｉnｔf(＂λ1=%.12e\ｎ＂,ｄ1）;ﻩpｒinｔｆ(＂λ５01=%.1２ｅ\n"，d５０1）;／/用反幂法求λｓchuｚhi(u);ｄs=fmｉfa(u,MaｔrixC,０）;ｐrintf("λs＝%.１2e＼ｎ",ds)；//用带原点平移的反幂法求λikfｏr(i=1;i＜＝3９;i++){ﻩa=d1＋（i＊（ｄ501-d1）)/40;ﻩﻩasｓigｎment(MatrixC）;ﻩｃhｕzhi(u）;ﻩﻩd=fｍifa（u,MａtriｘC，a);ﻩﻩprinｔｆ("与μ%02ｄ=%＋.12e最接近的特征值λｉ%02d=%+．1２e\n",i,a,i,d); ﻩ}//求A的条件数d=fabｓ((d1/ds))；ﻩpriｎtf("Ａ的(谱范数）条件数coｎｄ＜A>2=%.12e\n",ｄ)；//求ｄetAﻩaｓsignment（ＭaｔｒｉｘC）;LU(MatrixC);ﻩa=1;fｏr(i=0;i＜=５00;ｉ+＋）{a*=MatrixＣ[2][i];｝ｐrintｆ（"行列式deｔＡ＝％.12e\n",a);/／测试不同迭代初始向量对λ１计算结果的影响。
昆明理工数值分析大作业插值法数值分析作业
5 / 11
y(i)=s; end
%所得近似值给对应返回值
问题一: 求解及画图程序： xj=[0.4 0.55 0.65 0.8 0.95 1.05]; yj=[0.41075 0.57815 0.69675 0.9 1 1.25382]; xs=(0.3:0.01:1.1); y0=lagrange(xj,yj,xs); plot(xs,y0) >> hold on xs=[0.596 0.99]; y0=lagrange(xj,yj,xs) y0 = 0.6257 1.0542 求解结果; x y 问题一插值函数绘图： * 所求点 插值节点 0.596 0.6257
function y=newton(f,x0,x); m=length(x);n=length(x0); for i=1:m p=zeros(1,n-1); suma=zeros(1,n); suma(1)=f(1); sumax=0; for k=1:(n-1); p(k)=x(i)-x0(k); end for j=2:n; q=1; for s=1:(j-1); q=q*p(s); end suma(j)=f(j)*q; end for h=1:n; sumax=sumax+suma(h); end y(i)=sumax; end
������ =0 ������ =0 ������≠������
������ − ������������ ������������ − ������������
即：
������������ ������ =
������ ������ =0
������������
������ ������ +1 ������ ������−������ ������ ������ ’������ +1 ������ ������
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
数值分析大作业学号：*********专业：机械工程学生姓名：***2014年10月摘要：在自然科学与工程技术中，很多问题的解决常常归结为求解线性方程组Ax=b 。
随着计算机的发展，利用计算机这个强有力的计算工具去求解线性方程组是一个非常实用的问题。
在求解大型线性方程组时，直接法在多次消元，回代的过程中，四则运算的误差累计与传播无法控制，致使计算结果的精度就无法保证，特别是求解大型稀松矩阵时，还要对系数矩阵进行分解。
而迭代法相对于直接法而言，具有保持迭代矩阵不变的特点，计算程序一般也比较简单，且对于许多问题收敛速度比较快。
比较常用的迭代法有雅克比迭代法、高斯一塞德尔迭代法和逐次超松弛迭代法等，本次研究目的是通过求解一个线性方程组来比较它们的迭代效果，验证一些已有的结论。
1．数学原理1.1雅可比迭代法将线性方程组的系数矩n n ij R a A ⨯∈=)(分解为A=D+L+U ，其中D 是由A 的主对角元素构成的对角矩阵，L 是由A 的严格下三角部分构成的严格下三角矩阵， U 是由A 的严格上三角部分构成的严格上三角矩阵，即,2211⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=nn a a a D.0000,0000,1223113121,21323121⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=--n n n n n n n n a a a a a a U a a a a a a L若系数矩阵A 的对角元素),,2,1(0n i a ii =≠，则矩阵D 非奇异，取M=D ，N=-（L+U)，则J J g x G b D x U L D x +=++-=--11)(， 因而，构造的迭代法为：.),(,11)()1(b D g U L D G g x G x J J J k J k --+=+-=+=1.2高斯-赛得尔迭代法将线性方程组的系数矩n n ij R a A ⨯∈=)(分解为A=D+L+U 。
若系数矩阵A 的对角元素不等于0，则矩阵D 非奇异，取M=L+D ，N=-U ，则()()G G g x G b D L Ux D L x +=+++-=--11因而，构造的迭代法为：()().,,11)()1(b D L g U D L G g x G x G G G k G k --++=+-=+=1.3逐次超松弛迭代法线性方程组的系数矩n n ij R a A ⨯∈=)(分解为A=D+L+U 。
取()()[],11,1U D N L D M ωωωωω--=+=则()()[]()S S g x G b L D x U D L D x +=++--+=--111ωωωωω因而，构造的迭代法为：()()[]().,1,11)()1(b L D g U D L D G g x G x S G S k S k --++=--+=+=ωωωωω2.程序设计2.1雅可比迭代法求解雅可比迭代法MATLAB 程序如下： %majacobi.mfunction x=majacobi(A,b,x0,ep,N)%用途：用Jacobi 迭代法解线性方程组Ax=b%格式;x=majacobi(A,b,x0,ep,N) A 为系数矩阵，b 为右端向量， %x0为初始向量（默认零向量），ep 为精度（默认1e-6), %N 为最大迭代次数（默认500次），x 返回近似解向量 n=length(b);if nargin<5,N=500;end if nargin<4,ep=1e-6;endif nargin<3,x0=zeros(n,1);end x=zeros(n,1);k=0; while k<N for i=1:nx(i)=(b(i)-A(i,[1:i-1,i+1:n])*x0([1:i-1,i+1:n]))/A(i,i); endif norm(x-x0,inf)<ep,break; end x0=x ;k=k+1; endif k==N,Warning('100');end disp(['k=',num2str(k)])在MATLAB 命令窗口执行创建的m 文件 >> edit majacobi>> A=[10.9,1.2,2.1,0.9;1.2,11.2,1.5,2.5;2.1,1.5,9.8,1.3;0.9,2.5,1.3,12.3]A =10.9000 1.2000 2.1000 0.90001.2000 11.2000 1.50002.50002.1000 1.5000 9.8000 1.30000.9000 2.5000 1.3000 12.3000>> b=[-7.0,5.3,10.3,24.6]';>> x=majacobi(A,b)k=17x =-0.99860.00711.00321.96562.2高斯-赛得尔迭代法求解高斯-赛得尔迭代法MATLAB程序如下：%maseidel.mfunction x=maseidel (A,b,x0,ep,N)%用途：用Gauss-Seidel迭代法解线性方程组Ax=b%格式：x=maseidel （A,b,x0,ep,N) A为系数矩阵，b为右端向量，%x0为初始向量（默认零向量），ep为精度（默认1e-6),%N为最大迭代次数（默认500次），x返回近似解向量n=length(b);if nargin<5,N=500;endif nargin<4,ep=1e-6;endif nargin<3,x0=zeros(n,1);endx=zeros(n,1);k=0;while k<Nfor i=1:nif i==1x(1)=(b(1)-A(1,2:n)*x0(2:n))/A(1,1);else if i==nx(n)=(b(n)-A(n,1:n-1)*x(1:n-1))/A(n,n);elsex(i)=(b(i)-A(i,1:i-1)*x(1:i-1)-A(i,i+1:n)*x0(i+1:n))/A(i,i);endendendif norm(x-x0,inf)<ep,break;endx0=x;k=k+1;endif k==N,Warning('100');enddisp(['k=',num2str(k)])在MATLAB命令窗口执行创建的m文件：>> edit maseidel>> A=[10.9,1.2,2.1,0.9;1.2,11.2,1.5,2.5;2.1,1.5,9.8,1.3;0.9,2.5,1.3,12.3]A =10.9000 1.2000 2.1000 0.90001.2000 11.2000 1.50002.50002.1000 1.5000 9.8000 1.30000.9000 2.5000 1.3000 12.3000>> b=[-7.0,5.3,10.3,24.6]';>> x=maseidel(A,b)k=7x =-0.99860.00711.00321.96562.3逐次超松弛迭代法求解逐次超松弛迭代法MATLAB程序如下：%masor.mfunction x=masor(A,b,omega,x0,ep,N)%用途：用SOR迭代法解线性方程组Ax=b%格式：x=maseidel （A,b,x0,ep,N) A为系数矩阵，b为右端向量，%omega为松弛因子（默认1.5),x0为初始向量（默认零向量），ep为精度（默认1e-6),%N为最大迭代次数（默认500次），x返回近似解向量n=length(b);if nargin<6,N=500;endif nargin<5,ep=1e-6;endif nargin<4,x0=zeros(n,1);endif nargin<3,omega=1.5;endx=zeros(n,1);k=0;while k<Nfor i=1:nif i==1x1(1)=(b(1)-A(1,2:n)*x0(2:n))/A(1,1);else if i==nx1(n)=(b(n)-A(n,1:n-1)*x(1:n-1))/A(n,n);elsex1(i)=(b(i)-A(i,1:i-1)*x(1:i-1)-A(i,i+1:n)*x0(i+1:n))/A(i,i);endendx(i)=(1-omega)*x0(i)+omega*x1(i);endif norm(x0-x,inf)<ep,break;endk=k+1;x0=x;endif k==N,Warning('100');enddisp(['k=',num2str(k)])取不同的松弛因子进行试验，在MATLAB命令窗口执行创建的m文件：>> edit masor>> A=[10.9,1.2,2.1,0.9;1.2,11.2,1.5,2.5;2.1,1.5,9.8,1.3;0.9,2.5,1.3,12.3]A =10.9000 1.2000 2.1000 0.90001.2000 11.2000 1.50002.50002.1000 1.5000 9.8000 1.30000.9000 2.5000 1.3000 12.3000>> b=[-7.0,5.3,10.3,24.6]';>> x=masor(A,b,0.4)k=33x =-0.99860.00711.00321.9656取其它的松弛因子得：取松弛因子为0.8 >> x=masor(A,b,0.8) 得到k=12；x = -0.9986，0.0071，1.0032，1.9656。
取松弛因子为1.2 >> x=masor(A,b,1.2) 得到k=11；x =-0.9986，0.0071，1.0032，1.9656。
取松弛因子为1.6 >> x=masor(A,b,1.6) 得到k=33；x =-0.9986，0.0071，1.0032，1.9656。
三．结果分析和讨论由程序解得的结果可知：雅可比迭代法在第17次迭代达到收敛，高斯-塞德尔迭代法在第7次迭代达到收敛。
从此例可以看出，高斯-塞德尔迭代法比雅可比迭代法收敛速度快(达到同样的精度所需迭代次数少)。
由雅可比迭代公式可知，在迭代的每一步计算过程中使用的是旧的分量来计算，没有用到迭代出的新的分量。