第一题 矩阵变换和矩阵方程求解
一.程序要求: 已知
A=12.384120 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.718280 -0.037585
-0.113584 2.189736 2.031454 4.101011 19.897919 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.784165 -1.056781 2.101023 -3.111223 -0.103458 14.713846 3.123789 -2.213474
3.067813 1.112348 0.336993 -0.718280 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
b=(2.1874369 33.992318 -25.173417 0.84671695 1.784317 -86.612343 1.1101230 4.719345,-5.6784392);
1.用Householder变换，把A化为三对角阵B(并打印B).
2.用超松弛法求解Bx=b(取松弛因子ω=1.4,X(0)=0,迭代9次)
3.用列主元素消去法求解Bx=b。
二．解题算法
1.HOUSEHOLDER算法步骤：
⑴ 令A0=A, aij(1)=aij，已知Ar-1 即Ar-1=(aij(r))
⑵ Sr=( (air(r ))2)1/2
⑶ α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
⑷ yr=Ar-1ur/αr
⑸ kr= urTyr/2αr
⑹ qr=yr-krur
⑺ Ar=Ar-1-(qrurT+urqrT)， r=(1,2,…，n-2)
2. SOR解题算法
其基本思想是在高斯方法已求出x(m),x(m-1)的基础上，经过重新组合的新序列，而此新序列收敛速度加快。其算式是：
xi(m)=(1-ω)xi(m-1)+ω( bijxi(m)+ xj(m-1)+gi)
其中ω是超松弛因子，当ω>1时，可以加快收敛速度。
3. 列主元素消去法程序算法
对矩阵作恰当的调整，选取绝对值尽量大的元素作为主元素。然后进行行变换，把矩阵化为上三角阵，再进行回代，求出方程的解。
三．程序清单
#include <math.h>
#define N 9
main()
{
float e[N+1],f[N+1],x1[N+1],A0[N][N],B[N][N],y[N],u[N],w[N],x0[N];
float x[N],q[N],q1[N],u1[N],ma,a1[N],s,s2,k,b1[N];
static float b[N]={2.1874369,33.992318,-25.173417,0.8471695,1.784317,
-86.612343,1.1101230,4.719345,-5.56784392};
int i,j,r,sign;
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.784163,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.742384,-0.784165,-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}};

