数值分析实验报告_2

数值分析实验报告姓名周茹学号专业数学与应用数学指导老师李建良X Double n方程组的解%=司,%2,…Xj11double1方程组的阶n nInt1double*Gaussdouble Al[N+2][N+2],double Bl[N+2]int n{//N阶方程组Ax=B的高斯消去法求解int ijk;double1;double*X=new double[N+2];double A[N+2][N+2];double B[N+2];for i=l;i=N;i++；for j=l j=N;j++；A[i][j]=Al[i][j]fori=l;i=N;i++B[i]=Bl[i];fork=1;k=n;k++，//顺用消去{if A[k][k]==0break;for i=k+l;i=n;i++{l=A[i][k]/A[k][k];A[i][k]-=A[k][k]*l;B[i]-=B[k]*l;for j=k+l;j=n;j++A[i][j]-=A[k][j]*l;}}X[n]=B[n]/A[n][n];for int〃回代求解{double s=0;for intj=i+l;j=n;j++s=s A[i][j]*X[j];+X[i]=B[i]-s/A[i][i];return X;2,逆幕法求解模块变量引用表变量名数学符号类型数量说明A Doublen*n矩阵A，A=%km Double1000M迭代序列mkDouble nz加V doublen矩阵阶数N nInt1放大倍数t tDouble1tEPS Double1精度J=IOS逆幕法求解按模最小特征值void inverse_power{//int ijk=l;double t=le7;double*z=new double[N+2];double*v=new double[N+2];double ans=-100000;double tmp=-100;fori=l;i=N;i++v[i]=l;for i=l；i=N;i++//A数组元素扩大t倍forj=l;j=N;j++A[i][j]=A[i][j]*t;z=GaussA,v,N;m

[0]=-100000;m[l]=maxz;for i=l;i=N;i++v[i]=z[i]/m[l];whilefabsm[k-1]-m[k-2]=EPS{//判断收敛是否达到要求z=GaussA,v,N;〃高斯消去法解Azk=Vk-l m[k]=maxz;for i=l;i=N;i++v[i]=z[i]/m[k];m[k]=l/m[k];k++;}coutsetiosflagsios::fixedsetw10setprecision7;cout迭代步数为:k-lendl;fori=l;i=k-l;i++coutnm[i]=m[i]endl;cout”按模最小特征殖为“m[k-1]/t;}A,反思与思考误差分析在数值运算中是一个很重要又很复杂的问题，因为每一步的运算都会产生误差，误差产生可能来自于算法自身的不稳定，也有可能来自于计算机中的机器误差，一个具体问题中，往往有上千万次的加减乘除运算，有可能导致误差的累计，最终导致结果的错误，在本题中，误差来源于原始数据数量级跨度较大，所求答案□过小，序列□的单步变化量□远小于设定精度口，导致迭代刚开始就结束了算法，产生了错误的答案针对一个数值算法，我们要思考算法对某些特殊数据的作运算时的正确性，并针对此改进算法本次试验，让我知道了对数值实验的结果分析是一件十分细致又复杂的工作，虽然本题中对结果分析较为简单，但是实际应用中并不一定是这样正确性分析也只是数值实验中的一部分，仅做好这一点也是不够的，我们要从各个方面思考实验中可能出现的问题，这样有助于对知识的理解，也有助于开拓新思路，增长知识误差分析在数值运算中是一个很重要又很复杂的问题，因为每一步的运算都会产生误差，误差产生可能来自于算法自身的不稳定，也有可能来自于计算机中的机器误差，一个具体问题中，往往有上千万次的加减乘除运算，有可能导致误差的累计，最终导致结果的错误，在本题中，误差来源于原始数据数量级跨度较大，所求答案□过小，序列□的单步变化量□远小于设定精度口，导致迭代刚开始就结束了算法，产生了错误的答案针对一个数值算法，我们要思考算法对某些特殊数据的作运算时的正确性，并针对此改进算法本次试验，让我知道了对数值实验的结果分析是一件十分细致又复杂的工作，虽然本题中对结果分析较为简单，但是实际应用中并不一定是这样正确性分析也只是数值实验中的一部分，仅做好这一点也是不够的，我们要从各个方面思考实验中可能出现的问题，这样有助于对知识的理解，也有助于开拓新思路，增长知识线性方程组的数值实验

一、课题名字求解双对角线性方程组

二、问题描述考虑一种特殊的对角线元素不为零的双对角线性方程组（以n=7为例）dia、d a6223d a55d6_di写出一般的n（奇数）阶方程组程序（不要用消元法，因为不用它可以十分方便的解出这个方程组）三.摘要本文提出解三对角矩阵的一种十分简便的方法——追赶法，该算法适用于任意三对角方程组的求解、引言对于一般给定的口，我们可以用高斯消去法求解但是高斯消去法过程复杂繁琐对于特殊的三对角矩阵，如果A是不可约的弱对角占优矩阵，可以将□分解为口，再运用追赶法求解

2....l/

六、程序结构设计L数据流向图矩阵分解为、求解向量得到的解A UL Yx

2.符号引数学符号类型个数作用用表变量名A double1初始给出的矩阵（〃山U double1将A分解后得到的矩阵（乙山L double1将A分解后得到的矩阵Gi为…yJy double1中间变量/i fi…fJf double1初始给出的矩阵X double1在最后存放方程的根X1%2…XnJn nint1矩阵的阶数

3、N-S流程图开始输入、口、口、、N ab c[]f[]计算出、卜y,将其存放在b口、c口和f口v

4、模块构图和分程序1求解以、1和y的过程fori=0;iN-l;i++c[i]=c[i]/b[i];b[i+l]=b[i+l]-a[i]*c[i];ifi==Oelsef[i]=f[i]-a[i-l]*f[i-l]/b[i];2求解工1的过程f[N-l]=f[N-l]-a[N-2]*f[N-2]/b[N-l];fori=N-2;i=0;i—；f[i]=f[i]-c[i]*f[i+U

5、计算程序#includeiostream.h#includemath.h#includestdlib.hvoid mainintN,i;cout“输入对角线元素介数:\n;cinN;double*a=new double[N-l];double*b=new double[N];double*c=new double[N-l];double*f=new double[N];coutinput arraya[]\n;〃初始值为下三角中的元素fori=0;iN-l;i++cina[i];coutinput arrayb[]v vVT;〃初始值为对角元素fori=0;iN;i++cinb[i];coutinput arrayc[],\n,;〃初始值为上三角中的元素fori=0;iN-l;i++cinc[i];coutinput arrayf[]\n;〃初始值为f矩阵中的元素fori=0;iN;i++cinf[i];fori=0;iN-l;i++c[i]=c[i]/b[i];b[i+l]=b[i+l]-a[i]*c[i];ifi==Of[i]=f[i]/b[i];elsef[N-l]=f[N-l]-a[N-2]*f[N-2]/b[N-l];fori=N-2;i=0;i—；f[i]=f[i]-c[i]*f[i+Ucout输出计算结果,\n;，，；COUt fori=0;iN;i++coutf[i]Hcout\n;exitO;

七、计算结果与讨论

1.理论结论在验证程序的正确性时，我先使用的课本上给出的例题2-11Xi-120%2—1—0%310X420_%5_-1将a[i]、b[i]、c[i]、f[i]对应的数依次输入，运行程序口得到的结果课本上的答案一致程序的正确性得以验证仔细观察题目，发现本次课题上所给出的矩阵与标准的三对角矩阵略有不同但是，进一步观察计算过程，发现这并不影响计算结果的正确性因此，这个算法同样能够用来解答本次课题输入与课题中所给的矩阵形式一致的矩阵，并运行程序前人对角线元素个数:Input arrayat]L62000Input arrayb[]L523694Input arrayc[]400285Input arrayf[]L4283渝出计算结果65Cl,

0.6,-

0.8,

2.84198,

0.537037,-

0.0277778,

1.25Press anykey tocontinue.得到正确结果

2、讨论这个追赶法能够解出课题中所给方程的解，并且大大简化了计算的形式但是这个算法是在A为三对角矩阵的基础上编写的对于课题中所给的不对称的双对角矩阵还能找到更为简便的算法本次试验我明白数值计算过程中针对一个题目往往有多个解法，我们在进行数值实验时应多加考虑，找到较为便捷的算法，这对程序的设计，题目计算时间的缩短都大有益处非线性方程求解的数值实验

一、课题名字用迭代法求解方程的根

二、问题描述选取一个函数，证明数值验证迭代公式□,n=O l...z z在单根处是立方收敛的，但在多重根处则是线性收敛的

四、本文给出解非线性方程的一种迭代法，其迭代公

三、摘要式即为题目中所给出的公式,研究其在单根处和多根处的收敛性

五、引言反黑

六、计算结果与讨论在验证程序的正确性时，我采用了书上的例题口,A矩阵有特征值2,3,6,选定口，精度口，运行结果如图口，结果正确，然后更换数据，计算本题结果，运行结果如图口，本题中所给矩阵较为特殊（上三角矩阵），由高等代数知识可知，上三角矩阵的特征值为对角元素，对于本题来说,A的特征值对角元素为-

0.0001,

3.737,

0.006,

4.567,程序解得的答案符合实际，可是迭代步数k为2,这意味着迭代刚开始，序列□就已经收敛到-

0.000L明显不正确通过对A矩阵的观察与分析，我们发现，口-O.OOOL相对于其他矩阵元素绝对值太小，甚至比精度□都要小，因此，我猜想，□过小导致第二步迭代时，变化量□的数量级已经比□的数量级要小很多，在迭代终止条件判断时，错误的终止了迭代但如果重新选择精度，第一，无法确定当口为什么数量级时能够作为终止条件，第二，□过小可能会导致□与□之间的差值小于机器误差，同样会使迭代步数无法正确显示，那么，我们可以考虑将矩阵的元素同时扩大t倍，特征值同时也会扩大t倍,在输出结果时，将答案除以t,精度□不用改变，此处，我选取的参数t=D程序修改后的输出结果如图:口，同时，我们可以将中间量序列□输出，查看收敛过程，Inti]=-

988.7312504L

[2]=-

1000.1905560k

[3]=-

999.9968239卜

[4]=-

1000.0000529

[5]=-

999.

9999991.可以发现迭代步数不等于2,此收敛过程也是正常的，印证了此前的猜想，问题解决!

七、程序主要模块代码本程序采用C++编写,IDE采用Visual Studio2012L高斯消去解方程组模块变量引用表变量名数学符号类型说明A,A1Doublen*n系数矩阵A为顺序消去后的矩阵Al AB1,B Doublen方程组右端项b=（4也，…,b）Tn为顺序消去后的B1B。