北航数值分析第一次大作业(高斯gauss lu分解)
北航数值分析大作业 第一题 幂法与反幂法
数 值 分 析（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.要求计算矩阵的最⼤最⼩特征值，通过幂法求得模最⼤的特征值，进⾏⼀定判断即得所求结果；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的系数可能不等于零。
北航数值分析大作业第一题幂法与反幂法
《数值分析》计算实习题目第一题：1. 算法设计方案（1）1λ，501λ和s λ的值。
1)首先通过幂法求出按模最大的特征值λt1，然后根据λt1进行原点平移求出另一特征值λt2，比较两值大小，数值小的为所求最小特征值λ1，数值大的为是所求最大特征值λ501。
2)使用反幂法求λs ，其中需要解线性方程组。
因为A 为带状线性方程组，此处采用LU 分解法解带状方程组。
（2）与140k λλμλ-5011=+k 最接近的特征值λik 。
通过带有原点平移的反幂法求出与数k μ最接近的特征值 λik 。
（3）2cond(A)和det A 。
1）1=nλλ2cond(A)，其中1λ和n λ分别是按模最大和最小特征值。
2）利用步骤（1）中分解矩阵A 得出的LU 矩阵，L 为单位下三角阵,U 为上三角阵，其中U 矩阵的主对角线元素之积即为det A 。
由于A 的元素零元素较多，为节省储存量，将A 的元素存为6×501的数组中，程序中采用get_an_element()函数来从小数组中取出A 中的元素。
2.全部源程序#include <stdio.h>#include <math.h>void init_a();//初始化Adouble get_an_element(int,int);//取A 中的元素函数double powermethod(double);//原点平移的幂法double inversepowermethod(double);//原点平移的反幂法int presolve(double);//三角LU 分解int solve(double [],double []);//解方程组int max(int,int);int min(int,int);double (*u)[502]=new double[502][502];//上三角U 数组double (*l)[502]=new double[502][502];//单位下三角L 数组double a[6][502];//矩阵Aint main(){int i,k;double lambdat1,lambdat2,lambda1,lambda501,lambdas,mu[40],det;init_a();//初始化Alambdat1=powermethod(0);lambdat2=powermethod(lambdat1);lambda1=lambdat1<lambdat2?lambdat1:lambdat2;lambda501=lambdat1>lambdat2?lambdat1:lambdat2;presolve(0);lambdas=inversepowermethod(0);det=1;for(i=1;i<=501;i++)det=det*u[i][i];for (k=1;k<=39;k++){mu[k]=lambda1+k*(lambda501-lambda1)/40;presolve(mu[k]);lambda[k]=inversepowermethod(mu[k]);}printf("------------所有特征值如下------------\n");printf("λ=%1.11e λ=%1.11e\n",lambda1,lambda501);printf("λs=%1.11e\n",lambdas);printf("cond(A)=%1.11e\n",fabs(lambdat1/lambdas));printf("detA=%1.11e \n",det);for (k=1;k<=39;k++){printf("λi%d=%1.11e ",k,lambda[k]);if(k % 3==0) printf("\n");} delete []u;delete []l;//释放堆内存return 0;}void init_a()//初始化A{int i;for (i=3;i<=501;i++) a[1][i]=a[5][502-i]=-0.064;for (i=2;i<=501;i++) a[2][i]=a[4][502-i]=0.16;for (i=1;i<=501;i++) a[3][i]=(1.64-0.024*i)*sin(0.2*i)-0.64*exp(0.1/i); }double get_an_element(int i,int j)//从A中节省存储量的提取元素方法{if (fabs(i-j)<=2) return a[i-j+3][j];else return 0;}double powermethod(double offset)//幂法{int i,x1;double beta=0,prebeta=-1000,yita=0;for (i=1;i<=501;i++)u[i]=1,y[i]=0;//设置初始向量u[]for (int k=1;k<=10000;k++){yita=0;for (i=1;i<=501;i++) yita=sqrt(yita*yita+u[i]*u[i]);for (i=1;i<=501;i++) y[i]=u[i]/yita;for (x1=1;x1<=501;x1++){u[x1]=0;for (int x2=1;x2<=501;x2++)u[x1]=u[x1]+((x1==x2)?(get_an_element(x1,x2)-offset):get_an_element(x1,x2))*y[x2];} prebeta=beta;beta=0;for (i=1;i<=501;i++) beta=beta+ y[i]*u[i];if (fabs((prebeta-beta)/beta)<=1e-12) {printf("offset=%f lambda=%f err=%e k=%d\n",offset,(beta+offset),fabs((prebeta-beta)/beta),k);break;};//输出中间过程,包括偏移量,误差,迭代次数}return (beta+offset);}double inversepowermethod(double offset)//反幂法{int i;double u[502],y[502];double beta=0,prebeta=0,yita=0;for (i=1;i<=501;i++)u[i]=1,y[i]=0; //设置初始向量u[]for (int k=1;k<=10000;k++){yita=0;for (i=1;i<=501;i++) yita=sqrt(yita*yita+u[i]*u[i]);for (i=1;i<=501;i++) y[i]=u[i]/yita;solve(u,y);prebeta=beta;beta=0;for (i=1;i<=501;i++) beta=beta+ y[i]*u[i];beta=1/beta;if (fabs((prebeta-beta)/beta)<=1e-12) {printf("offset=%f lambda=%f err=%e k=%d\n",offset,(beta+offset),fabs((prebeta-beta)/beta),k);break;};//输出中间过程,包括偏移量,误差,迭代次数}return (beta+offset);}int presolve(double offset)//三角LU分解{int i,k,j,t;double sum;for (k=1;k<=501;k++)for (j=1;j<=501;j++){u[k][j]=l[k][j]=0;if (k==j) l[k][j]=1;} //初始化LU矩阵for (k=1;k<=501;k++){for (j=k;j<=min(k+2,501);j++){sum=0;for (t=max(1,max(k-2,j-2)) ; t<=(k-1) ; t++)sum=sum+l[k][t]*u[t][j];u[k][j]=((k==j)?(get_an_element(k,j)-offset):get_an_element(k,j))-sum;}if (k==501) continue;for (i=k+1;i<=min(k+2,501);i++){sum=0;for (t=max(1,max(i-2,k-2));t<=(k-1);t++)sum=sum+l[i][t]*u[t][k];l[i][k]=(((i==k)?(get_an_element(i,k)-offset):get_an_element(i,k))-sum)/u[k][k];}}return 0;}int solve(double x[],double b[])//解方程组{int i,t;double y[502];double sum;y[1]=b[1];for (i=2;i<=501;i++){sum=0;for (t=max(1,i-2);t<=i-1;t++)sum=sum+l[i][t]*y[t];y[i]=b[i]-sum;}x[501]=y[501]/u[501][501];for (i=500;i>=1;i--){sum=0;for (t=i+1;t<=min(i+2,501);t++)sum=sum+u[i][t]*x[t];x[i]=(y[i]-sum)/u[i][i];}return 0;}int max(int x,int y){return (x>y?x:y);}int min(int x,int y){return (x<y?x:y);}3.计算结果结果如下图所示:部分中间结果：给出了偏移量(offset)，误差(err)，迭代次数(k)4.讨论迭代初始向量的选取对计算结果的影响,并说明原因使用u[i]=1（i=1,2,...,501）作为初始向量进行迭代，可得出以上结果。
北航数值分析全部三次大作业
北航数值分析全部三次大作业第一次大作业是关于解线性方程组的数值方法。
我们被要求实现各种常用的线性方程组求解算法，例如高斯消元法、LU分解法和迭代法等。
我首先学习了这些算法的原理和实现方法，并借助Python编程语言编写了这些算法的代码。
在实验中，我们使用了不同规模和条件的线性方程组进行测试，并比较了不同算法的性能和精度。
通过这个作业，我深入了解了线性方程组求解的原理和方法，提高了我的编程和数值计算能力。
第二次大作业是关于数值积分的方法。
数值积分是数值分析中的重要内容，它可以用于计算曲线的长度、函数的面积以及求解微分方程等问题。
在这个作业中，我们需要实现不同的数值积分算法，例如矩形法、梯形法和辛普森法等。
我学习了这些算法的原理和实现方法，并使用Python编写了它们的代码。
在实验中，我们计算了不同函数的积分值，并对比了不同算法的精度和效率。
通过这个作业，我深入了解了数值积分的原理和方法，提高了我的编程和数学建模能力。
第三次大作业是关于常微分方程的数值解法。
常微分方程是数值分析中的核心内容之一，它可以用于描述众多物理、化学和生物现象。
在这个作业中，我们需要实现不同的常微分方程求解算法，例如欧拉法、龙格-库塔法和Adams法等。
我学习了这些算法的原理和实现方法，并使用Python编写了它们的代码。
在实验中，我们解决了一些具体的常微分方程问题，并比较了不同算法的精度和效率。
通过这个作业，我深入了解了常微分方程的原理和方法，提高了我的编程和问题求解能力。
总的来说，北航数值分析课程的三次大作业非常有挑战性，但也非常有意义。
通过这些作业，我在数值计算和编程方面得到了很大的提升，也更加深入地了解了数值分析的理论和方法。
虽然这些作业需要大量的时间和精力，但我相信这些努力将会对我未来的学习和工作产生积极的影响。
北航数值分析大作业一
北京航空航天大学数值分析大作业一学院名称自动化专业方向控制工程学号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 值。
数值分析实验报告
%消元过程
fori=k+1:n
m=A(i,k)/A(k,k);
forj=k+1:n
A(i,j)=A(i,j)-m*A(k,j);
end
b(i)=b(i)-m*b(k);
end
det=det*A(k,k);
end
det=det*A(n,n);
%回代过程
ifabs(A(n,n))<1e-10
flag='failure';return;
*x=(x0,x1….,xn),插值节点
*y=(y0,y1,…,yn);被插函数f(x)在插值节点处的函数值
*t求插值函数Pn(x)在t处的函数值
*返回值 插值函数Pn(x)在t处的函数值
*/
procedureNewton
forj=0to n
d1jyj;
endfor
forj=1to n
fori=j to n
[n,m]=size(A);nb=length(b)
%当方程组行与列的维数不相等时，停止计算，并输出出错信息
ifn~=m
error('The row and columns of matrix A must beepual!');
return;
end
%当方程组与右端项的维数不匹配时，停止计算，并输出错误信息
clear
fprintf('gauss-seidel迭代法')
x1_(1)=0;
x2_(1)=0;
x3_(1)=0;
fori=1:9
x1_(i+1)=7.2+0.1*x2_(i)+0.2*x3_(i);
数值分析实验报告---高斯消去法 LU分解法
数值分析上机报告①高斯消去法利用高斯消去法的matlab程序源代码：A=[10,-7,0,1;-3,2.099999,6,2;5,-1,5,-1;2,1,0,2];b=[8;5.900001;5;1];x=A\bc=det(A)上述程序中A表示计算实习题1中线性方程组的系数矩阵，b表示线性方程组右边的矩阵，x表示线性方程组的解。
C所输出的是系数矩阵A的行列式的值。
程序运行结果：②列主元的高斯消去法利用列主元的高斯消去法matlab程序源代码：首先建立一个gaussMethod.m的文件，用来实现列主元的消去方法。
function x=gaussMethod(A,b)%高斯列主元消去法，要求系数矩阵非奇异的，%n = size(A,1);if abs(det(A))<= 1e-8error('系数矩阵是奇异的');return;end%for k=1:nak = max(abs(A(k:n,k)));index = find(A(:,k)==ak);if length(index) == 0index = find(A(:,k)==-ak);end%交换列主元temp = A(index,:);A(index,:) = A(k,:);A(k,:) = temp;temp = b(index);b(index) = b(k); b(k) = temp;%消元过程for i=k+1:nm=A(i,k)/A(k,k);%消除列元素A(i,k+1:n)=A(i,k+1:n)-m*A(k,k+1:n);b(i)=b(i)-m*b(k);endend%回代过程x(n)=b(n)/A(n,n);for k=n-1:-1:1;x(k)=(b(k)-A(k,k+1:n)*x(k+1:n)')/A(k,k);endx=x';end然后调用gaussMethod函数，来实现列主元的高斯消去法。
数值分析第一次大作业
《数值分析》计算实习报告第一题院系：机械工程及自动化学院_学号： _____姓名： _ ______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(C语言详细注释)
《数值分析》计算实习作业《一》北航第一题 设有501501⨯的矩阵⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=501500499321a bc b a b cc b a b ccb a bc c b a b c b a A其中.064.0,16.0);501,2,1(64.0)2.0sin()024.064.1(1.0-===--=c b i e i i a i i 矩阵的特征值)501,,2,1( =i i λ满足||min ||,501150121i i s λλλλλ≤≤=<<<试求1. 5011，λλ和s λ的值2. 的与数4015011λλκλμ-+=k 最接近的特征值)39,,2,1( =K κλi3. 的（谱范数）条件数2)A (cond 和行列式A det 要求1. 算法的设计方案（A 的所有零元素都不能存储）2. 全部源程序（详细注释）。
变量为double ，精度-1210=ε，输出为e 型12位有效数字3. 特征值s 5011，，λλλ和)39,,2,1( =K κλi 以及A cond det ,)A (2的值4. 讨论迭代初始向量的选取对计算结果的影响，并说明原因解答：1. 算法设计对于s λ满足||min ||5011i i s λλ≤≤=，所以s λ是按模最小的特征值，直接运用反幂法可求得。
对于5011，λλ，一个是最大的特征值，一个是最小的特征值，不能确定两者的绝对值是否相等，因此必须首先假设||||5011λλ≠，然后运用幂法，看能否求得一个特征值，如果可以求得一个，证明A 是收敛的，求得的结果是正确的，然后对A 进行带原点平移的幂法，偏移量是前面求得的特征值，可以求得另一个特征值，最后比较这两个特征值，较大的特征值是501λ，较小的特征值就是1λ。
如果在假设的前提下，无法运用幂法求得按模最大的特征值，即此时A 不收敛，则需要将A 进行带原点平移的幂法，平移量可以选取1，再重复上述步骤即可求得两个特征值。
数值分析实验报告---高斯消去法 LU分解法
数值分析实验报告---高斯消去法 LU分解法实验一：高斯消去法一、实验目的1. 掌握高斯消去法的原理2. 用高斯消去法解线性方程组3. 分析误差二、实验原理高斯消去法(又称为高斯-约旦消去法)是一种利用矩阵消元的方法，将线性方程组化为改进的阶梯形式，从而解出线性方程组的解的方法。
具体而言，高斯消去法将线性方程组的系数矩阵化为一个上三角矩阵，再利用回带法求解线性方程组的解。
三、实验内容1.1、用高斯消去法解线性方程组在具体实验中，我们将使用高斯消去法来解决下述的线性方程组。
5x+2y+z=102x+6y+2z=14x-y+10z=25为了使用高斯消去法来解这个方程组，首先需要将系数矩阵A进行变换，消除A矩阵中第一列中的下角元素，如下所示:1, 2/5, 1/50, 28/5, 18/50, 0, 49/28接着使用回代法来计算该方程组的解。
回代法的过程是从下往上进行的，具体步骤如下：第三个方程的解：z=49/28；第二个方程的解： y=（14-2z-2x）/6；第一个方程的解： x=（10-2y-z）/5。
1.2、分析误差在使用高斯消去法求解线性方程组时，一般会出现截断误差，导致得到的解与真实解之间存在一些误差。
截断误差的大小和矩阵的维数有关。
为了估计截断误差，我们使用矩阵B来生成误差，在具体实验中，我们将使用下面的矩阵：我们来计算该矩阵的行列式，如果方程组有唯一解，则行列性不为0。
本例中，行列式的值是 -1，因此方程组有唯一解。
然后我们计算真实解和高斯消去法得到的解之间的误差，具体公式如下所示：误差 = 真实解的范数 - 高斯消去法得到的解的范数其中，范数的定义如下：||x||1=max{|xi|}; ||x||2=sqrt{(|x1|^2 + |x2|^2 + ... + |xn|^2)}四、实验步骤1、将高斯消去法的每一个步骤翻译成代码，并保存为一个独立的函数。
2、将代码上传至 Python 交互式环境，并使用高斯消去法来解线性方程组。
北航数值分析大作业第一题
1 算法方案 1.1 λ1，λ501，λs 的计算
(1) (2) (3) (4) (5) 将矩阵 A[501][501]以压缩存储后的形式 C[5][501]输入 使用一次幂法得到按模最大的特征值 矩阵向左平移 λm 距离（A-λmI） ，再使用一次幂法得到按模最大的特征值 s，则 λm1=s-λm1 比较 λm1 和 λm2 的大小与正负，得到 λ 和 λ501 对 A 使用一次反幂法得到按模最小的特征值 λs
while (e>=pow(10,-12)); return 1/be;//返回 1/be2 作为矩阵 m[5][501]的按模最小向量 } //333333333333333333333333333333333333333333333333333333333333333333333333 33333333333333333333333333333333333333333333333333333333333333333333333 double det(double c[1+r+s][q]) { int max3(int a,int b,int c); int fmax2(int a,int b); int fmin2(int a,int b); int i,j,k,t; double sum,det=1; for(k=1;k<=q;k++) { for(j=k;j<=fmin2(k+s,q);j++)//求 ukj { sum=0; for(t=max3(1,k-r,j-s);t<=k-1;t++) { sum=sum+c[k-t+s][t-1]*c[t-j+s][j-1]; } c[k-j+s][j-1]=c[k-j+s][j-1]-sum; }
北航数值分析第一次大作业
b2[i-1]=b[i-1]-sum3; } x[n-1]=b2[n-1]/C[s][n-1]; for(i=n-1;i>=1;i--) { double sum4=0; for(int t=i+1;t<=min(i+s,n);t++) { sum4+=C[i-t+s][t-1]*x[t-1]; } x[i-1]=(b2[i-1]-sum4)/C[s][i-1]; } } /*反幂法*/ double FMF(double C[m][n]) { LU(C); for(int k=1;k<=n;k++) u[k-1]=1; /*为迭代初始向量赋值*/ beta1=beta2=0; do { ent=0; for(int i=1;i<=n;i++) ent+=u[i-1]*u[i-1]; ent=sqrt(ent); for(i=1;i<=n;i++) y[i-1]=u[i-1]/ent; HD(C,y,u); beta1=beta2; beta2=0; for(i=1;i<=n;i++) { beta2+=y[i-1]*u[i-1]; } }while(fabs(1/beta2-1/beta1)/fabs(1/beta2)>1.0e-12); return 1/beta2; } /*求 detA*/ double det(double C[m][n]) { LU(C); double detA=1; for(int j=1;j<=n;j++)
数值分析第一次作业
姓名：吴少波 学号：SY1105513
一、算法的设计方案 1.将带状矩阵 A 压缩为矩阵 C 存储。先用幂法算出 A 按模最大的特征值，记为 maxLambda， 再 将 其 平 移 ，用 带 原点 平 移 的 幂 法求 A-maxLambdaI 按模 最 大的 特 征 值 ， 记为 p1 ， 记 p2=p1+maxLambda，比较 maxLambda 和 p2 的大小，大的为λ 501，小的为λ 1。 用反幂法求解λ s 时，其中需解方程 Auk=yk-1，先把矩阵 A LU 分解（不列主元） ，再在每次循环 迭代时回代求解。 2.将 A 平移μ k（k=1,2，…，39）个单位，用带原点平移的反幂法求与μ k（k=1,2，…，39） 最接近的 39 个特征值。 3.cond(A)2=│maxLambda / λ s│ A 的行列式等于把 A LU 分解后 A 所有对角线上元素的乘积。 二、源程序（VC6.0 环境下的 C 语言） #include<stdio.h> #include<stdlib.h> #include<math.h> #include<malloc.h> #define m 5 #define n 501 #define r 2 #define s 2 double C[m][n]; double u[n]; double y[n]; double ent,beta1,beta2; void YS(); /*将带状矩阵 A 压缩为 C*/ int max(int a,int b); /*两数求较大的一个*/ int min(int a,int b); /*两数求较小的一个*/ double MF(double C[m][n]); /*幂法*/ double FMF(double C[m][n]); /*反幂法*/ void LU(double C[m][n]); /*LU 分解*/ void HD(double C[m][n],double b[n],double x[n]); /*回代过程*/ double det(double C[m][n]); /*求 detA*/ double Move_MF(double C[m][n],double maxLambda); /*带原点平移的幂法*/ double Move_FMF(double C[m][n],double p); /*带原点平移的反幂法*/ /**主函数**/ void main() { /*定义变量*/ double maxLambda=0,minLambda=0,condA,detA,Lambda1,Lambda501,p1,p2,Mu_k,Lambdaik; /*算第一题*/
北航研究生数值分析编程大作业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、最大特征值λ501，最小特征值λ1，按模最小特征值λs的计算方法首先使用一次幂法运算可以得到矩阵的按模最大的特征值λ，λ必为矩阵A的最大或最小特征值，先不做判断。
对原矩阵A进行一次移项，即（A-λI），在进行一次幂法运算，可以得到另一个按模最大特征值λ0。
比较λ和λ的大小，较大的即为λ501，较小的即为λ1。
对矩阵A进行一次反幂法运算，即可得到按模最小特征值λs。
2、A与μk 值最接近的特征值λik的计算方法首先计算出k所对应的μk 值，对原矩阵A进行一次移项，即（A-μkI），得到一个新的矩阵，对新矩阵进行一次反幂法运算，即可得到一个按模最小特征值λi 。
则原矩阵A与μk值最接近的特征值λik=λi+μk。
3、A的（谱范数）条件数cond（A）2的计算方法其中错误!未找到引用源。
矩阵A的按模最大和按模最小特征值。
（二）程序编写思路。
由于算法要求A的零元素不存储，矩阵A本身为带状矩阵，所以本题的赋值，LU分解，反幂法运算过程中，均应采用Doolittle分解法求解带状线性方程组的算法思路。
幂法、反幂法和LU分解均是多次使用，应编写子程序进行反复调用。
二、源程序：#include<stdio.h>#include<iostream>#include<stdlib.h>#include<math.h>#include<float.h>#include<iomanip> /*头文件*//*定义全局变量*/#define N 502 /*取N为502，可实现从1到501的存储，省去角标变换的麻烦*/ #define epsilon 1.0e-12 /*定义精度*/#define r 2 /*r，s为带状矩阵的半带宽，本题所给矩阵二者都是2*/ #define s 2double c[6][N]; /*定义矩阵存储压缩后的带状矩阵*/double fuzhi(); /*赋值函数*/void LUfenjie(); /*LU分解程序*/int max(int a,int b); /*求两个数字中较大值*/int min(int a,int b); /*求两个数字中较小值*/double mifa(); /*幂法计算程序*/double fanmifa(); /*反幂法计算程序*/double fuzhi() /*赋值程序，按行赋值，行从1到5，列从1到501，存储空间的第一行第一列不使用，角标可以与矩阵一一对应，方便书写程序*/{int i,j;i=1;for(j=3;j<N;j++){c[i][j]=-0.064;}i=2;for(j=2;j<N;j++){c[i][j]=0.16;}i=3;for(j=1;j<N;j++){c[i][j]=(1.64-0.024*j)*sin(0.2*j)-0.64*exp(0.1/j);}i=4;for(j=1;j<N-1;j++){c[i][j]=0.16;}i=5;for(j=1;j<N-2;j++){c[i][j]=-0.064;}return(c[i][j]);}int max(int a,int b){ return((a>b)?a:b);}int min(int a,int b){ return((a<b)?a:b);}void LUfenjie() /*LU分解程序，采用的是带状矩阵压缩存储后的LU分解法*/{double temp;int i,j,k,t;for(k=1;k<N;k++){ for(j=k;j<=min(k+s,N-1);j++){temp=0;for(t=max(1,max(k-r,j-s));t<=(k-1);t++){temp=temp+c[k-t+s+1][t]*c[t-j+s+1][j];}c[k-j+s+1][j]=c[k-j+s+1][j]-temp;}for(i=k+1;i<=min(k+r,N-1);i++){temp=0;for(t=max(1,max(i-r,k-s));t<=(k-1);t++){temp=temp+c[i-t+s+1][t]*c[t-k+s+1][k];}c[i-k+s+1][k]=(c[i-k+s+1][k]-temp)/c[s+1][k];}}}double mifa() /*幂法计算程序*/ {double u0[N],u1[N];double temp,Lu,beta=0,beta0;int i,j;for(i=1;i<N;i++) /*选取迭代初始向量*/{u0[i]=1;}do{beta0=beta;temp=0;for(i=1;i<N;i++){temp=temp+u0[i]*u0[i]; }Lu=sqrt(temp);for(i=1;i<N;i++){u1[i]=u0[i]/Lu;}for(i=1;i<N;i++){temp=0;for(j=max(i-1,1);j<=min(i+2,N-1);j++){temp=temp+c[i-j+s+1][j]*u1[j]; }u0[i]=temp;} //新的u0temp=0;for(i=1;i<N;i++){temp=temp+u1[i]*u0[i]; }beta=temp;}while(fabs(beta-beta0)/fabs(beta)>=epsilon); /*迭代运行条件判断*/return(beta);}double fanmifa() /*反幂法计算程序*/{double u0[N],u1[N],u2[N];double temp,Lu,beta=0,beta0;int i,t;for(i=1;i<N;i++) /*选取迭代初始向量*/{u0[i]=1;}do{beta0=beta;temp=0;for(i=1;i<N;i++){temp=temp+u0[i]*u0[i]; }Lu=sqrt(temp);for(i=1;i<N;i++){u1[i]=u0[i]/Lu;u2[i]=u1[i];}fuzhi();LUfenjie();/*带状矩阵压缩存储并进行LU分解后，求解线性方程组得到迭代向量u k，即程序中的u0*/for(i=2;i<N;i++){ temp=0;for(t=max(1,i-r);t<=(i-1);t++){temp=temp+c[i-t+s+1][t]*u2[t];}u2[i]=u2[i]-temp;}u0[N-1]=u2[N-1]/c[s+1][N-1];for(i=N-2;i>=1;i--){ temp=0;for(t=i+1;t<=min(i+s,N-1);t++){temp=temp+c[i-t+s+1][t]*u0[t];}u0[i]=(u2[i]-temp)/c[s+1][i];}temp=0;for(i=1;i<N;i++){temp=temp+u1[i]*u0[i]; }beta=temp;beta=1/beta; /*beta即为所求特征值，可直接返回*/}while(fabs(beta-beta0)/fabs(beta)>=epsilon); /*迭代运行条件判断*/return(beta);}void main(){double u[40]; /*定义数组，存放k值运算得到的μk值*/double lambda1,lambda501,lambdak,a,b,d,cond,det;int i,j,k;fuzhi();a=mifa(); /*幂法计算按模最大值*/fuzhi();d=fanmifa(); /*反幂法计算按模最小值*/fuzhi();for(j=1;j<N;j++){c[3][j]=c[3][j]-a;}b=mifa()+a; /*移项后幂法计算按模最大值*/if(a>b) /*比较两个按模最大值大小，并相应输出最大特征值λ501和最小特征值λ1*/ {lambda1=b;lambda501=a;printf("矩阵A最小的特征值lambda1=%13.11e\n",lambda1);printf("矩阵A最大的特征值lambda501=%13.11e\n",lambda501);}else{lambda1=a;lambda501=b;printf("矩阵A最小的特征值lambda1=%13.11e\n",lambda1);printf("矩阵A最大的特征值lambda501=%13.11e\n",lambda501);}printf("矩阵A按模最小特征值lambdas=%13.11e\n",d); /*输出按模最小特征值λs*/for(k=1;k<40;k++) /*对每一个进行移项反幂法运算，求出最接近μk的特征值并输出*/ {u[k]=(lambda501-lambda1)*k/40+lambda1;fuzhi();for(j=1;j<N;j++){c[3][j]=c[3][j]-u[k];}lambdak=fanmifa()+u[k];i=k;printf("矩阵A最接近uk特征值lambdak%d=%13.11e\n",i,lambdak);}cond=fabs(a/d);printf("A的条件数=%13.11e\n",cond); /*计算A条件数并输出*/fuzhi(); /*计算A的行列式值并输出*/LUfenjie();det=1;for(i=1;i<N;i++){det=det*c[3][i];}printf("行列式的值detA=%13.11e\n",det);}三、程序的运行结果：四、初始向量的选取对计算结果的影响：（一）选取形式不变，数值变换1、取u0为[0.5,0.5………..0.5]，运行结果如下：2、取u0为[50,50………..50]，运行结果如下：从运行结果来看，此类初始向量的选取对结果不会产生影响，即使选成0，结果也不变化。
北航数值分析报告第一次大作业(幂法反幂法)
一、问题分析与算法描述1. 问题的提出：〔1〕用幂法、反幂法求矩阵的按摸最大和最小特征值，并求出相应的特征向量。
其中要求：迭代精度达到。
〔2〕用带双步位移的QR法求上述的全部特征值，并求出每一个实特征值相应的特征向量。
2. 算法的描述：(1) 幂法幂法主要用于计算矩阵的按摸为最大的特征值和相应的特征向量。
其迭代格式为：终止迭代的控制选用。
幂法的使用条件为实矩阵A具有n个线性无关的特征向量，其相应的特征值满足不等式或幂法收敛速度与比值或有关，比值越小，收敛速度越快。
(2) 反幂法反幂法用于计算实矩阵A按摸最小的特征值，其迭代格式为：每迭代一次都要求解一次线性方程组。
当k足够大时，，可近似的作为矩阵A的属于的特征向量。
比值越小，收敛的越快。
反幂法要求矩阵A非奇异。
(3) 带双步位移的QR分解法QR方法适用于计算一般实矩阵的全部特征值，尤其适用于计算中小型实矩阵的全部特征值。
本算例中采用带双步位移的QR方法，可加速收敛，其迭代格式为：二、计算结果与分析1. 计算结果：(1) 幂法：初始条件：最大迭代次数L=1000；向量计算结果：第1次迭代结果：最大特征值：0.00000e+000第2次迭代结果：最大特征值：2.48910e+000 相对误差：1.00000e+000 第3次迭代结果：最大特征值：1.67719e+000 相对误差：第4次迭代结果：最大特征值：-2.10960e+000 相对误差：1.79503e+000 第5次迭代结果：最大特征值：-6.13203e-001 相对误差：2.44030e+000 ……第794次迭代结果：最大特征值：-1.97638e+000 相对误差：最大特征值：-1.97638e+000 相对误差：********************最终迭代结果***************特征值：-1.97638e+000 相对误差：迭代次数：795(2) 反幂法：初始条件：最大迭代次数L=1000；向量运行结果：第1次迭代结果：最大特征值：1.07542e+000第2次迭代结果：最大特征值：-3.66550e+000 相对误差：1.29339e+000 第3次迭代结果：最大特征值：1.22709e+001 相对误差：1.29871e+000 第4次迭代结果：最大特征值：-1.03421e+000 相对误差：1.28650e+001 第5次迭代结果：最大特征值：相对误差：……第995次迭代结果：最大特征值：相对误差：第996次迭代结果：最大特征值：相对误差：最大特征值：相对误差：第998次迭代结果：最大特征值：相对误差：第999次迭代结果：最大特征值：相对误差：第1000次迭代结果：最大特征值：相对误差：******************************超过最大设定迭代次数，迭代失败！(3) 带双步位移的QR法：初始条件：最大迭代次数L=1000；向量运行结果：全部特征值：特征向量〔经谱X数归一化〕：实特征值对应特征向量：-0.062705 -0.022368 0.304372 0.064466 0.521833 -0.157024 0.136942 -0.218108 0.250264 -0.043064 -0.228688 -0.184632 -0.072871 0.124721 0.029070 0.102566 -0.136358 0.167727 0.085747 0.546165 实特征值对应特征向量：-0.018001 0.019652 0.273447 0.070528 0.274896 -0.144015 0.048385 0.376439 -0.583051 -0.054008 -0.168682 -0.113430 -0.034709 0.009204 0.472291 0.125664 -0.190617 0.113145 0.046278 0.059871 实特征值对应特征向量：0.106861 0.087709 -0.024967 -0.020897 0.064302 0.034047 0.535143 0.046383 0.028832 0.003479-0.097276 -0.383801 0.089445 -0.039560 -0.036928 -0.021330 0.014811 0.705836 -0.108904 0.082022 实特征值对应特征向量：-0.055201 0.003399 0.242191 0.102847 0.372470 -0.372826 0.113953 0.240659 -0.310401 -0.076590 -0.244632 -0.192549 -0.077259 0.263328 0.201662 0.154166 -0.407814 0.186782 0.094649 0.173302 实特征值对应特征向量：0.427828 -0.546801 0.007822 -0.382580 0.025199 0.012788 0.033241 0.005389 -0.004065 0.043524 -0.032112 -0.044233 0.135395 -0.006564 0.001214 0.020165 0.011678 0.050001 -0.585765 0.013115 实特征值对应特征向量：0.236032 -0.139250 -0.008143 0.638527 -0.009049 -0.002911 -0.001307 0.003054 0.006515 -0.030134 0.012712 0.011368 -0.018792 -0.001753 -0.005749 -0.014290 -0.005292 -0.014591 0.717590 0.001369 实特征值对应特征向量：-0.227404 -0.048154 0.022615 0.297305 0.070372 0.039927 0.078503 0.015822 -0.012182 0.605334 -0.083616 -0.106270 -0.573963 -0.019907 0.003839 0.051362 0.036567 0.115613 0.332707 0.036954 实特征值对应特征向量：-0.027768 -0.051081 -0.159642 -0.054573 -0.084441 0.118378 0.029553 0.211088 0.203867 0.0486272. 结果分析以上三种方法中，幂法计算共进展了795次迭代才达到收敛，计算量较大，收敛性不好；反幂法计算结果未能收敛，通过进一步分析发现，这是因为反幂法迭代程序未考虑按模最小特征值为复数的情况，造成迭代失败。
北航硕士研究生数值分析大作业一
数值分析—计算实习作业一学院：17系专业：精密仪器及机械姓名：张大军学号：DY14171142014-11-11数值分析计算实现第一题报告一、算法方案算法方案如图1所示。
（此算法设计实现完全由本人独立完成）图1算法方案流程图二、全部源程序全部源程序如下所示#include <iostream.h>#include <iomanip.h>#include <math.h>int main(){double a[501];double vv[5][501];double d=0;double r[3];double uu;int i,k;double mifayunsuan(double *a,double weiyi);double fanmifayunsuan(double *a,double weiyi);void yasuo(double *A,double (*C)[501]);void LUfenjie(double (*C)[501]);//赋值语句for(i=1;i<=501;i++){a[i-1]=(1.64-0.024*i)*sin(0.2*i)-0.64*exp(0.1/i);}//程序一：使用幂方法求绝对值最大的特征值r[0]=mifayunsuan(a,d);//程序二：使用幂方法求求平移λ[0]后绝对值最大的λ，得到原矩阵中与最大特征值相距最远的特征值d=r[0];r[1]=mifayunsuan(a,d);//比较λ与λ-λ[0]的大小，由已知得if(r[0]>r[1]){d=r[0];r[0]=r[1];r[1]=d;}//程序三：使用反幂法求λr[2]=fanmifayunsuan(a,0);cout<<setiosflags(ios::right);cout<<"λ["<<1<<"]="<<setiosflags(ios::scientific)<<setprecision(12)<<r[0]<<endl;cout<<"λ["<<501<<"]="<<setiosflags(ios::scientific)<<setprecision(12)<<r[1]<<endl;cout<<"λ[s]="<<setiosflags(ios::scientific)<<setprecision(12)<<r[2]<<endl;//程序四：求A的与数u最接近的特征值for(k=1;k<40;k++){uu=r[0]+k*(r[1]-r[0])/40;cout<<"最接近u["<<k<<"]"<<"的特征值为"<<setiosflags(ios::scientific)<<setprecision(12)<<fanmifayunsuan(a,uu)<<endl;}//程序五：谱范数的条件数是绝对值最大的特征值除以绝对值最小的特征值的绝对值cout<<"cond(A)2="<<fabs(r[0]/r[2])<<endl;//程序六：A的行列式的值就是A分解成LU之U的对角线的乘积yasuo(a,vv);LUfenjie(vv);uu=1;for(i=0;i<501;i++){uu=uu*vv[2][i];}cout<<"Det(A)="<<uu<<endl;return 1;}double mifayunsuan(double *a,double weiyi){int i,k;double b=0.16;double c=-0.064;double ee,w,v1,v2,mm,sum;double u[501];double y[505]={0};for(i=0;i<501;i++)u[i]=1;//给u赋初值if (weiyi!=0){for (i=0;i<501;i++)a[i]-=weiyi;}ee=1;k=0;//使得初始计算时进入循环语句while(ee>1e-12){mm=0;for(i=0;i<501;i++){mm=mm+u[i]*u[i];}w=sqrt(mm);for(i=0;i<501;i++){y[i+2]=u[i]/w;//注意此处编程与书上不同，之后会解释它的巧妙之处1 }for(i=0;i<501;i++){u[i]=c*y[i]+b*y[i+1]+a[i]*y[i+2]+b*y[i+3]+c*y[i+4];//1显然巧妙之处凸显出来}sum=0;for(i=0;i<501;i++){sum+=y[i+2]*u[i];}v1=v2;v2=sum;//去除特殊情况，减少漏洞if(k==0){k++;}else{ee=fabs(v2-v1)/fabs(v2);}}if (weiyi!=0){for (i=0;i<501;i++)a[i]+=weiyi;}//还原A矩阵return (v2+weiyi);}double fanmifayunsuan(double *a,double weiyi){int i,k;double b=0.16;double c=-0.064;double ee,w,v1,v2,mm,sum;double u[501];double y[501];double C[5][501];void yasuo(double *A,double (*C)[501]);void LUfenjie(double (*C)[501]);void qiuU(double (*C)[501],double *y,double *u);//把A阵压缩到C阵中for(i=0;i<501;i++)u[i]=1;//给u赋初值if (weiyi!=0){for (i=0;i<501;i++)a[i]-=weiyi;}yasuo(a,C);LUfenjie(C);ee=1;k=0; //使得初始计算时进入循环语句while(ee>1e-12){mm=0;for(i=0;i<501;i++){mm=mm+u[i]*u[i];}w=sqrt(mm);for(i=0;i<501;i++){y[i]=u[i]/w;}qiuU(C,y,u);sum=0;for(i=0;i<501;i++){sum+=y[i]*u[i];}v1=v2;v2=sum;//去除特殊情况，减少漏洞if(k==0){k++;}else{ee=fabs(1/v2-1/v1)/fabs(1/v2);}}if (weiyi!=0){for (i=0;i<501;i++)a[i]+=weiyi;}//还原A矩阵return (1/v2+weiyi);}void yasuo(double *A,double (*C)[501]){double b=0.16;double c=-0.064;int i;for(i=0;i<501;i++){C[0][i]=c;C[1][i]=b;C[2][i]=A[i];C[3][i]=b;C[4][i]=c;}}void LUfenjie(double (*C)[501]){int k,t,j;int r=2,s=2;double sum;int minn(int ,int );int maxx(int ,int );for(k=0;k<501;k++){for(j=k;j<=minn(k+s,501-1);j++){if(k==0)sum=0;else{sum=0;for(t=maxx(k-r,j-s);t<k;t++){sum=sum+C[k-t+s][t]*C[t-j+s][j];}}C[k-j+s][j]=C[k-j+s][j]-sum;}for(j=k+1;j<=minn(k+r,501-1);j++){if(k<501-1){if(k==0)sum=0;else{sum=0;for(t=maxx(j-r,k-s);t<k;t++){sum=sum+C[j-t+s][t]*C[t-k+s][k];}}C[j-k+s][k]=(C[j-k+s][k]-sum)/C[s][k];}}}}void qiuU(double (*C)[501],double *y,double *u){int i,t;double b[501];double sum;int r=2,s=2;int minn(int ,int );int maxx(int ,int );for(i=0;i<501;i++){b[i]=y[i];}for(i=1;i<501;i++){sum=0;for(t=maxx(0,i-r);t<i;t++){sum=sum+C[i-t+s][t]*b[t];}b[i]=b[i]-sum;}u[500]=b[500]/C[s][500];for(i=501-2;i>=0;i--){sum=0;for(t=i+1;t<=minn(i+s,500);t++){sum=sum+C[i-t+s][t]*u[t];}u[i]=(b[i]-sum)/C[s][i];}}int minn(int x,int y){int min;if(x>y)min=y;elsemin=x;return min;}int maxx(int b,int c){int max;if(b>c){if(b>0)max=b;elsemax=0;}else{if(c>0)max=c;elsemax=0;}return max;}三、特征值以及的值λ[1]=-1.070011361502e+001 λ[501]=9.724634098777e+000λ[s]=-5.557910794230e-003最接近u[1]的特征值为-1.018293403315e+001最接近u[2]的特征值为-9.585707425068e+000最接近u[3]的特征值为-9.172672423928e+000最接近u[4]的特征值为-8.652284007898e+000最接近u[5]的特征值为-8.0934********e+000最接近u[6]的特征值为-7.659405407692e+000最接近u[7]的特征值为-7.119684648691e+000最接近u[8]的特征值为-6.611764339397e+000最接近u[9]的特征值为-6.0661********e+000最接近u[10]的特征值为-5.585101052628e+000最接近u[11]的特征值为-5.114083529812e+000最接近u[12]的特征值为-4.578872176865e+000最接近u[13]的特征值为-4.096470926260e+000最接近u[14]的特征值为-3.554211215751e+000最接近u[15]的特征值为-3.0410********e+000最接近u[16]的特征值为-2.533970311130e+000最接近u[17]的特征值为-2.003230769563e+000最接近u[18]的特征值为-1.503557611227e+000最接近u[19]的特征值为-9.935586060075e-001最接近u[20]的特征值为-4.870426738850e-001最接近u[21]的特征值为2.231736249575e-002最接近u[22]的特征值为5.324174742069e-001最接近u[23]的特征值为1.052898962693e+000最接近u[24]的特征值为1.589445881881e+000最接近u[25]的特征值为2.060330460274e+000最接近u[26]的特征值为2.558075597073e+000最接近u[27]的特征值为3.080240509307e+000最接近u[28]的特征值为3.613620867692e+000最接近u[29]的特征值为4.0913********e+000最接近u[30]的特征值为4.603035378279e+000最接近u[31]的特征值为5.132924283898e+000最接近u[32]的特征值为5.594906348083e+000最接近u[33]的特征值为6.080933857027e+000最接近u[34]的特征值为6.680354092112e+000最接近u[35]的特征值为7.293877448127e+000最接近u[36]的特征值为7.717111714236e+000最接近u[37]的特征值为8.225220014050e+000最接近u[38]的特征值为8.648666065193e+000最接近u[39]的特征值为9.254200344575e+000cond(A)2=1.925204273902e+003 Det(A)=2.772786141752e+118四、现象讨论在大作业的程序设计过程当中，初始向量的赋值我顺其自然的设为第一个分量为1，其它分量为0的向量，计算结果与参考答案存在很大差别，计算结果对比如下图2所示（左侧为正确结果，右侧为错误结果），导致了我花了很多的时间去检查程序算法。
北航数值分析大作业一
北航数值分析大作业一————————————————————————————————作者: ————————————————————————————————日期：ﻩ《数值分析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);/／测试不同迭代初始向量对λ１计算结果的影响。
北航数值分析第一次大作业
一、算法的设计方案1、求矩阵最大特征值，最小特征值与按模最小特征值的方法首先用幂法求出矩阵A 的一个特征值λ，则其必为最大特征值与最小特征值二者其一，之后对矩阵A 进行一次移项，即A-λI ，然后再次用幂法求出另一个按模最大特征值，再比较这两个值的大小，则较大的为矩阵A 的最大特征值，较小的为矩阵A 的最小特征值。
用反幂法可以求得矩阵的按模最小特征值λs 2、求矩阵A 与k μ最接近的特征值k i λ可以先对矩阵A 进行移项，即A-k μI ，对这个移项后的矩阵用反幂法求出按模最小的特征值，然后再加上k μ，就求出所要求的k i λ。
3、求矩阵A 的条件数cond(A)2和行列式detA由于矩阵A 是非奇异的实对称矩阵，所以可以用以下公式方便地求出矩阵A 的条件数cond(A)2=sλλ501对于矩阵A 行列式的求法也比较简单。
由于在用反幂法的过程中对A 进行了Doolittle LU 分解，所以detA=detL*detU ，而detL=1，detU 可以用对角线元素相乘方便地算出，所以detA 就是U 阵对角线元素的乘积。
4、几点说明由于A 中的零元素都不存储，所以在存储矩阵的时候采用书上26页的压缩存储方式。
在反幂法中采用LU 分解求解带状线性方程组的算法来求解每一次迭代的方程组，由于每一次方程左边的系数都相同，所以只要进行一次LU 分解即可。
因为幂法，反幂法，LU 分解，求最大值与最小值在程序编写的过程中多次用到，所以这几项作为子函数单独进行编写。
二、源程序如下：#include "stdio.h"#include "math.h"# define s 2# define r 2# define N 501double c[5][N]={0};double lameda[40];double max(double x,double y);double min(double x,double y);double mifa();double fanmifa();void LUfenjie();void main(){int i=0,j=0;/*============对数组进行赋值==============*/ for(j=3;j<=N;j++)c[0][j-1]=-0.064;for(j=2;j<=N;j++)c[1][j-1]=0.16;for(j=1;j<=N;j++)c[2][j-1]=(1.64-0.024*j)*sin(0.2*j)-0.64*exp(0.1/j);for(j=1;j<=N-1;j++)c[3][j-1]=0.16;for(j=1;j<=N-2;j++)c[4][j-1]=-0.064;/*========幂法求最大和最小特征值==============*/ double a=mifa();for(j=1;j<=N;j++)c[2][j-1]-=a;double b=mifa()+a;double lameda501=max(a,b);double lameda1=min(a,b);printf("矩阵A最大的特征值=%13.11e\n",lameda501);printf("矩阵A最小的特征值=%13.11e\n",lameda1);/*========反幂法求绝对值最小特征值===========*/for(j=1;j<=N;j++)c[2][j-1]+=a;double lamedas=fanmifa();printf("矩阵A按模最小的特征值=%13.11e\n",lamedas);/*========求条件数和行列式的值===========*/double detA=1;for(j=1;j<=N;j++)detA*=c[2][j-1];printf("矩阵A的行列式=%13.11e\n",detA);double condA=fabs(lameda501/lamedas);printf("矩阵A的条件数=%13.11e\n",condA);/*========反幂法求与uk最接近的特征值========*/for(int k=1;k<40;k++){ for(j=3;j<=N;j++)c[0][j-1]=-0.064;for(j=2;j<=N;j++)c[1][j-1]=0.16;for(j=1;j<=N;j++)c[2][j-1]=(1.64-0.024*j)*sin(0.2*j)-0.64*exp(0.1/j);for(j=1;j<=N-1;j++)c[3][j-1]=0.16;for(j=1;j<=N-2;j++)c[4][j-1]=-0.064;for(j=1;j<=N;j++)c[2][j-1]-=(lameda1+k*(lameda501-lameda1)/40);lameda[k]=fanmifa()+(lameda1+k*(lameda501-lameda1)/40);printf("矩阵A最接近u%d的特征值=%13.11e\n",k,lameda[k]);}}double max(double x,double y) //求两数中的最大值{double z;z=x>y ? x:y;return(z);}double min(double x,double y) //求两数中的最小值{double z;z=x<y ? x:y;return(z);}double mifa() //幂法求按模最大特征值{double u[N]={0};double sum=0;double zero=0;double y[N]={0};double b1=0,b2=0;int i=0,j=0;for(i=0;i<N;i++)u[i]=1;do{ sum=0;for(i=0;i<N;i++)sum+=u[i]*u[i];for(i=0;i<N;i++)y[i]=u[i]/sqrt(sum);for(i=0;i<N;i++){ u[i]=0;for(j=max(0,i-2);j<=min(i+2,N-1);j++){u[i]+=c[i-j+s][j]*y[j];}}b2=b1;zero=0;for(i=0;i<N;i++)zero+=y[i]*u[i];b1=zero;}while((fabs(b1-b2)/fabs(b1))>1e-12);return(b1);}double fanmifa() //反幂法求按模最小特征值{double u[N]={0};double sum=0;double zero=0;double y[N]={0};double b[N]={0};double b1=0,b2=0;int i=0,j=0,t=0;for(i=0;i<N;i++)u[i]=1;LUfenjie();do{ sum=0;for(i=0;i<N;i++)sum+=u[i]*u[i];for(i=0;i<N;i++){b[i]=u[i]/sqrt(sum);y[i]=b[i];}for(i=1;i<N;i++)for(t=max(0,i-r);t<i;t++)y[i]-=c[i-t+s][t]*y[t];u[N-1]=y[N-1]/c[s][N-1];for(i=N-2;i>=0;i--){ for(t=i+1;t<=min(i+s,N-1);t++)y[i]-=c[i-t+s][t]*u[t];u[i]=y[i]/c[s][i];}b2=b1;zero=0;for(i=0;i<N;i++)zero+=b[i]*u[i];b1=1/zero;}while((fabs(b1-b2)/fabs(b1))>1e-12);return(b1);}void LUfenjie() //对矩阵做LU分解{ int k=0;int j=0;int i=0;int t=0;for(k=0;k<N;k++){ for(j=k;j<=min(k+s,N-1);j++)for(t=max(max(0,k-r),j-s);t<k;t++)c[k-j+s][j]-=c[k-t+s][t]*c[t-j+s][j];for(i=k+1;i<=min(k+r,N-1);i++){ for(t=max(max(0,i-r),k-s);t<k;t++)c[i-k+s][k]-=c[i-t+s][t]*c[t-k+s][k];c[i-k+s][k]/=c[s][k];}}}三、运行结果如下四、初始向量对计算结果的影响在本程序的编写中，为了方便起见，所以迭代向量的初值选为u=[1 1 1 ...1]。
北航数值分析计算实习报告一
北京航空航天大学《数值分析》计算实习报告第一大题学 院：自动化科学与电气工程学院 专 业： 控制科学与工程 学 生 姓 名： 学 号： 教 师： 电 话： 完 成 日 期： 2015年11月6日北京航空航天大学Beijing University of Aeronautics and Astronautics实习题目：第一题 设有501501⨯的实对称矩阵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 λ，并且有1.求1λ，501λ和s λ的值。
2.求A 的与数4015011λλλμ-+=kk 最接近的特征值)39,,2,1(⋅⋅⋅=k k i λ。
3.求A 的(谱范数)条件数2)A (cond 和行列式detA 。
说明：1.在所用的算法中，凡是要给出精度水平ε的，都取12-10=ε。
2.选择算法时，应使矩阵A 的所有零元素都不储存。
3.打印以下内容：（1）全部源程序；（2）特征值），，39,...,2,1(,s 5011=k k i λλλλ以及A det ,)A (cond 2的值。
4.采用e 型输出实型数，并且至少显示12位有效数字。
一、算法设计方案 1、求1λ，501λ和s λ的值。
由于||min ||,501150121i i s λλλλλ≤≤=≤⋅⋅⋅≤≤，可知绝对值最大特征值必为1λ和501λ其中之一，故可用幂法求出绝对值最大的特征值λ，如果λ=0，则1λ=λ，否则501λ=λ。
将矩阵A 进行一下平移：I -A A'λ= （1）对'A 用幂法求出其绝对值最大的特征值'λ，则A 的另一端点特征值1λ或501λ为'λ+λ。
s λ为按模最小特征值，||min ||5011i i s λλ≤≤=，可对A 使用反幂法求得。
北航数值分析第一次大作业(幂法反幂法)
一、问题分析及算法描述1. 问题的提出：（1）用幂法、反幂法求矩阵A =[a ij ]20×20的按摸最大和最小特征值，并求出相应的特征向量。
其中 a ij ={sin (0.5i +0.2j ) i ≠j 1.5cos (i +1.2j ) i =j要求：迭代精度达到10−12。
（2）用带双步位移的QR 法求上述的全部特征值，并求出每一个实特征值相应的特征向量。
2. 算法的描述：(1) 幂法幂法主要用于计算矩阵的按摸为最大的特征值和相应的特征向量。
其迭代格式为：{ 任取非零向量u 0=(h 1(0)，⋯，h n (0))T|h r (k−1)|=max 1≤j≤n |h r (k−1)| y ⃑ k−1=u ⃑ k−1|h r (k−1)| u ⃑ k =Ay ⃑ k−1=(h 1(k )，⋯，h n (k ))T βk =sgn (h r (k−1))h r (k ) (k =1，2，⋯) 终止迭代的控制选用≤ε。
幂法的使用条件为n ×n 实矩阵A 具有n 个线性无关的特征向量x 1，x 2，⋯，x n ，其相应的特征值λ1，λ2，⋯，λn 满足不等式|λ1|>|λ2|≥|λ3|≥⋯≥|λn |或λ1=λ2=⋯=λm|λ1|>|λm+1|≥|λm+2|≥⋯≥|λn |幂法收敛速度与比值|λ2λ1|或|λm+1λ1|有关，比值越小，收敛速度越快。
(2) 反幂法反幂法用于计算n ×n 实矩阵A 按摸最小的特征值，其迭代格式为：{任取非零向量u 0∈R nηk−1=√u ⃑ k−1T u ⃑ k−1 y ⃑ k−1=u ⃑ k−1ηk−1⁄ Au ⃑ k =y ⃑ k−1 βk =y ⃑ k−1u ⃑ k (k =1，2，⋯) 每迭代一次都要求解一次线性方程组Au ⃑ k =y ⃑ k−1。
当k 足够大时，λn ≈1βk ，y ⃑ k−1可近似的作为矩阵A 的属于λn 的特征向量。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
一、问题分析及算法描述编写程序，分别用列主元的Gauss 消去法和LU 分解法求解下面线型代数方程组AX=b 的解，其中A 为N ×N 矩阵，N=50，其中第i(i ≥1)行、第j(i ≥1)列元素a ij =1i+j −1，右端向量b 的第i(i ≥1)个分量为b i = 10i+j −1N j=1.列主元素Gauss 消去过程中，要用到两种初等行变换。
第一种，交换两行的位置；第二种，用一个数乘某一行加到另一行上。
在第k 次消元之前，先对增广矩阵 A (k),b (k) 作第一种行变换，使得a ik (k)中绝对值最大的元素交换到第k 行的主对角线位置上，然后再使用第二种行变换进行消元。
如此往复，最后得到一个上三角系数矩阵，并回代求解解向量。
由于每次消元前选取了列主元素，因此与顺序Guass 消元法相比，可提高数值计算的稳定性，且其计算量与顺序Guass 消元法相同。
列主元的Gauss 消去法要求系数矩阵A 非奇异。
LU 分解法，即通过一系列初等行变换将系数矩阵A 分解成一个下三角矩阵L 与一个上三角矩阵U 的乘积，进一步通过求解两个三角矩阵得出解向量。
若L 为单位下三角矩阵，U 是上三角矩阵，则称为Doolittle 分解；若L 为下三角矩阵，U 是单位上三角矩阵，则称为Crout 分解。
若系数矩阵A 的前n-1阶顺序主子式不为零，则Doolittle\Crout 分解具有唯一性。
若在每步行变换中选取主元，可提高数值计算稳定性。
本算例中采用选主元的Doolittle 分解。
通过分析可知，本算例中待求解线型方程组系数矩阵为非奇异矩阵，且其前n-1阶顺序主子式不为零。
方程组的解向量为x = 10,10,⋯,10 T 。
满足列主元高斯消去法以及LU 分解法的基本使用条件。
为了验证上述两种方法对本算例的适用性，笔者利用Microsoft Visual C++6.0编写了该算例的列主元高斯消去法以及LU 分解法的程序代码，并进行了运算求解。
二、计算结果及分析1.高斯求解程序运行结果：**********************************各行主元大小*********************************第1行主元大小：1.00000e+000 第2行主元大小：8.33333e-002 第3行主元大小：6.48148e-003 第4行主元大小：1.13636e-003 第5行主元大小：1.30273e-004 第6行主元大小：7.84929e-006 第7行主元大小：-5.19257e-007 第8行主元大小：3.74111e-008 第9行主元大小：-2.11897e-009 第10行主元大小：1.27171e-010 第11行主元大小：-7.11093e-012 第12行主元大小：-2.06444e-013 第13行主元大小：1.41512e-014 第14行主元大小：-1.24316e-015 第15行主元大小：3.50355e-016 第16行主元大小：3.99615e-016 第17行主元大小：3.02521e-016 第18行主元大小：6.93669e-017 第19行主元大小：1.11654e-016***********************************方程组的解**********************************x1：9.999996x2：10.000578 x3：9.977938x4：10.360343 x5：6.904310x6：25.281558 x7：-33.905604 x8：76.093315 x9：-10.091680 x10：-76.424166 x11：108.556590 x12：53.247546 x13：-123.827690 x14：116.991522 x15：-164.520188 x16：291.762560 x17：-166.575110 x18：-8.426663 x19：73.980824 x20：-9.3859802.LU求解程序运行结果：***********************************方程组的解**********************************x1：9.999996x2：10.000578 x3：9.977938x4：10.360343 x5：6.904310x6：25.281558 x7：-33.905604 x8：76.093315 x9：-10.091680 x10：-76.424166 x11：108.556590 x12：53.247546 x13：-123.827690 x14：116.991522 x15：-164.520188 x16：291.762560 x17：-166.575110 x18：-8.426663 x19：73.980824 x20：-9.385980两种数值计算方法结果均与实际解相比产生了较大的误差，这是由于在消元过程中产生了绝对值很小的主元值，以列主元Guass消元法为例，最小主元值为6.93669e-017，这样就造成行乘数绝对值非常大，出现较小数加不到较大数中的现象，舍入误差的积累很大，运算结果完全失真。
为了提高运算精度，笔者采用long double精度，但误差减小程度不显著。
三、结论本方程组在进行高斯消元过程中出现小列主元情况，因而是病态矩阵。
单纯通过提高运算精度并不能有效改善其病态程度，可采用平衡方法、残差校正法对该病态矩阵进行求解。
四、附录（VC++程序代码）1.列主元Guass消去法：#include<stdio.h>#include<math.h>#include<stdlib.h>#define N 20#define e 1e-20//主元绝对值下限void main(){long double A[N][N];long double b[N];long double x[N];long double m,S,temp,max;int i,j,k,col;//************************原始方程的计算和输出************************//方程组原始条件的代入for(i=0;i<N;i++){b[i]=0;for(j=0;j<N;j++){A[i][j]=1/((long double)(i)+(long double)(j)+1);b[i]=b[i]+10*A[i][j];}}printf("*************************************************************\n");printf("**********************高斯列主元消去法***********************\n");printf("*************************************************************\n");//输出原方程组系数矩阵printf("原方程组系数矩阵A：\n");for(i=0;i<N;i++){printf(" 第%d行\n",i+1);for(j=0;j<N;j++){if(j%5==0)printf("第%d列至第%d列：\n",j+1,j+5);printf("%1.5e ",A[i][j]);if((j+1)%5==0)printf("\n");}}//输出原方程组右端列向量printf("原方程组右端列向量b：\n");for(i=0;i<N;i++)printf("第%d行：%1.5e\n",i+1,b[i]);//****************************消去过程********************************* printf("**********************各行主元大小***************************\n"); for(k=0;k<N-1;k++){//查找最大元素所在的行max=fabs(A[k][k]);col=k;for(i=k;i<N;i++){if(max<fabs(A[i][k])){max=fabs(A[i][k]);col=i;}}printf("第%d行主元大小：%1.5e\n",k+1,A[col][k]);//判断主元大小if(max<=e){printf("主元过小，造成舍入误差过大，高斯求解失效！\n");exit(0);}//换行过程for(j=k;j<N;j++){temp=A[col][j];A[col][j]=A[k][j];A[k][j]=temp;}temp=b[col];b[col]=b[k];b[k]=temp;//消元过程for(i=k+1;i<N;i++){m=A[i][k]/A[k][k];for(j=k;j<N;j++)A[i][j]=A[i][j]-m*A[k][j];for(j=0;j<k;j++)A[i][j]=0;b[i]=b[i]-m*b[k];}}//***************************消元结果********************************** printf("****************************消元结果*************************\n");//高斯消元后新的系数矩阵printf("消元后系数矩阵A'：\n");for(i=0;i<N;i++){printf(" 第%d行\n",i+1);for(j=0;j<N;j++){if(j%5==0)printf("第%d列至第%d列：\n",j+1,j+5);printf("%1.5e ",A[i][j]);if((j+1)%5==0)printf("\n");}}//高斯消元后新的右端列向量printf("原方程组右端列向量b'：\n");for(i=0;i<N;i++)printf("第%d行：%1.5e\n",i+1,b[i]);//*****************************回代过程******************************** x[N-1]=b[N-1]/A[N-1][N-1];for(k=N-2;k>=0;k--){S=b[k];for(j=k+1;j<N;j++){S=S-A[k][j]*x[j];}x[k]=S/A[k][k];}//*****************************输出结果******************************** printf("*************************方程组的解**************************\n");for(i=0;i<N;i++)printf("x%d：%f\n",i+1,x[i]);}2.LU分解法：#include<stdio.h>#include<math.h>void main(){long double a[N][N],l[N][N],u[N][N];long double b[N],x[N],y[N],s[N],temp;long double max;int i,j,k,t;int M[N];//*************************原始方程的计算和输出************************ //方程组原始条件的代入for(i=0;i<N;i++){b[i]=0;for(j=0;j<N;j++){a[i][j]=1/((long double)(i)+(long double)(j)+1);b[i]=b[i]+10*a[i][j];}}printf("*************************************************************\n");printf("**********************选主元的Doolittle消去法*****************\n");printf("*************************************************************\n");//输出原方程组系数矩阵printf("原方程组系数矩阵A：\n");for(i=0;i<N;i++){printf(" 第%d行\n",i+1);for(j=0;j<N;j++){if(j%5==0)printf("第%d列至第%d列：\n",j+1,j+5);printf("%1.5e ",a[i][j]);if((j+1)%5==0)printf("\n");}}//输出原方程组右端列向量printf("原方程组右端列向量b：\n");for(i=0;i<N;i++)printf("第%d行：%1.5e\n",i+1,b[i]);//*****************************分解过程******************************** for(k=0;k<N;k++){for(i=k;i<N;i++)s[i]=a[i][k];for(t=0;t<=k-1;t++)s[i]-=l[i][t]*u[t][k];}max=0;for(i=k;i<N;i++){if(max<fabs(s[i])){max=fabs(s[i]);M[k]=i;}}if(M[k]!=k){for(t=0;t<=k-1;t++){temp=l[k][t];l[k][t]=l[M[k]][t];l[M[k]][t]=temp;}for(t=k;t<N;t++){temp=a[k][t];a[k][t]=a[M[k]][t];a[M[k]][t]=temp;}temp=s[k];s[k]=s[M[k]];s[M[k]]=temp;}u[k][k]=s[k];for(j=k+1;j<N;j++){u[k][j]=a[k][j];for(t=0;t<=k-1;t++)u[k][j]-=l[k][t]*u[t][j];}for(i=k+1;i<N;i++)l[i][k]=s[i]/u[k][k];}//输出LU矩阵printf("LU矩阵：\n");for(i=0;i<N;i++){printf(" 第%d行\n",i+1);for(j=0;j<N;j++){if(j%5==0)printf("第%d列至第%d列：\n",j+1,j+5);if(j<i)printf("%1.5e ",l[i][j]);if(j>=i)printf("%1.5e ",u[i][j]);if((j+1)%5==0)printf("\n");}}//求Qbfor(k=0;k<N-1;k++){t=M[k];temp=b[k];b[k]=b[t];b[t]=temp;}//求解y和xy[0]=b[0];for(i=1;i<N;i++){y[i]=b[i];for(t=0;t<=i-1;t++)y[i]-=l[i][t]*y[t];}x[N-1]=y[N-1]/u[N-1][N-1];for(i=N-2;i>=0;i--){x[i]=y[i];for(t=i+1;t<N;t++)x[i]=x[i]-u[i][t]*x[t];x[i]=x[i]/u[i][i];}//输出结果printf("*************************方程组的解**************************\n");for(i=0;i<N;i++)printf("x%d：%f\n",i+1,x[i]);}。