《数值分析》上机实习报告

《数值分析》上机实习报告

数值分析程序设计作业

一.housholder变换.求方程的解

二.求特征值

一,程序要求:

（1）用Household变换，把A化为三对角阵B(并打印B).

（2）用超松弛法求解Bx=b(取松弛因子ω=1.4,X(0)=0,迭代9次)

（3）用列主元消去法求解Bx=b。

(4) 用对分法求B(即A)的小于23且最接近23的特征值的近似值(误差小于0.1)

(5) 用反幂法求B的该特征值的更精确近似值及相应的特征向量.

二，解题算法

1：Household 法

（1） 令A0=A, aij(1)=aij，已知Ar-1 即Ar-1=(aij(r))

（2） Sr=( (air(r ))2)1/2

（3） αr=Sr2+|a(r)r+1,r|Sr

ur=[0,0,a(r)r+1,r+Sign(a(r)r+1,r)Sr,a(r)r+2,r,…,anr(r)]T

（4） yr=Ar-1ur/αr

（5） kr= urTyr/2αr

（6） qr=yr-krur

（7） Ar=Ar-1-(qrurT+urqrT) r=(1,2,……n-2)

2：超松弛法

其基本思想是迭代，即在高斯方法已求出x(m),x(m-1)的基础上，组合新的序列，从而加快收敛速度。其算式是：

xi(m)=(1-ω)xi(m-1)+ω( bijxi(m)+ xj(m-1)+gi)

其中ω是超松弛因子，当ω>1时，可以加快收敛速度。

3：高斯消去法：

对矩阵作恰当的调整，选取绝对值尽量大的元素作为主元素。然后把矩阵化

为上三角阵，再进行回代，求出方程的解。

对分法和反幂法求解特征值及特征向量

4. 对分法:对于对称的三对角方阵,其特征值必落在[-‖B‖, ‖B‖]中,其中‖B‖是其任意范数.通过同号数的计算,可知在[0, ‖B‖]中特征根的个数,在不断对分,就可得特征根所在的区间,如在[lk,uk]中有一根λk,且[lk,uk]很小,就可取lk+uk/2作为λk的近似值.

5. 反幂法:基于乘幂法中求A-1不太容易,因此不直接用A-1xk=xk+1,而使用Axk+1=xk.用求解线性方程组求出xk+1,再标准化,就可得到绝对值最小的特征值与相应的特征向量,这就是反代法(也称反幂法).

三,程序清单

#include "math.h"

#include "io.h"

《数值分析》上机实习报告

#include "stdio.h"

#include "string.h"

#define N 9

float sign(float x);

void guss(float a[N][N]);

float fanmi(float a[N][N]);

float sor(float a[N][N],float b[N]);

float duifen();

void housholder();

int f(float x);

float e[N],o[N];

float approx;

float

b[N]={2.1874369,33.992318,-25.173417,0.84671695,1.784317,-86.612343,1.1101230,4.719345,-5.6784392};

float a[N][N]=

{

12.38412,2.115237,-1.061074,1.112336,-0.113584,0.718719,1.742382,3.067813,-2.031743,

2.115237,19.141823,-3.125432,-1.012345,2.189736,1.563849,-0.784165,1.112348,3.123124,

-1.061074,-3.125432,15.567914,3.123848,2.031454,1.836742,-1.056781,0.336993,-1.010103,

1.112336,-1.012345,3.123848,27.108437,4.101011,-3.741856,2.101023,-0.71828,-0.037585,

-0.113584,2.189736,2.031454,4.101011,19.897918,0.431637,-3.111223,2.121314,1.784317,

0.718719,1.563849,1.836742,-3.741856,0.431637,9.789365,-0.103458,-1.103456,0.238417,

1.742382,-0.784156,-1.056781,2.101023,-3.111223,-0.103458,14.713846,3.123789,-2.213474,

3.067813,1.112348,0.336993,-0.71828,2.121314,-1.103456,3.123789,30.719334,4.446782,

-2.031743,3.123124,-1.010103,-0.037585,1.784317,0.238417,-2.213474,4.446782,40.00001

《数值分析》上机实习报告

};

