幂法 数值分析

int i, j; //生成上带宽矩阵 a_U 的非 0 元素 for ( i = 0; i < s; ++i ) {
a_U[i] = new double [n-1-i]; } for ( j = 0; j < n-1; ++j ) {
a_U[0][j] = 0.16; } for ( j = 0; j < n-2; ++j ) {
a_U[1][j] = -0.064; }
//生成 A 的对角线向量 a_D
for ( i = 0; i < n; ++i )
{
a_D[i]
=
(1.64
-
0.024*(i+1))*sin( 0.2*(i+1) )-0.64*exp( 0.1/(i+1) )
+ u_max;
}
//生成下带宽矩阵 a_L 的非 0 元素
cout<<"x["<<i<<"]="<<x[i]<<endl<<endl; }
//释放动态分配的内存空间 delete u1; delete x; delete u; delete c; delete a_L; delete a_U; delete y; delete u0; delete a_D; }
for ( j = 0; j < s - i; ++j ) c[i][j] = 0;
} //将 A 的下对角元素赋给 C for ( i = 1; i < r + 1; ++i ) {
for ( j = 0; j < n-i; ++j ) c[s+i][j] = a_L[i-1][j];
for ( j = n - i; j < n; ++j) c[s+i][j] = 0;
cout<<endl;
//求 A 的条件数
double a_cond = cond( u_max, u_min );
cout<<"A
Hale Waihona Puke Baidu
的
条
件
数
为:"<<a_cond<<endl<<endl;
//求与数 u[k]最接近的 A 的特征值 double *u, *x, *u1, cc; cc = (a_max - a_min)/40; const int k = 39; u = new double[k]; x = new double[k]; u1 = new double[n]; for ( i = 0; i < k; ++i )
一、算法设计方案：
int i = 0, j = 0;
1. 求出 A 的压缩后的矩阵 C，并对 C 作 LU 分 解；
2. 根据反幂法求出λ s 的值； 若矩阵 A 存在相等的模最小的特征值，则当 这些特征值完全相等时，反幂法可以求出此 特征值和此特征值对应的特征向量。当这些 特征值互为相反数时，由反幂法的推导过程 可知：此时亦能求出这些特征值的模，但求 不出与这些特征值对应的特征向量。综上所 述，反幂法可以求出矩阵 A 模最小的特征值。 同理可知：幂法可以求出矩阵 A 模最大的特 征值。
a_D = new double [n]; u0 = new double [n]; y = new double [n]; a_U = new double *[s]; a_L = new double *[r]; c = new double *[r+s+1];
//生成压缩矩阵 for ( i = 0; i != r+s+1; ++i )
{ u[i] = a_min + (i + 1)*cc; cout<<u[i]<<endl; x[i] = 0.0; arr_ger( a_U, a_D, a_L, s, r, n, -1*u[i] ); arr_compress( c, a_D, a_U, a_L, s, r, n); for ( j = 0; j < n; ++j ) { u0[j] = 5; } u0[0] = 1; x[i] = inv_power( c, u0, y, s, r, n );//y 为
int i, j, k; //将 A 的主对角元赋给 C for ( j = 0; j < n; ++j )c[s][j] = a_D[j]; //将 A 的上对角元素赋给 C for ( i = 0; i < s; ++i ) {
for ( j = s - i, k = 0; j < n; ++j, ++k ) c[i][j] = a_U[s-i-1][k];
//将带状矩阵 A 压缩为矩阵 C arr_compress( c, a_D, a_U, a_L, s, r, n);
//求出 A 最大的特征值 double a_max = power( c, u0, s, r, n ) fabs(u_max);//a_max 为 A 最大特征值 cout<<"A 的 最 大 特 征 值 为:"<<a_max<<endl; cout<<endl;
int i; double A = 1.0; for ( i = 0; i < n; ++i ) {
A *= c[s][i]; } return A; }
//6.此函数为求出模最大的特征值 #include "math.h" //其中,u 为特征向量,y 亦为特征向量,u=Ay double eig_max( double *y, double *u, int n ) {
//将带状矩阵 A 压缩为矩阵 C arr_compress( c, a_D, a_U, a_L, s, r, n);
//求出 A 最小的特征值
double a_min = inv_power( c, u0, y, s, r, n ) fabs(u_max);//a_min 为 A 最小特征值
cout<<"A 的最小特征值为:"<<a_min<<endl;
//将带状矩阵 A 压缩为矩阵 C arr_compress( c, a_D, a_U, a_L, s, r, n);
LU_decomp( c, s, r, n );
cout<<"A
的
行
列
式
为:"<<setiosflags( ios::scientific )<<setprecision(1
2)<<deta( c, s, r, n )<<endl;
为:"<<u_min<<endl; cout<<endl;
//求出 A 最大的特征值
//创建初始向量 u0 for ( i = 0; i < n; ++i) u0[i] = 5; u0[0] = 1.0;
//生成矩阵 A 的非零元素 arr_ger( a_U, a_D, a_L, s, r, n, fabs(u_max));
} }
//3.此函数生成矩阵 A 的非零元素 #include "math.h" #include "iostream.h" #include "head.h" //其中,s 为上半带宽,r 为下半带宽,n 为矩阵 a 的 维数 //a_U 为 a 的上三角部分,a_D 为 a 的对角线部 分,a_L 为 a 的下三角部分 void arr_ger( double **a_U, double *a_D, double **a_L, int s, int r, int n,double u_max ) {
3. 根据幂法求出矩阵 A 模最大的特征值λ m； 4. 计算出矩阵 B = A + |λ m|I 的模最大的特征值
μ m 和模最小的特征值μ s。则有： λ 501 = μ m - |λ m| λ 1 = μ s - |λ m|
5. 求出μ k ； μ k =λ m + k(λ 501 –λ 1 )/40
return fabs( max/min ); }
//5.此函数求 A 的行列式 #include "iostream.h" #include "iomanip.h" //c 为 LU 分解后的压缩矩阵 double deta( double **c, int s, int r, int n ) {
用反幂法计算出矩阵 Bk 的模最小的特征值。 Bk = A -μ k I
(其中，k = 1,2,…,39)。 6. 求出 A 的条件数 cond(A)2 以及 A 的行列式
detA，其中： cond(A)2 = |λ m|/|λ s|； detA = det(LU) = detU。
二、全部源程序：
//1.主程序 #include "stdafx.h" #include "iostream.h" #include "iomanip.h" #include "head.h" #include "math.h"
c[i] = new double [n];
for ( i = 0; i < s; ++i ) {
a_U[i] = new double [n-1-i]; }
//求 A 的行列式 double u_max = 0, u_min = 0;//u_max 为 A 模
最大的特征值,u_min 为 A 模最小的特征值 //生成矩阵 A 的非零元素 arr_ger( a_U, a_D, a_L, s, r, n, u_max);
double b0 = 1, b1 = 1, norm; double *y0; // y = new double[n]; y0 = new double[n]; int i, j = 0; LU_decomp( c, s, r, n );
for ( i = 0; i < r; ++i ) {
a_L[i] = new double [n-i-1]; for ( j = 0; j < n-i-1; ++j ) {
a_L[i][j] = a_U[i][j]; } } }
//4.求矩阵条件数 #include "math.h" double cond( double max, double min ) {
void main( ) {
double *a_D, *u0, *y; //u0 为迭代初始值 double **a_L, **a_U, **c;
//下面的 n 代表 A 阵维数,s 代表上半带宽,r 代表下半带宽,可以随意设置
const int n = 501; const int s = 2; const int r = 2;
cout<<endl;
//求出 A 的模最大的特征值 u_max = m_max( c, a_U, a_D, a_L, s, r, n ); cout<<"A 的 模 最 大 特 征 值
为:"<<u_max<<endl; cout<<endl;
// 求出 A 的模最小的特征值 u_min = m_min( c, a_U, a_D, y, a_L, s, r, n ); cout<<"A 的 模 最 小 特 征 值
double b = 0; int i; for ( i = 0; i < n; ++i ) {
b += y[i]*u[i]; } return (b); }
//7.此函数实现反幂法求模最大的特征值 #include "head.h"
#include "math.h" #include "iostream.h" //其中,y 代表特征向量,u 代表初始非零向量,c 代 表压缩后的 a 阵 //s 代表 a 的上半带宽,r 代表 a 的下半带宽,n 代 表 a 的维数 double inv_power( double **c, double *u, double *y, int s, int r, int n ) {
特征向量
x[i] = u[i] + x[i] ;
//求出 A 最小的特征值 //创建初始向量 u0 for ( i = 0; i < n; ++i) u0[i] = 5; u0[0] = 1.0;
//生成矩阵 A 的非零元素 arr_ger( a_U, a_D, a_L, s, r, n, fabs(u_max));
//2.此函数将带状稀疏矩阵 a 压缩为 c //其中,s 为上半带宽,r 为下半带宽,n 为矩阵 a 的 维数 //a_U 为 a 的上三角部分,a_D 为 a 的对角线部 分,a_L 为 a 的下三角部分 void arr_compress( double **c, double *a_D, double **a_U, double **a_L, int s, int r, int n ) {