问题描述
逆幂法用于计算非奇异阵 A A A的按模最小的特征值和特征向量。
算法思想
实际计算中,可以构造以下迭代格式:
{
L
y
ˉ
k
=
z
k
−
1
U
y
k
=
y
ˉ
k
m
k
=
max
(
y
k
)
i
z
k
=
y
k
/
m
k
(
1
<
=
i
<
=
n
)
\begin{cases} L{{\bar y}_k} = {z_{k - 1}}\\ U{y_k} = {{\bar y}_k}\\ {m_k} = \max {({y_k})_i}\\ {z_k} = {y_k}/{m_k}\\ \end{cases} \quad(1 < = i < = n)
⎩⎪⎪⎪⎨⎪⎪⎪⎧Lyˉk=zk−1Uyk=yˉkmk=max(yk)izk=yk/mk(1<=i<=n)
其中
L
L
L是下三角矩阵,
U
U
U是上三角矩阵,
A
−
1
{A^{ - 1}}
A−1的最大特征值是
1
λ
n
\frac{1}{{{\lambda _n}}}
λn1,若
0
<
∣
λ
n
∣
<
∣
λ
n
+
1
∣
,
k
→
∞
0 < |{\lambda _n}| < |{\lambda _{n + 1}}|,k \to \infty
0<∣λn∣<∣λn+1∣,k→∞
{
m
k
→
1
λ
n
z
k
→
x
n
/
max
(
x
n
)
\begin{cases} {m_k} \to \frac{1}{{{\lambda _n}}}\\ {z_k} \to {x_n}/\max ({x_n}) \end{cases}
{mk→λn1zk→xn/max(xn)
测试数据
已知特征值
12
12
12,初始迭代向量
(
1
,
1
,
1
)
(1,1,1)
(1,1,1)和矩阵
[
4
1
4
1
10
1
4
1
10
]
\quad \begin{bmatrix} 4&1&4\\ 1&{10}&1\\ 4&1&{10} \end{bmatrix}
⎣⎡41411014110⎦⎤
C语言代码
头文件
#ifndef NPOWERWAY_H_INCLUDED
#define NPOWERWAY_H_INCLUDED
#endif // NPOWERWAY_H_INCLUDED
#include <stdio.h>
#include <stdlib.h>
#include<math.h>
#define MaxSize 100
#define e 1e-6
double A[MaxSize][MaxSize];
double L[MaxSize][MaxSize];
double R[MaxSize][MaxSize];
double T[MaxSize][MaxSize];
double Z[MaxSize];
double Z1[MaxSize];
double Z2[MaxSize];
double Z3[MaxSize];
int n;//矩阵尺寸大小
double m;
void InitMatrix()
{
int i,j;
for(i=0;i<n;i++)//计算A-m*E
for(j=0;j<n;j++)
{
if(i==j)
A[i][j]=A[i][j]-m;
}
//
for(i=0;i<n;i++)
for(j=0;j<n;j++)
{
if(i==j)
{
L[i][j]=1;
R[i][j]=A[i][i];
}
if(i<j)
L[i][j]=0;
if(i>j)
R[i][j]=0;
}
}
void Doolittle()
{
int i,j,k;
for(k=0;k<n;k++)
{
for(i=k;i<n;i++)//此处i是从k开始循环,若从0
{//开始循环,那么L上半部分非0
double sum1=0;
for(j=0;j<k;j++)
sum1=sum1+L[k][j]*R[j][i];
R[k][i]=A[k][i]-sum1;
}
for(i=k+1;i<n;i++)//此处i是从k+1开始循环,若从0
{
double sum2=0;
for(j=0;j<k;j++)
sum2=sum2+L[i][j]*R[j][k];
L[i][k]=(A[i][k]-sum2)/R[k][k];
}
}
}
void npowerway1()
{
int i,j,k;
double max,sum,sum1;
for(k=1;k<100;k++)
{
for(i=0;i<n;i++)
Z1[i]=Z[i];//储存上一次迭代结果
for(i=0;i<n;i++)
{//计算Y(k)一拔,本程序中Z2
sum1=0;
for(j=0;j<i;j++)
sum1=sum1+L[i][j]*Z2[j];
sum1=Z1[i]-sum1;
Z2[i]=sum1/L[i][i];
}
for(i=n-1;i>=0;i--)
{//计算Y(k),本程序中Z3
sum=0;
for(j=n-1;j>i;j--)
sum=sum+R[i][j]*Z3[j];
sum=Z2[i]-sum;
Z3[i]=sum/R[i][i];
}
for(i=0;i<n;i++)
Z[i]=Z3[i];
max=Z[0];
for(i=0;i<n;i++) //寻找Z的优势分量
{
if(Z[i]>max)
max=Z[i];
}
printf("第%d次迭代后特征值为:%lf\n",k,1/max+m);
printf("第%d次迭代后归一化Z向量为:\n",k);
for(i=0;i<n;i++) //Z向量归一化处理
{
Z[i]=Z[i]/max;
printf("%lf\n",Z[i]);
}
printf("\n");
int flag=0;
for(i=0;i<n;i++)
if(fabs(Z[i]-Z1[i])<e)
flag++;
if(flag==n)
{
printf("第%d次迭代满足精度!",k);
printf("满足精度的特征值为: %lf\n",1/max+m);
printf("该特征值对应的特征向量为:\n");
for(i=0;i<n;i++)
printf("%lf\n",Z[i]);
break;
}
}
}
void input()
{
int i,j;
printf("请输入矩阵行数:\n");
scanf("%d",&n);
printf("请输入已知特征值:\n");
scanf("%lf",&m);
printf("请输入初始迭代向量:\n");
for(i=0;i<n;i++)
scanf("%lf",&Z[i]);
printf("请输入待分解矩阵A:\n");
for(i=0;i<n;i++)
for(j=0;j<n;j++)
scanf("%lf",&A[i][j]);
}
void print()
{
int i,j;
printf("下三角矩阵L:\n");
for(i=0;i<n;i++)
{
for(j=0;j<n;j++)
printf("%lf ",L[i][j]);
printf("\n");
}
printf("上三角矩阵R:\n");
for(i=0;i<n;i++)
{
for(j=0;j<n;j++)
printf("%lf ",R[i][j]);
printf("\n");
}
printf("\n");
}
void npowerway()
{
input();
InitMatrix();
Doolittle();
npowerway1();
print();
}
源文件
#include <stdio.h>
#include <stdlib.h>
#include"npowerway.h"
#include<windows.h>
int main()
{
npowerway();
return 0;
}