数值分析实验报告

数值分析实验报告*名：**学号： ************ 专业：数学与应用数学指导老师：***线性方程组的数值实验一、课题名字：求解双对角线性方程组 二、问题描述考虑一种特殊的对角线元素不为零的双对角线性方程组（以n=7为例）⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡d a da d ad a da d a d 7665544332211⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡x x x x x x x 7654321=⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡b b b b b b b 7654321 写出一般的n （奇数）阶方程组程序（不要用消元法，因为不用它可以十分方便的解出这个方程组） 。

三、摘要本文提出解三对角矩阵的一种十分简便的方法——追赶法，该算法适用于任意三对角方程组的求解。

四、引言对于一般给定的d Ax =，我们可以用高斯消去法求解。

但是高斯消去法过程复杂繁琐。

对于特殊的三对角矩阵，如果A 是不可约的弱对角占优矩阵，可以将A 分解为UL ，再运用追赶法求解。

五、计算公式（数学模型）对于形如⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡---b ac b ac ba cb n nn n n 11122211.........的三对角矩阵UL A =，容易验证U 、L 具有如下形式：⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=u a u au a u n n U (3)3221， ⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=1 (1)11321l ll L 比较UL A =两边元素，可以得到⎪⎩⎪⎨⎧-===la b u c l bu iiii i i111 i=2, 3, ... ,n 考虑三对角线系数矩阵的线性方程组 f Ax = 这里()T nx xx x ...21=,()Tnffff ...21=令y Lx =，则有f Uy = 于是有()⎪⎩⎪⎨⎧-==--u yaf y u f y ii i ii11111*i=2， 3, ... ,n再根据y Lx =可得到⎪⎩⎪⎨⎧-==+xl yx y x i iiinn 1i=n-1, n-2, ... ,1六、程序结构设计1.数据流向图2、符号引用表3、N-S 流程图，4、模块构图和分程序（1）求解u i、l i和y i的过程for(i=0;i<N-1;i++){c[i]=c[i]/b[i];b[i+1]=b[i+1]-a[i]*c[i];if(i==0)f[i]=f[i]/b[i];elsef[i]=(f[i]-a[i-1]*f[i-1])/b[i];}（2）求解x i的过程f[N-1]=(f[N-1]-a[N-2]*f[N-2])/b[N-1];for(i=N-2;i>=0;i--)f[i]=f[i]-c[i]*f[i+1];5、计算程序#include<iostream.h>#include<math.h>#include<stdlib.h>void main(){int N,i;cout<<"输入对角线元素个数："<<'\n';cin>>N;double *a=new double[N-1];double *b=new double[N];double *c=new double[N-1];double *f=new double[N];cout<<"input array a[]"<<'\n'; //初始值为下三角中的元素for(i=0;i<N-1;i++)cin>>a[i];cout<<"input array b[]"<<'\n'; //初始值为对角元素for(i=0;i<N;i++)cin>>b[i];cout<<"input array c[]"<<'\n'; //初始值为上三角中的元素for(i=0;i<N-1;i++)cin>>c[i];cout<<"input array f[]"<<'\n'; //初始值为f矩阵中的元素for(i=0;i<N;i++)cin>>f[i];for(i=0;i<N-1;i++){c[i]=c[i]/b[i];b[i+1]=b[i+1]-a[i]*c[i];if(i==0)f[i]=f[i]/b[i];elsef[i]=(f[i]-a[i-1]*f[i-1])/b[i];}f[N-1]=(f[N-1]-a[N-2]*f[N-2])/b[N-1];for(i=N-2;i>=0;i--)f[i]=f[i]-c[i]*f[i+1];cout<<"输出计算结果："<<'\n';cout<<"("; for(i=0;i<N;i++)cout<<f[i]<<" ,";cout<<")"<<'\n'; exit(0);}七、计算结果与讨论1、理论结论在验证程序的正确性时，我先使用的课本上给出的例题⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡--------2112112112112⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡x x x x x 54321=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡00001 将a[i]、b[i]、c[i]、f[i]对应的数依次输入，运行程序得到的结果课本上的答案一致。

程序的正确性得以验证。

仔细观察题目，发现本次课题上所给出的矩阵与标准的三对角矩阵略有不同。

但是，进一步观察计算过程，发现这并不影响计算结果的正确性。

因此，这个算法同样能够用来解答本次课题。

输入与课题中所给的矩阵形式一致的矩阵，并运行程序：得到正确结果 2、讨论这个追赶法能够解出课题中所给方程的解，并且大大简化了计算的形式。

但是这个算法是在A 为三对角矩阵的基础上编写的。

对于课题中所给的不对称的双对角矩阵还能找到更为简便的算法。

本次试验我明白数值计算过程中针对一个题目往往有多个解法，我们在进行数值实验时应多加考虑，找到较为便捷的算法，这对程序的设计，题目计算时间的缩短都大有益处。

非线性方程求解的数值实验一、课题名字：用迭代法求解方程的根 二、问题描述选取一个函数，证明数值验证迭代公式()()()()()x fx f x f x f x f xxn n n nnnn ''2''122--=+ ，n=0, 1, ...在单根处是立方收敛的，但在多重根处则是线性收敛的。