void housholder()

{

float Ur[N],Yr[N],qr[N];

int i,j,l,k;

float sr1=0,sr=0.0,ar=0.0,*p,kr=0.0;

for(i=0;i<N-2;i++)

{

for(j=i+1;j<N;j++)

sr1+=a[j][i]*a[j][i];

sr=sqrt(sr1);

ar=sr1+fabs(a[i+1][i])*sr;

for(l=0;l<N;l++)

{

if (l<=i) Ur[l]=0.0;

else if(l==i+1) Ur[l]=a[l][i]+sign(a[l][i])*sr;

else Ur[l]=a[l][i];

}

for(k=0;k<N;k++) { Yr[k]=0; }

for(k=0;k<N;k++)

for(l=0;l<N;l++) Yr[k]+=a[k][l]/ar*Ur[l];

for(k=0;k<N;k++) kr+=Ur[k]*0.5/ar*Yr[k];

for(k=0;k<N;k++) qr[k]=Yr[k]-kr*Ur[k];

for(k=0;k<N;k++)

for(l=0;l<N;l++) a[k][l]-=(qr[k]*Ur[l]+Ur[k]*qr[l]);

sr1=0;kr=0;ar=0.0;

}

system("cls");

for(k=0;k<N;k++)

{

for(l=0;l<N;l++) printf("%.4f ,",a[k][l]);

printf("\n");

}

for(i=0;i<N;i++) o[i]=a[i][i];

for(i=0;i<N-1;i++) e[i]=a[i][i+1];

}

float sor(float a[N][N],float b[N]) /*超松驰法求解的子函数*/

{

float a1[N][N];

int i,j,k,m;

float x[N],temp[N][N];

float w=1.4,h=0,g=0;

《数值分析》上机实习报告

for(m=0;m<N;m++) x[m]=0;

for(i=0;i<N;i++)

for(j=0;j<N;j++) a1[i][j]=a[i][j];

for(i=0;i<N;i++)

for(j=0;j<N;j++) temp[i][j]=0.0;

for(i=0;i<N;i++)

for(j=0;j<N;j++) temp[i][j]=-a1[i][j]/a1[i][i];

for(i=0;i<N;i++) b[i]/=a1[i][i];

for(k=0;k<N;k++)

{

for(i=1;i<=N;i++)

{

for(j=1;j<=N;j++)

{

if(j<=i-1) h+=temp[i-1][j-1]*x[j-1];

else if(j==i) h+=0 ;

else h+=temp[i-1][j-1]*x[j-1];

}

g=(h+b[i-1])*w;

x[i-1]=(1-w)*x[i-1]+g;

h=0.0;g=0.0;

}

}

printf("超松驰法结果：:\n");

printf("\n");

for(i=0;i<N;i++)

printf("%f,",x[i]); /*输出超松驰法求解的特征向量*/

return 0;

}

float sign(float x) /*子函数符号函数*/

{

if (x>0.0) return 1.0;

else if (x<0.0) return -1.0;

else return 0.0;

}

void guss(float a[N][N]) /*高斯法求解的子函数*/

{

float a2[N][N];

int i,j,k;

float l[N]={0,},x[N],h=0;

float

b[N]={2.1874369,33.992318,-25.173417,0.84671695,1.784317,-86.612343,1.1101230,4.719345,-5.6784392};

for(i=0;i<N;i++)

《数值分析》上机实习报告

for(j=0;j<N;j++) a2[i][j]=a[i][j];

for(k=0;k<N;k++)

{

for(j=k+1;j<N;j++) l[j]=a2[j][k]/a2[k][k];

for( …… 此处隐藏：11912字，全部文档内容请下载后查看。喜欢就下载吧 ……