clrscr();
for(i=0;i<N;i++)
for(j=0;j<N;j++)
A0[i][j]=A[i][j];
for(r=0;r<N-2;r++)
{
s=0;
for(i=r+1;i<N;i++)
s=s+A[i][r]*A[i][r];
s=sqrt(s);
ma=s*s+fabs(A[r+1][r])*s;
if(A[r+1][r]>0)
sign=1;
if(A[r+1][r]<0)
sign=-1;
for(i=0;i<N;i++)
{if(i<=r)
u[i]=0;
else
{
if(i==r+1)
u[i]=A[r+1][r]+sign*s;
else
u[i]=A[i][r];
}
}
for(i=0;i<N;i++)
{
y[i]=0;
for(j=0;j<N;j++)
y[i]=y[i]+A[i][j]*u[j];
y[i]=y[i]/ma;
}
k=0;
for(i=0;i<N;i++)
k=k+u[i]*y[i];
k=k/(2*ma);
for(i=0;i<N;i++)
q[i]=y[i]-k*u[i];
for(i=0;i<N;i++)
for(j=0;j<N;j++)
A[i][j]=A[i][j]-q[i]*u[j]-q[j]*u[i];
}
for(i=0;i<N;i++)
for(j=0;j<N;j++)
B[i][j]=A[i][j];
for(i=0;i<N;i++){
for(j=0;j<N;j++)
printf("%8.4f",B[i][j]);
printf("\n"); }
for(i=0;i<N;i++)
{x0[i]=0;
x[i]=0;}
for(i=0;i<N;i++){
x[0]=-(B[0][1]/B[0][0])*x[1]+b[0]/B[0][0];
x[0]=1.4*x[0]-0.4*x0[0];
for(j=1;j<N-1;j++){
x[j]=-(B[j][j-1]/B[j][j])*x[j-1]-(B[j][j+1]/B[j][j])*x0[j+1]+b[j]/B[j][j];
x[j]=1.4*x[j]-0.4*x0[j];}
x[N-1]=-(B[N-1][N-2]/B[N-1][N-1])*x0[N-2]+b[N-1]/B[N-1][N-1];
x[N-1]=1.4*x[N-1]-0.4*x0[N-1];
for(j=0;j<N;j++)
x0[j]=x[j];}
printf("\n x:\n");
for(i=0;i<N;i++)
printf("x[%d]=%10.6f\n",i,x[i]);
for(i=0;i<N-1;i++){
b1[i+1]=b[i];
a1[i+2]=B[i+1][i];
e[i+1]=B[i][i+1];
b1[i+1]=B[i][i];}
b1[N]=B[N-1][N-1];
q1[0]=0;
u1[0]=0;
a1[1]=0;
f[N]=b[N-1];
b1[N]=B[N-1][N-1];
for(i=1;i<N;i++)
{
q1[i]=-e[i]/(b1[i]+a1[i]*q1[i-1]);
u1[i]=(f[i]-a1[i]*u1[i-1])/(b1[i]+a1[i]*q1[i-1]);
}
u1[N]=(f[N]-a1[N]*u1[N-1])/(b1[N]+a1[N]*q1[N-1]);
x1[N]=u1[N];
for(i=N-1;i>0;i--)
x1[i]=q1[i]*x1[i+1]+u1[i];
printf("\n x:\n");
for(i=1;i<N+1;i++)
printf("x[%d]=%10.6f\n",i,x1[i]);
}
四．运行结果:
1．以下是HOUSEHOLDER变换结果:
12.38412 -4.89308 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
-4.89308 25.39841 6.49410 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
0.00000 6.49410 20.61151 8.24392 -0.00000 -0.00000 -0.00000 0.00000 0.00000
0.00000 0.00000 8.24392 23.42285 -13.88008 0.00000 0.00000 0.00000 0.00000
0.000000.00000 -0.00000 -13.88008 29.69825 4.53450 0.00000-0.00000 0.00000
0.00000 0.00000 -0.00000 0.00000 4.53450 16.00611 4.88143 0.00000 -0.00000
0.00000 0.00000 -0.00000 0.00000 0.00000 4.88143 26.01332 -4.50362 0.00000
0.00000 0.00000 0.00000 0.00000 -0.00000 0.00000 -4.50362 21.25407 4.50448
0.00000 0.00000 0.00000 0.00000 0.00000 -0.00000 0.00000 4.50448 14.53409
2．以下是列主元素消去法的结果:
x[0]=1.075800
x[1]=2.275745
x[2]=-2.855516
x[3]=2.293101
x[4]=2.112639
x[5]=-6.423842
x[6]=1.357923
x[7]=0.634243
x[8]=-0.587266
3．以下是超松弛法结果:
x[0]=1.075800
x[1]=2.275745
x[2]=-2.855516
x[3]=2.293101
x[4]=2.112639
x[5]=-6.423842
x[6]=1.357923
x[7]=0.634243
x[8]=-0.587266
五，问题讨论
1.SOR方法的矩阵形式为：
X(m)=(E-ωL)-1((1-ω)E+ωR)x(m-1)+(E-ωL)-1ωg
若记 Lω=(E-ωL)-1((1-ω)E+ωR)，SOR收敛的充要条件是S(Lω)<1.且若A为对称正定阵，则当松弛因子ω满足0<ω<2时，SOR方法收敛。
此题中矩阵B是三对角阵，且对称正定，所以选择合适的松驰因子ω，收敛速度是很快的。其最佳松弛因子可用下式计算：
ω=2/[1+ 2 ]
此题根据要求选择ω=1.4。
2．列主元素消去法的优点：
① 由于选主元，使|lij|最小，这样去乘方程的每一系数时，系数的舍入误差不至扩大，并防止溢出与停机。
② 由于选主元，回代时作除数主元aii(i)的绝对值也是最大,这样扩大的误差也是最小.
③ 若detA≠0,则选主元后的主元aii(i) ≠0,这是因为若aii(i)=0则必有aji(i)=0 (j>i)这样按行列式的Laplace展开式就有detA=0而矛盾,因此用主元消去法中进行下去,不止中断停机.同时也由于选主元, aii(i)接近零的概率减少.
④ 列主元素消去法是稳定.


第二题 求特征值以及特征向量