三、摘要本文给出解非线性方程的一种迭代法，其迭代公式即为题目中所给出的公式，研究其在单根处和多根处的收敛性。

四、引言反幂五、计算结果与讨论在验证程序的正确性时，我采用了书上的例题⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=20101350144---A ，A 矩阵有特征值2,3,6，选定()Tv 1,1,10=，精度310-=ξ，运行结果如图，结果正确，然后更换数据,计算本题结果，运行结果如图，本题中所给矩阵较为特殊（上三角矩阵），由高等代数知识可知，上三角矩阵的特征值 为对角元素，对于本题来说，A 的特征值对角元素为-0.0001,3.737,0.006,4.567，程序解得的答案符合实际，可是迭代步数k 却 为2，这意味着迭代刚开始，序列⎭⎬⎫⎩⎨⎧k m 1就已经收敛到-0.0001，明显不正确。

通过对A 矩阵的观察与分析，我们发现，==11a n λ-0.0001,相对于其他矩阵元素绝对值太小，甚至比精度ξ都要小，因此，我猜想，2λ过小导致第二步迭代时，变化量12m m -的数量级已经比310-=ξ的数量级要小很多，在迭代终止条件判断时，错误的终止了迭代。

但如果重新选择精度，第一，无法确定当ξ为什么数量级时能够作为终止条件，第二，ξ过小可能会导致ξ与12m m -之间的差值小于机器误差，同样会使迭代步数无法正确显示，那么，我们可以考虑将矩阵的元素同时扩大t 倍，特征值同时也会扩大t 倍，在输出结果时，将答案除以t ，精度310-=ξ不用改变，此处，我选取的参数t=710程序修改后的输出结果如图：，同时，我们可以将中间量序列⎭⎬⎫⎩⎨⎧k m 1输出，查看收敛过程， ,可以发现迭代步数不等于2，此收敛过程也是正常的，印证了此前的猜想，问 题解决！六、程序主要模块代码本程序采用c++编写，IDE 采用Visual Studio 20121,高斯消去解方程组模块 变量引用表 变量名数学符号 类型 数量 说明A,A1 ()n n j i a A ⨯=,Double n*n 系数矩阵AA1为顺序消去后的矩阵A B1,B T n b b b b ),...,,(21=Double n 方程组右端项B1为顺序消去后的B X ()Tn x x x x ,...,,21=Double n 方程组的解 l ldouble1nnInt1方程组的阶double* Gauss (double A1[N+2][N+2],double B1[N+2],int n) {//N 阶方程组Ax=B 的高斯消去法求解 int i,j,k; double l; double* X=new double[N+2]; double A[N+2][N+2]; double B[N+2];for (i=1;i<=N;i++)for (j=1;j<=N;j++)A[i][j]=A1[i][j];for(i=1;i<=N;i++) B[i]=B1[i];for(k=1;k<=n;k++)//顺序消去{if (A[k][k]==0) break;for (i=k+1;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+1;j<=n;j++)A[i][j]-=A[k][j]*l;}}X[n]=B[n]/A[n][n];for (int i=n-1;i>=1;i--)//回代求解{double s=0;for (int j=i+1;j<=n;j++)s=s+A[i][j]*X[j];X[i]=(B[i]-s)/A[i][i];}return X;}2,逆幂法求解模块变量引用表void inverse_power(){//逆幂法求解按模最小特征值int i,j,k=1;double t=1e7;double *z=new double[N+2];double *v=new double[N+2];double ans=-100000;double tmp=-100;for(i=1;i<=N;i++) v[i]=1;for (i=1;i<=N;i++)//A数组元素扩大t倍for(j=1;j<=N;j++)A[i][j]=A[i][j]*t; z=Gauss(A,v,N); m[0]=-100000; m[1]=max(z);for (i=1;i<=N;i++) v[i]=z[i]/m[1];while(fabs(m[k-1]-m[k-2])>=EPS){//判断收敛是否达到要求 z=Gauss(A,v,N);//高斯消去法解Az(k)=V(k-1) m[k]=max(z); for (i=1;i<=N;i++) v[i]=z[i]/m[k]; m[k]=1/m[k]; k++;}cout<<setiosflags(ios::fixed)<<setw(10)<<setprecision(7); cout<<"迭代步数为："<<k-1<<endl; for(i=1;i<=k-1;i++) cout<<"m["<<i<<"]="<<m[i]<<endl;cout<<"按模最小特征值为："<<m[k-1]/t;}七、反思与思考误差分析在数值运算中是一个很重要又很复杂的问题，因为每一步的运算都会产生误差，误差产生可能来自于算法自身的不稳定，也有可能来自于计算机中的机器误差，一个具体问题中，往往有上千万次的加减乘除运算，有可能导致误差的累计，最终导致结果的错误，在本题中，误差来源于原始数据数量级跨度较大，所求答案n λ过小，序列⎭⎬⎫⎩⎨⎧k m 1的单步变化量111--k k m m 远小于设定精度310-=ξ，导致迭代刚开始就结束了算法，产生了错误的答案。