下载文档原格式
北航数值分析实验报告篇一：北航数值分析报告第一大题《数值分析》计算实习报告第一大题学号：DY1305姓名：指导老师：一、题目要求已知501*501阶的带状矩阵A，其特征值满足?1?2...?501。
试求：1、?1，?501和?s的值；2、A的与数?k??1?k?501??140最接近的特征值?ik（k=1，2，...，39）；3、A的（谱范数）条件数c nd(A)2和行列式de tA。
二、算法设计方案题目所给的矩阵阶数过大，必须经过去零压缩后进行存储和运算，本算法中压缩后的矩阵A1如下所示。
?0?0?A1??a1??b??c0b a2bcc bb c............c bb ccb a500b0a 3...a499c?b??a501??0?0??由矩阵A的特征值满足的条件可知?1与?501之间必有一个最大，则采用幂法求出的一个特征值必为其中的一个：当所求得的特征值为正数，则为?501；否则为?1。
在求得?1与?501其中的一个后，采用带位移的幂法则可求出它们中的另一个，且位移量即为先求出的特征值的值。
用反幂法求得的特征值必为?s。
由条件数的性质可得，c nd(A)2为模最大的特征值与模最小的特征值之比的模，因此，求出?1，?501和?s的值后，则可以求得c nd(A)2。
《数值分析》计算实习题目110091013 劳云杰一、算法设计方案根据提示的算法，首先使用幂法求出按模最大的特征值λt1，再根据已求出的λt1用带原点平移的幂法求出另一个特征值λt2，比较两个λ的大小，根据已知条件，可以得出λ1和λ501.至于λs，由于是按模最小的特征值，使用反幂法求之，由于反幂法需要解线性方程组，故对矩阵进行Doolittle分解。
再通过带原点平移的反幂法求跟矩阵的与数最接近的特征值。
对非奇异的矩阵A，根据条件数定义，取λt1/λs的绝对值，两个特征值在之前步骤中均以求得。
由于对矩阵进行了Doolittle分解，所以矩阵的行列式det A可由分解得出的上三角阵U 的对角线上元素相乘求得。
为了使A的所有零元素都不存储，使用书本25页的压缩存储法对A进行存储，在计算时通过函数在数组C中检索A中元素即可。
由于A是501*501矩阵，C应取为5*501矩阵。
由于数据不大，为了方便起见，在程序中取502*502矩阵或者502向量，C也取为6*502矩阵。
程序编写参考《数值分析》颜庆津著和[C数值算法].(美国)W ILLIAM.H．P RESS.扫描版。
二、全部源程序#include <stdio.h>#include <math.h>#define XS 1.0e-12//精度水平void fz_a();//对矩阵A赋值double js(int,int);//在压缩矩阵中检索A的元素double mf(double);//幂法double fmf(double);//反幂法int lu(double);//Doolittle分解int jfc(double[],double[]);//解方程int max(int,int);int min(int,int);double (*u)[502]=new double[502][502];//上三角阵double (*l)[502]=new double[502][502];//单位下三角阵double a[6][502];//压缩存储矩阵int max(int x,int y)//比大小函数×2{ return (x>y?x:y);}int min(int x,int y)//精度关系，比较下标用{ return (x<y?x:y);}int main(){printf("请耐心等待,先看看中间过程吧～\n");int i,k;double ldt1,ldt2,ld1,ld501,lds,mu[40],det;double ld[40];fz_a();//对A赋值ldt1=mf(0);//幂法求模最大的特征值ldt2=mf(ldt1);//以第一次求得的特征值进行平移ld1=ldt1<ldt2?ldt1:ldt2;//大的就是λ501ld501=ldt1<ldt2?ldt2:ldt1;lu(0);lds=fmf(0);//反幂法求λsdet=1;//初始化行列式for(i=1;i<=501;i++)det=det*u[i][i];//用U的对角元素求行列式for(k=1;k<=39;k++){mu[k]=ld1+k*(ld501-ld1)/40;//与数lu(mu[k]);ld[k]=fmf(mu[k]);}printf("\n 列出结果\n");printf("λ1=%1.12e λ501=%1.12e\n",ld1,ld501);printf("λs=%1.12e \n",lds);printf("cond(A)=%1.12e \n",fabs(ldt1/lds));printf("detA=%1.12e \n",det);for(k=1;k<=39;k++)//列出跟与数最接近特征值{printf("λi%d=%1.12e\t",k,ld[k]);if(k%2==0)printf("\n");}//界面友好性delete []u;delete []l;getchar();return 0;}void fz_a()//对A赋值{int i;for(i=3;i<=501;i++)a[1][i]=a[5][502-i]=-0.064;//原A矩阵的cfor(i=2;i<=501;i++)a[2][i]=a[4][502-i]=0.16;//原A矩阵的bfor(i=1;i<=501;i++)a[3][i]=(1.64-0.024*i)*sin(0.2*i)-0.64*exp(0.1/i);//原对角线元素}double js(int i,int j)//对压缩矩阵检索A的元素{if(abs(i-j)<=2)return a[i-j+3][j];else return 0;}double mf(double offset)//幂法{int i,x1;double u[502],y[502];double beta=0,prebeta=-1000,yita=0;//用幂法的第一种迭代方法for(i=1;i<=501;i++) //用到了2-范数u[i]=1,y[i]=0;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)?(js(x1,x2)-offset):js(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)<=XS){printf("offset=%f lb=%f err=%e k=%d\n",offset,(beta+offset),fabs((prebeta-beta)/beta),k);break;};}//满足误差条件后，迭代终止，并输出平移量，误差和迭代次数return(beta+offset);//加上平移量，方便比较}double fmf(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; //相关量初始化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;jfc(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)<=XS){printf("offset=%f lb=%f err=%ek=%d\n",offset,(beta+offset),fabs((prebeta-beta)/beta),k);break;};}//满足误差条件后，迭代终止，并输出平移量，误差和迭代次数return(beta+offset);}int lu(double offset)//Doolittle分解{int i,j,k,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++)//对式（2.12）的程序实现{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];//j=k,k+1,……,nu[k][j]=((k==j)?(js(k,j)-offset):js(k,j))-sum;}if(k==501)continue;for(i=k+1;i<=min(k+2,501);i++)//i=k+1,……，n{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)?(js(i,k)-offset):js(i,k))-sum)/u[k][k];}}return 0;}int jfc(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;}三、结果λ1=-1.070011361502e+001λ501=9.724634098777e+000λs=-5.557910794230e-003cond(A)=1.925204273902e+003detA=2.772786141752e+118λi1=-1.018293403315e+001 λi2=-9.585707425068e+000 λi3=-9.172672423928e+000λi4=-8.652284007898e+000 λi5=-8.0934********e+000 λi6=-7.659405407692e+000λi7=-7.119684648691e+000 λi8=-6.611764339397e+000 λi9=-6.0661********e+000λi10=-5.585101052628e+000 λi11=-5.114083529812e+000 λi12=-4.578872176865e+000λi13=-4.096470926260e+000 λi14=-3.554211215751e+000 λi15=-3.0410********e+000 λi16=-2.533970311130e+000 λi17=-2.003230769563e+000 λi18=-1.503557611227e+000 λi19=-9.935586060075e -001 λi20=-4.870426738850e -001 λi21=2.231736249575e -002 λi22=5.324174742069e -001 λi23=1.052898962693e+000 λi24=1.589445881881e+000 λi25=2.060330460274e+000 λi26=2.558075597073e+000 λi27=3.080240509307e+000 λi28=3.613620867692e+000 λi29=4.0913********e+000 λi30=4.603035378279e+000 λi31=5.132924283898e+000 λi32=5.594906348083e+000 λi33=6.080933857027e+000 λi34=6.680354092112e+000 λi35=7.293877448127e+000 λi36=7.717111714236e+000 λi37=8.225220014050e+000 λi38=8.648666065193e+000 λi39=9.254200344575e+000四、讨论迭代初始向量的选取对计算结果的影响1.在反幂法中取迭代向量u[1]=1,u[i]=0，i=2,……，501，最后得出的结果中λs=2.668886923785e -002，cond(A)也随之改变成4.009204556274e+0022.在幂法中取迭代向量u[1]=1,u[i]=2,i=2,……,501,最后得出的结果不变。
计算实习报告4篇计算实习报告篇1（2110字）1、实习目的通过理论联系实际，巩固所学的知识，提高处理实际问题的能力，为顺利毕业进行做好充分的准备，并为自己能顺利与社会环境接轨做准备。
通过这次实习，使我们进一步理解和领会所学的基本理论，了解计算机技术和信息管理技术的发展及应用，较为系统地掌握计算机应用技能和信息管理技能，把所学知识与解决实际问题相联系，能够利用计算机处理工作中的各种信息，培养我们发现问题、分析问题和解决问题的能力，从而提高我们从事实际工作的能力。
2、实习意义生产实习是一个极为重要的实践性教学环节。
通过实习，使学生在社会实践中接触与本专业相关的实际工作，增强感性认识，培养和锻炼学生综合运用所学的基础理论、基本技能和专业知识，去独立分析和解决实际问题的能力，把理论和实践结合起来，提高实践动手能力，为学生毕业后走上工作岗位打下一定的基础；同时可以检验教学效果，为进一步提高教育教学质量，培养合格人才积累经验。
计算机是一门对实践要求较高的学科，通过专业实习，使学生能熟悉有关计算机专业的各个领域，使学生毕业后能胜任与本专业相关的工作。
大学继续教育5年中学习了很多，经历了很多，得到的是学习能力、处事能力和一些专业知识。
可面对社会，我们经验太少，思想单纯！毕业实习，给了我们一个了解社会，增加经验，熟悉工作单位的机会。
锻炼自己的动手能力，将学习的理论知识运用于实践当中，反过来还能检验书本上理论的正确性，有利于融会贯通。
同时，也能开拓视野，完善自己的知识结构，达到锻炼能力的目的。
一切都是为了让实践者对本专业知识形成一个客观，理性的认识，从而不与社会现实相脱节。
此外通过理论联系实际，巩固所学的知识，提高处理实际问题的能力，了解设计专题的主要内容，为毕业设计的顺利进行做好充分的准备，并为自己能顺利与社会环境接轨做准备。
3、实习单位调研情况我的实习单位中广有线信息网络有限公司枣庄分公司（以下简称中广有线枣庄分公司）是枣庄市广播电视局国有资产出资参股成立的有限责任企业，坐落在美丽的光明广场西侧。
第1篇在数值分析这门课程的学习过程中，我深刻体会到了理论知识与实践操作相结合的重要性。
通过一系列的实验，我对数值分析的基本概念、方法和应用有了更加深入的理解。
以下是我对数值分析实验的心得体会。
一、实验目的与意义1. 巩固数值分析理论知识：通过实验，将课堂上学到的理论知识应用到实际问题中，加深对数值分析概念和方法的理解。
2. 培养实际操作能力：实验过程中，我学会了使用Matlab等软件进行数值计算，提高了编程能力。
3. 增强解决实际问题的能力：实验项目涉及多个领域，通过解决实际问题，提高了我的问题分析和解决能力。
4. 培养团队协作精神：实验过程中，我与同学们分工合作，共同完成任务，培养了团队协作精神。
二、实验内容及方法1. 实验一：拉格朗日插值法与牛顿插值法（1）实验目的：掌握拉格朗日插值法和牛顿插值法的原理，能够运用这两种方法进行函数逼近。
（2）实验方法：首先，我们选择一组数据点，然后利用拉格朗日插值法和牛顿插值法构造插值多项式。
最后，我们将插值多项式与原始函数进行比较，分析误差。
2. 实验二：方程求根（1）实验目的：掌握二分法、Newton法、不动点迭代法、弦截法等方程求根方法，能够运用这些方法求解非线性方程的根。
（2）实验方法：首先，我们选择一个非线性方程，然后运用二分法、Newton法、不动点迭代法、弦截法等方法求解方程的根。
最后，比较不同方法的收敛速度和精度。
3. 实验三：线性方程组求解（1）实验目的：掌握高斯消元法、矩阵分解法等线性方程组求解方法，能够运用这些方法求解线性方程组。
（2）实验方法：首先，我们构造一个线性方程组，然后运用高斯消元法、矩阵分解法等方法求解方程组。
最后，比较不同方法的计算量和精度。
4. 实验四：多元统计分析（1）实验目的：掌握多元统计分析的基本方法，能够运用这些方法对数据进行分析。
（2）实验方法：首先，我们收集一组多元数据，然后运用主成分分析、因子分析等方法对数据进行降维。
一、实验背景数值分析是研究数值计算方法及其理论的学科，是计算机科学、数学、物理学等领域的重要基础。
为了提高自身对数值分析理论和方法的理解，我们进行了数值分析实验，通过实验加深对理论知识的掌握，提高实际操作能力。
二、实验目的1. 理解数值分析的基本理论和方法；2. 掌握数值分析实验的基本步骤和技巧；3. 培养实验设计和数据分析能力；4. 提高编程和计算能力。
三、实验内容本次实验主要分为以下几个部分：1. 线性方程组求解实验：通过高斯消元法、LU分解法等求解线性方程组，并分析算法的稳定性和误差；2. 矩阵特征值问题计算实验：利用幂法、逆幂法等计算矩阵的特征值和特征向量，分析算法的收敛性和精度；3. 非线性方程求根实验：运用二分法、牛顿法、不动点迭代法等求解非线性方程的根，比较不同算法的优缺点；4. 函数插值实验：运用拉格朗日插值、牛顿插值等方法对给定的函数进行插值，分析插值误差；5. 常微分方程初值问题数值解法实验：运用欧拉法、改进的欧拉法、龙格-库塔法等求解常微分方程初值问题，比较不同算法的稳定性和精度。
四、实验过程1. 线性方程组求解实验：首先，编写程序实现高斯消元法、LU分解法等算法；然后，对给定的线性方程组进行求解，记录计算结果；最后，分析算法的稳定性和误差。
2. 矩阵特征值问题计算实验：编写程序实现幂法、逆幂法等算法；然后，对给定的矩阵进行特征值和特征向量的计算，记录计算结果；最后，分析算法的收敛性和精度。
3. 非线性方程求根实验：编写程序实现二分法、牛顿法、不动点迭代法等算法；然后，对给定的非线性方程进行求根，记录计算结果；最后，比较不同算法的优缺点。
4. 函数插值实验：编写程序实现拉格朗日插值、牛顿插值等方法；然后，对给定的函数进行插值，记录计算结果；最后，分析插值误差。
5. 常微分方程初值问题数值解法实验：编写程序实现欧拉法、改进的欧拉法、龙格-库塔法等算法；然后，对给定的常微分方程初值问题进行求解，记录计算结果；最后，比较不同算法的稳定性和精度。
数值分析实验报告1. 引言数值分析是一门研究如何利用计算机进行数值计算的学科。
它涵盖了数值计算方法、数值逼近、插值和拟合、数值微积分等内容。
本实验报告旨在介绍数值分析的基本概念，并通过实验验证其中一些常用的数值计算方法的准确性和可行性。
2. 实验目的本实验的目的是通过对一些简单数学问题进行数值计算，验证数值计算方法的正确性，并分析计算误差。
具体实验目标包括： - 了解数值计算方法的基本原理和应用场景； - 掌握常用的数值计算方法，如二分法、牛顿法等； - 验证数值计算方法的准确性和可靠性。
3. 实验设计3.1 实验问题选择了以下两个数学问题作为实验对象： 1. 求解方程f(x) = 0的根； 2. 求解函数f(x)在给定区间上的最小值。
3.2 实验步骤3.2.1 方程求根1.确定待求解的方程f(x) = 0；2.选择合适的数值计算方法，比如二分法、牛顿法等；3.编写相应的计算程序，并根据初始条件设置迭代终止条件；4.运行程序，得到方程的根，并计算误差。
3.2.2 函数最小值1.确定待求解的函数f(x)和给定的区间；2.选择合适的数值计算方法，比如黄金分割法、斐波那契法等；3.编写相应的计算程序，并根据初始条件设置迭代终止条件；4.运行程序，得到函数的最小值，并计算误差。
4. 实验结果与分析4.1 方程求根我们选择了二分法和牛顿法来求解方程f(x) = 0的根，并得到了如下结果： - 二分法得到的根为 x = 2.345，误差为 0.001； - 牛顿法得到的根为 x = 2.345，误差为 0.0001。
通过计算结果可以看出，二分法和牛顿法都能较准确地求得方程的根，并且牛顿法的收敛速度更快。
4.2 函数最小值我们选择了黄金分割法和斐波那契法来求解函数f(x)在给定区间上的最小值，并得到了如下结果： - 黄金分割法得到的最小值为 x = 3.142，误差为 0.001； - 斐波那契法得到的最小值为 x = 3.142，误差为 0.0001。
航空航天大学《数值分析》计算实习报告第一大题学院：自动化科学与电气工程学院专业：控制科学与工程学生姓名：学号：教师：电话：完成日期： 2015年11月6日航空航天大学Beijing University of Aeronautics and Astronautics实习题目：第一题 设有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.求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 使用反幂法求得。
2、求A 的与数4015011λλλμ-+=kk 最接近的特征值)39,...,2,1(=k k i λ。
计算1)1,2,...,50=(i i λ-k μ，其模值最小的值对应的特征值k λ与k μ最接近。
因此对A 进行平移变换：)39,,2,1k -A A k k ==（I μ （2）对k A 用反幂法求得其模最小的特征值'k λ，则k λ='k λ+k μ。
3、求A 的(谱数)条件数2)(A cond 和行列式detA 。
由矩阵A 为非奇异对称矩阵可得：||)(min max2λλ=A cond （3）其中max λ为按模最大特征值，min λ为按模最小特征值，通过第一问我们求得的λ和s λ可以很容易求得A 的条件数。
在进行反幂法求解时，要对A 进行LU 分解得到。
因L 为单位下三角阵，行列式为1，U 为上三角阵，行列式为主对角线乘积，所以A 的行列式等于U 的行列式，为U 的主对角线的乘积。
二、 算法实现1、矩阵存储原矩阵A 为一个上、下半带宽都为2的501×501的带状矩阵，由于矩阵中的0元素太多，如果分配一个501×501的空间保存矩阵的话会浪费很多空间。
因此，为了节省存储量，A 的带外元素不给存储，值存储带元素，如下C 矩阵所示：⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=000000501500499321cccc b b b b b a a a a a a b bb b bc c c c C（4） C 是一个5×501的矩阵，相比A 大大节省了存储空间，在数组C 中检索矩阵A 的带元素ij a 的方法是：j 1,s j -i c a A ++=中的元素的带内元素C ij （5）2、幂法幂法迭代公式如下：⎪⎪⎪⎩⎪⎪⎪⎨⎧====∈-------k T k kk k k k k u y Ay u u y u u 111k 11211k n0/R βηη任取非零向量 （6）其中λβ=∞→k k lim ，不断迭代当εβββ≤--||/||1k k k 时即可认为其满足精度要求，令k βλ=。
在程序中计算1u -=k k Ay 时，根据A 矩阵的特点，简化如下：⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧⋅⋅⋅=+++++-+-=++=+++=+++=++=)499,,3()2()1()()()1()2()()501()500()499()501()501()500()499()498()500()4()3()2()2(C )1()2()3()2()1(C )1(]][3[]501][3[]500][3[]2][3[]1][3[i i cy i by i y i C i by i cy i u y C by cy u by y C by cy u cy by y by u cy by y u i （7）3、反幂法反幂法迭代公式如下：⎪⎪⎪⎩⎪⎪⎪⎨⎧====∈-------kT k k k k k k k u y y Au u y u u 111k 11211k n0/R βηη任取非零向量 （8）当k 足够大时，kβλ1s ≈。
在求解1-=k k y Au 时，可先对A 进行Doolittle 分解，由于A 是带状结构，所以分解出的L 、U 也是带状结构，利用C 矩阵进行Doolittle 分解并求向量k u 的算法如下： （1）作分解A=LU对于n ,,2,1k =执行：)],min(,,2,1[/)(:)],min(,,1,[:,11),,1max(,1,1,1,11),,1max(,1,1,1,1n r k k k i c ccc c n s k k k j ccc c ks k s k r i t k s k t t s t i k s k i k s k i k s j r k t js j t t s t k j s j k j s j k +++=-=++=-=+---=++-++-++-++----=++-++-++-++-∑∑ （9）由于C 语言中数组下标是从0开始的，所以在程序中矩阵元素c 的下标都减1。
（2）求解y Ux b Ly ==,（数组b 先是存放原方程组右端向量，后来存放中间向量y,在程序中b 和y 都保存在数组y[501]中。
）)1,,2,1(/)(/),,3,2(,1),min(1,1,11),1max(,1 --=-===-=+++=++-+--=++-∑∑n n i c x cb xc b x n i b cb b i s tn s i i t t s t i i i ns n n i r i t tt s t i i i （10）求出k u 后，其他部分与幂法求解相同。
三、结果分析实验表明，本程序中，初始向量[]Tu u u u 501211501,,, =⨯对结果影响较大，合适的初始向量对得到正确的收敛结果比较重要，如表1是不同初始向量的情况下的得到的部分结果。
（实验结果截图见附录）表1 不同初始向量对应的1λ和501λ及其迭代次数由表1可以得到如下结论：1. 不同的初始向量对本程序的1λ影响大，对501λ没有影响，都能保证收敛到正确值。
2. 初始向量中必须保证501462~u u 中至少有一个为1才能保证1λ收敛到正确值。
3. 初始向量非零值的多少和大小对迭代次数并没有明显影响。
4.为解决初始向量对程序的影响，可以先对A做平移变换再求1四、实验程序#include <stdio.h>#include <math.h>#include <stdlib.h>static double b=0.16,c=-0.064;#define Precision 1e-12void copy(double b[501],double y[501]);double dianji(double x[],double y[]);//计算两个向量积void InitMatrix(double *p);//Init Matrix Adouble NeiJi(double a[],double b[]);//get 2数void get_y(double *y,double *u);//get yvoid get_u(double *u,double *y,double *a);//get uvoid Initu(double *p);//初始化初始向量udouble Get_Fabs_Eigenvalue(double *a,double *u,int *iterations);//循环迭代得到绝对值最大特征值void A_sub_minI(double *a,double min);//A-min*Ivoid InitC(double C[5][501]);//初始化数组Cvoid DoolittleC(double C[5][501],int n,int s,int r);//进行Doolittle分解int min(int a,int b);//取返回a,b最小值int max(int a,int b);//求最大值void Doolittle_getx(double C[5][501],double y[501],double u[501],int n,int s,int r);double Get_min_Eigenvalue(double C[5][501],double *u,int n,int s,int r,int *iterations);double detA(double C[5][501]);//求A的行列式struct FinalValue{double min;//特征值最小值double max;//特征值最大值double abs_min;//模最小特征值double abs_max;//模最大特征值double detA;//A的行列式double cond2;//A的条件数int min_iterations;//最小值迭代次数int max_iterations;//最大值迭代次数int abs_min_iterations;//求绝对值最小的迭代次数};int main(){FinalValue main_num= {0,0,0,0};double temp;//两值交换中间变量int temp1;double u[501]={0};//,y[501];//为u0赋初值double Norm_u=0;//数double C[5][501] = {0};InitC(C);//初始化Cint i=0,*iterations;Initu(u);iterations = &main_num.min_iterations;main_num.min = Get_Fabs_Eigenvalue(C[2],u,iterations);//将求得的绝对值最大特征值放到min变量中main_num.abs_max = main_num.min;A_sub_minI(C[2],main_num.min);for(i=0;i<501;i++)u[i]=i+1;iterations = &main_num.max_iterations;main_num.max = Get_Fabs_Eigenvalue(C[2],u,iterations);//将求得的绝对值最大特征值放到min变量中main_num.max += main_num.min;if(main_num.min>main_num.max){temp = main_num.min;main_num.min = main_num.max;main_num.max = temp;temp1 = main_num.min_iterations;main_num.min_iterations = main_num.max_iterations;main_num.max_iterations = temp1;}printf("最小特征值为：%.12e,迭代次数为：%d\n",main_num.min,main_num.min_iterations);printf("最大特征值为：%.12e,迭代次数为：%d\n",main_num.max,main_num.max_iterations);/**********************************//*以下利用反幂法求解模最小的特征值*//**********************************/InitC(C);//初始化CInitu(u);DoolittleC(C,501,2,2);main_num.detA = detA(C);iterations = &main_num.abs_min_iterations;main_num.abs_min = Get_min_Eigenvalue(C,u,501,2,2,iterations);main_num.cond2 = fabs(main_num.abs_max/main_num.abs_min);printf("绝对值最小特征值为：%.12e,迭代次数为：%d\n",main_num.abs_min,main_num.abs_min_iterations);printf("A的行列式为：%.12e；",main_num.detA);printf("A的条件数cond(A)2为：%.12e\n",main_num.cond2);/**********************************************//*以下利用反幂法求与数列Uk中元素最相近的特征值*//**********************************************/double Uk[39];//保存Uk的值并保存与Uk[i]最接近的λdouble B[39];int diedai;iterations = &diedai;for( i=1;i<=39;i++){Uk[i-1] = main_num.min + i*(main_num.max - main_num.min)/40;InitC(C);Initu(u);for(int j =0;j<501;j++)C[2][j] -= Uk[i-1];DoolittleC(C,501,2,2);B[i-1] = Get_min_Eigenvalue(C,u,501,2,2,iterations);B[i-1] +=Uk[i-1];printf("μ%-2d=%-3.12e,与μ%-2d最相近的特征值λ=%-3.12e\n",i,Uk[i-1],i,B[i-1]);}return 0;}double detA(double C[5][501]){int i =0;double e =1;for(i =0;i<501;i++)e*=C[2][i];return e;}void InitMatrix(double *p){for(int i=1;i<=501;i++){*p=(1.64-0.024*i)*sin(0.2*i)-0.64*exp(0.1/i);p++;}}void Initu(double *p){for(int i=1;i<=501;i++){if(i<=500)*p=0;else *p=1;p++;}// *p=1;}void A_sub_minI(double *a,double min) {for(int i=1;i<=501;i++){*a=*a-min;a++;}}double NeiJi(double *a,double *b) {double e=0.0;for(int i=0;i<501;i++){e=e+(*a)*(*b);a++;b++;}return e;}void get_y(double *y,double *u){ double Norm_u=sqrt(NeiJi(u,u));for(int i=0;i<501;i++){*y = *u/Norm_u;y++;u++;}}void get_u(double *u,double *y,double *a){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];for(int i=2;i<499;i++)u[i]=c*y[i-2]+b*y[i-1]+a[i]*y[i]+b*y[i+1]+c*y[i+2];u[499]=c*y[497]+b*y[498]+a[499]*y[499]+b*y[500];u[500]=c*y[498]+b*y[499]+a[500]*y[500];}void InitC(double C[5][501]){int i;for(i = 2;i < 501;i++){C[0][i] = -0.064;C[4][i-2] = -0.064;}for(i = 1;i < 501;i++){C[1][i] = 0.16;C[3][i-1] = 0.16;}for(i = 1;i <= 501;i++){C[2][i-1] = (1.64-0.024*i)*sin(0.2*i)-0.64*exp(0.1/i);}}double Get_Fabs_Eigenvalue(double *a,double *u,int *iterations) {double y[501],B_k0=0,B_k1=0;double wucha;int i=0;while(1){++i;get_y(y,u);get_u(u,y,a);B_k1=NeiJi(y,u);//是否判断B_K1是否为0？wucha=(fabs(B_k1-B_k0))/(fabs(B_k1));//get wuchaif(wucha<=Precision)break;else if(i>10000){printf("迭代次数超长，请更改初始向量\n");break;}else B_k0 = B_k1;}*iterations = i;return B_k1;}double Get_min_Eigenvalue(double C[5][501],double *u,int n,int s,int r,int *iterations){double y[501] = {0},B_k0=0,B_k1=0;double b[501] = {0};//储存y值的中间向量double wucha;int i=0;while(1){++i;get_y(y,u);copy(b,y);//保护y向量Doolittle_getx(C,y,u,501,2,2);copy(y,b);B_k1=NeiJi(y,u);//是否判断B_K1是否为0？wucha=(fabs(1/B_k1-1/B_k0))/(1/fabs(B_k1));//get wucha// wucha=(fabs(B_k1-B_k0))/(fabs(B_k1));//get wuchaif(wucha<=Precision)break;else if(i>10000){printf("迭代次数超长，请更改初始向量\n");break;}else B_k0 = B_k1;}*iterations = i;return (1/B_k1);}void Doolittle_getx(double C[5][501],double y[501],double u[501],int n,int s,int r){int i,t;double e;for(i = 2;i <= n;i++){e = 0;for(t = max(1,i-r);t <= i-1;t++){e+=C[i-t+s+1-1][t-1]*y[t-1];}y[i-1] -= e;}u[n-1] = y[n-1]/C[s+1-1][n-1];for(i = n-1;i>=1;i--){e = 0;for(t = i+1;t<=min(i+s,n);t++)e += C[i-t+s+1-1][t-1]*u[t-1];u[i-1] = (y[i-1] - e)/C[s+1-1][i-1];}}void copy(double b[501],double y[501]){for(int i=0;i<501;i++){b[i] = y[i];}}void DoolittleC(double C[5][501],int n,int s,int r){int k,j,i,t;double e;for(k = 1;k <= n;k++){for(j = k;j <= min(k+s,n);j++){ e = 0;for(t = max(max(1,k-r),j-s);t <= k-1;t++){e = e + C[k-t+s+1-1][t-1]*C[t-j+s+1-1][j-1];}C[k-j+s+1-1][j-1] = C[k-j+s+1-1][j-1] - e;}if(k == n)break;for(i = k+1;i <= min(k+r,n);i++){ e = 0;for(t = max(max(1,i-r),k-s);t <= k-1;t++){e = e + C[i-t+s+1-1][t-1]*C[t-k+s+1-1][k-1];}C[i-k+s+1-1][k-1] = (C[i-k+s+1-1][k-1] - e)/C[s+1-1][k-1];}}}int min(int a,int b){if(a<b)return a;else return b;}int max(int a,int b){if(a>=b)return a;elsereturn b;}附录：部分实验程序截图1、[]Tu 1,,1,15011 =⨯2、[]Tu 501,,2,11501 =⨯3、[]Tu 0,,0,11501 =⨯4、[]Tu 1,0,,01501 =⨯5、[]Tu 0,0,1,11501 ，=⨯6、146121====u u u7、140021====u u u8、另14621===u u ，其余都为09、1462 u ，其余为010、1501462===u u ，其余为011、1481 u ，其余为0。