一.程序要求：
用对分法求B(即A)的小于23且最接近23的特征值的近似值(误差小于0.1)在用反幂法求B的该特征值的更精确近似值及相应的特征向量.
二.程序算法
1. 对分法:对于对称的三对角方阵,其特征值必落在[-‖B‖, ‖B‖]中,其中‖B‖是其任意范数.通过同号数的计算,可知在[0, ‖B‖]中特征根的个数,在不断对分,就可得特征根所在的区间,如在[lk,uk]中有一根λk,且[lk,uk]很小,就可取lk+uk/2作为λk的近似值.
2.反幂法:基于乘幂法中求A-1不太容易,因此不直接用A-1xk=xk+1,而使用Axk+1=xk.用求解线性方程组求出xk+1,再标准化,就可得到绝对值最小的特征值与相应的特征向量,这就是反代法(也称反幂法).
三.程序清单
1. 对分法:
#include<stdio.h>
#include<math.h>
#define N 9
int F(float X)
{int n=0,i;
double s[N],w;
float B[N][N]={{12.384120,-4.893077,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000},
{-4.893077,25.398416,6.494097,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000},
{0.000000,6.494097,20.611499,8.243925,0.000000,0.000000,0.000000,0.000000,0.000000},
{0.000000,0.000000,8.243925,23.422838,-13.880071,0.000000,0.000000,0.000000,0.000000},
{0.000000,0.000000,0.000000,-13.880071,29.698278,4.534502,0.000000,0.000000,0.000000},
{0.000000,0.000000,0.000000,0.000000,4.534502,16.006117,4.881435,0.000000,0.000000},
{0.000000,0.000000,0.000000,0.000000,0.000000,4.881435,26.013315,-4.503635,0.000000},
{0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,-4.503635,21.254061,4.504498},
{0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,4.504498,14.534122}};

s[0]=B[0][0]-X;
s[1]=(B[1][1]-X)-B[0][1]*B[0][1]/(B[0][0]-X);
for(i=2;i<N;i++)
{w=s[i-1]*s[i-2];
if(w!=0)
s[i]=B[i][i]-X-B[i-1][i]*B[i-1][i]/s[i-1];
if(s[1]==0)
s[i]=B[i][i]-X;
if(s[i-1]==0)
s[i]=-1;
}
for(i=0;i<N;i++)
if(s[i]>=0)
n++;
return(n);
}

void main()
{int F(float X);
int M=0;
double X1=0.0,X2=23.0,T1;
while(fabs(X2-X1)>0.1)
{T1=(X1+X2)/2.0;
if((F(T1)-F(X2))>=1)
X1=T1;
if((F(T1)-F(X2))==0)
X2=T1;
if((F(T1)-F(X2))>=0)
T1=(X1+X2)/2.0;
M++;
}
printf("T1=%f\n",T1);
printf("M=%d\n",M);
}
2.反幂法:
#include<stdio.h>
#include<math.h>
#define N 9
main()
{float B[N][N]={{12.384120,-4.893077,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000},
{-4.893077,25.398416,6.494097,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000},
{0.000000,6.494097,20.611499,8.243925,0.000000,0.000000,0.000000,0.000000,0.000000},
{0.000000,0.000000,8.243925,23.422838,-13.880071,0.000000,0.000000,0.000000,0.000000},
{0.000000,0.000000,0.000000,-13.880071,29.698278,4.534502,0.000000,0.000000,0.000000},
{0.000000,0.000000,0.000000,0.000000,4.534502,16.006117,4.881435,0.000000,0.000000},
{0.000000,0.000000,0.000000,0.000000,0.000000,4.881435,26.013315,-4.503635,0.000000},
{0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,-4.503635,21.254061,4.504498},
{0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,4.504498,14.534122}};

int i,j,k;
float max,T1,T2,q[N],u[N],y[N],x[N];
float z[N]={1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0};
T1=21.876953; /* 对分法所求得的特征值*/
for(i=0;i<N;i++)
B[i][i]=B[i][i]-T1;
for(k=1;k<=N+1;k++)
{q[0]=-B[0][1]/B[0][0];
u[0]=z[0]/B[0][0];
for(i=1;i<N-1;i++)
q[i]=-B[i][i+1]/(B[i][i]+B[i][i-1]*q[i-1]);
for(i=1;i<N;i++)
u[i]=(z[i]-B[i][i-1]*u[i-1])/(B[i][i]+B[i][i-1]*q[i-1]);
y[N-1]=u[N-1];
for(i=7;i>=0;i--)
y[i]=q[i]*y[i+1]+u[i];
max=0;
for(i=0;i<N;i++)
if(fabs(max)<=fabs(y[i]))
max=y[i];
for(i=0;i<N;i++)
z[i]=y[i]/max;
}
T2=T1+1.0/max;
printf("T2=%f\n",T2);
for(i=0;i<N;i++)
x[i]=z[i];
for(i=0;i<N;i++)
printf("x[%d]=%9.6f\n",i,z[i]);
}
四.运行结果
1．对分法求强特征值为:
T1=21.918365
2．用反幂法求特征值为:
T2=21.918293
特征向量为:
x[0]=0.157189
x[1]=-0.306283
x[2]=0.282570
x[3]=0.286065
x[4]=0.198838
x[5]=0.534490
x[6]=0.462646
x[7]=1.000000
x[8]=0.610016
五.问题讨论
1.对分法简单可靠,数值稳定性较高,对于求少量几个特征值特别适宜.但收敛速度较慢.
2.反幂法是结合对分法使用的,所以近似值较恰当,其收敛速度是惊人的.只要迭代两次就可得到较满意的结果.但在运用中需把一般矩阵化为Hessebberg阵才可计算.