%已知矩阵:A
%迭代步数:M
%求得的矩阵特征值:l
A = hess(A)
for i=1:M
N = size(A)
n = N(1,1)
u = A(n,n)
[q,r]=qr(A-u*eye(n,n))
A = r*q+u*eye(n,n)
l = diag(A)
end
------------------------------------
A=[0 5 0 0 0 01 0 4 0 0 00 1 0 3 0 00 0 1 0 2 00 0 0 1 0 10 0 0 0 1 0]
A =
0 5 0 0 0 0
1 0 4 0 0 0
0 1 0 3 0 0
0 0 1 0 2 0
0 0 0 1 0 1
0 0 0 0 1 0
>>rqrtz(A,50)
ans =
-3.2030
3.2030
-1.8837
1.8837
-0.6167
0.6167
>>eig(A)
ans =
-3.3243
3.3243
-1.8892
-0.6167
1.8892
0.6167
c++求矩阵特征值方法：函数：
Matrix_EigenValue(double *K1,int n,int LoopNumber,double Error1,double *Ret)
K1：n阶方阵
n：方阵K1的阶数
LoopNumber：在误差无法保证能得到结果时运算的最大次数
Error1：误差控制变量
Ret：返回的一个n*2的矩阵。矩阵的每一行是求得的特征值，第一列代表特征值实数部分，第二列代表虚数部分
函数成功返回True，失败返回False
特别说明：
Matrix_Hessenberg：把n阶方阵K1化为上三角Hessenberg矩阵,其中A储存上三角Hessenberg矩阵
源代码:
bool Matrix_EigenValue(double *K1,int n,int LoopNumber,double Error1,double *Ret)
{
int i,j,k,t,m,Loop1
double b,c,d,g,xy,p,q,r,x,s,e,f,z,y,temp,*A
A=new double[n*n]
Matrix_Hessenberg(K1,n,A)
m=n
Loop1=LoopNumber
while(m!=0)
{
t=m-1
while(t>0)
{
temp=abs(A[(t-1)*n+t-1])
temp+=abs(A[t*n+t])
temp=temp*Error1
if (abs(A[t*n+t-1])>temp)
{
t--
}
else
{
break
}
}
if (t==m-1)
{
Ret[(m-1)*2]=A[(m-1)*n+m-1]
Ret[(m-1)*2+1]=0
m-=1
Loop1=LoopNumber
}
else if(t==m-2)
{
b=-A[(m-1)*n+m-1]-A[(m-2)*n+m-2]
c=A[(m-1)*n+m-1]*A[(m-2)*n+m-2]-A[(m-1)*n+m-2]*A[(m-2)*n+m-1]
d=b*b-4*c
y=sqrt(abs(d))
if (d>0)
{
xy=1
if (b<0)
{
xy=-1
}
Ret[(m-1)*2]=-(b+xy*y)/2
Ret[(m-1)*2+1]=0
Ret[(m-2)*2]=c/Ret[(m-1)*2]
Ret[(m-2)*2+1]=0
}
else
{
Ret[(m-1)*2]=-b/2
Ret[(m-2)*2]=-b/2
Ret[(m-1)*2+1]=y/2
Ret[(m-2)*2+1]=-y/2
}
m-=2
Loop1=LoopNumber
}
else
{
if (Loop1<1)
{
return false
}
Loop1--
j=t+2
while (j<m)
{
A[j*n+j-2]=0
j++
}
j=t+3
while (j<m)
{
A[j*n+j-3]=0
j++
}
k=t
while (k<m-1)
{
if (k!=t)
{
p=A[k*n+k-1]
q=A[(k+1)*n+k-1]
if (k!=m-2)
{
r=A[(k+2)*n+k-1]
}
else
{
r=0
}
}
else
{
b=A[(m-1)*n+m-1]
c=A[(m-2)*n+m-2]
x=b+c
y=b*c-A[(m-2)*n+m-1]*A[(m-1)*n+m-2]
p=A[t*n+t]*(A[t*n+t]-x)+A[t*n+t+1]*A[(t+1)*n+t]+y
q=A[(t+1)*n+t]*(A[t*n+t]+A[(t+1)*n+t+1]-x)
r=A[(t+1)*n+t]*A[(t+2)*n+t+1]
}
if (p!=0 || q!=0 || r!=0)
{
if (p<0)
{
xy=-1
}
else
{
xy=1
}
s=xy*sqrt(p*p+q*q+r*r)
if (k!=t)
{
A[k*n+k-1]=-s
}
e=-q/s
f=-r/s
x=-p/s
y=-x-f*r/(p+s)
g=e*r/(p+s)
z=-x-e*q/(p+s)
for (j=kj<mj++)
{
b=A[k*n+j]
c=A[(k+1)*n+j]
p=x*b+e*c
q=e*b+y*c
r=f*b+g*c
if (k!=m-2)
{
b=A[(k+2)*n+j]
p+=f*b
q+=g*b
r+=z*b
A[(k+2)*n+j]=r
}
A[(k+1)*n+j]=q
A[k*n+j]=p
}
j=k+3
if (j>m-2)
{
j=m-1
}
for (i=ti<j+1i++)
{
b=A[i*n+k]
c=A[i*n+k+1]
p=x*b+e*c
q=e*b+y*c
r=f*b+g*c
if (k!=m-2)
{
b=A[i*n+k+2]
p+=f*b
q+=g*b
r+=z*b
A[i*n+k+2]=r
}
A[i*n+k+1]=q
A[i*n+k]=p
}
}
k++
}
}
}
delete []A
return true
}
#include<iostream.h>#include<stdlib.h>
#include<math.h>
#include<iomanip.h>
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
//class Matrix定义矩阵类
const int Max_xy = 20//矩阵的最大维数
class Matrix
{
private:
double data[Max_xy][Max_xy]
unsigned x, y //x,y
public:
Matrix() //默认构造函数
Matrix(const Matrix &source) //拷贝构造函数
void creat() //输入矩阵
void init()
void transpose() //矩阵转置
void show() //输入此矩阵
double mode() const //求一维矩阵的长度
void check_shiduichen() //检查是否为是对称矩阵
void creat_unit(unsigned i) //生成i行单位矩阵
void set_x(unsigned xx) //设置行数
void set_y(unsigned yy) //设置列数
unsigned get_x() //得到行数
unsigned get_y() //得到列数
void shucheng(double changshu) //数乘运算
void setdata(unsigned i, unsigned j, double source)//定位输入数据
double getdata(unsigned i, unsigned j) //定位得到数据
void sturm() //求特征值
unsigned sturm_s(double m) //计算sturm系列的同好数
Matrix operator = (const Matrix &right)
friend Matrix &operator + (const Matrix &left, const Matrix &right) //重载+号
friend Matrix &operator - (const Matrix &left, const Matrix &right) //重载-号
friend Matrix &operator * (const Matrix &left, const Matrix &right) //重载乘号
friend ostream &operator <<(ostream &os, const Matrix &source) //重载输出
friend void Householder(Matrix &source) //用Householde矩阵将实对称矩阵化为三对角矩阵
}
Matrix temp_Matrix //全局变量Matrix
//===================================================================
//--------------------默认构造函数
Matrix::Matrix()
{
init()
}
//----------------------------拷贝构造函数
Matrix::Matrix(const Matrix &source)
{
init()
x = source.x
y = source.y
for(unsigned i = 0i <xi++)
for(unsigned j = 0j <yj++)
data[i][j] = source.data[i][j]
}
//------------------------------------------初始化矩阵元素
void Matrix::init()
{
x = y = 0
for(unsigned i = 0i <Max_xyi++)
for(unsigned j = 0j <Max_xyj++)
data[i][j] = 0
}
//------------------------------矩阵转置
void Matrix::transpose()
{
double temp
int t
for(unsigned i = 0i <Max_xyi++)
for(unsigned j = 0j <= ij++)
{
temp = data[i][j]
data[i][j] = data[j][i]
data[j][i] = temp
}
t = x
x = y
y = t
}
//--------------------------------------求一维矩阵的长度
double Matrix::mode() const
{
double s = 0
unsigned i, j
if(x == 1)
for(i = 0, j = 0j <yj++)
s += data[i][j] * data[i][j]
else if(y == 1)
for(i = 0, j = 0i <xi++)
s += data[i][j] * data[i][j]
else
{
cout <<"\n不是一维的!"
exit(0)
}
s = sqrt(s)
return (s)
}
//----------------------------------------重载=号
Matrix Matrix::operator = (const Matrix &source)
{
x = source.x
y = source.y
for(unsigned i = 0i <xi++)
for(unsigned j = 0j <yj++)
data[i][j] = source.data[i][j]
return *this
}
//-------------------------------------------重载+号
Matrix &operator + (const Matrix &left, const Matrix &right)
{
if(left.x != right.x || left.y != right.y)
{
cout <<"\n维数不相等,不能相加!"
exit(0)
}
for(unsigned i = 0i <left.xi++)
for(unsigned j = 0j <left.yj++)
temp_Matrix.data[i][j] = right.data[i][j] + left.data[i][j]
temp_Matrix.x = right.x
temp_Matrix.y = right.y
return temp_Matrix
}
//---------------------------------------------重载+号
Matrix &operator - (const Matrix &left, const Matrix &right)
{
if(left.x != right.x || left.y != right.y)
{
cout <<"\n维数不相等,不能相减!"
exit(0)
}
for(unsigned i = 0i <left.xi++)
for(unsigned j = 0j <left.yj++)
temp_Matrix.data[i][j] = left.data[i][j] - right.data[i][j]
temp_Matrix.x = right.x
temp_Matrix.y = right.y
return temp_Matrix
}
//----------------------------------------重载乘号
Matrix &operator * (const Matrix &left, const Matrix &right)
{
if(left.y != right.x)
{
cout <<"\n两个矩阵相乘错误."
exit(0)
}
temp_Matrix.init()
unsigned i, j, k
for(i = 0i <left.xi++)
for(j = 0j <right.yj++)
for(k = 0k <left.yk++)
temp_Matrix.data[i][j] += left.data[i][k] * right.data[k][j]
temp_Matrix.x = left.x
temp_Matrix.y = right.y
return temp_Matrix
}
//-------------------------------------------输入矩阵
void Matrix::creat()
{
cout <<"输如行列式:"
cout <<"\n行数:"
cin >>x
cout <<"列数:"
cin >>y
for(unsigned i = 0i <xi++)
{
cout <<"输入第" <<i + 1 <<"行:"
for(unsigned j = 0j <yj++)
cin >>data[i][j]
}
}
//----------------------------------------------输出矩阵
void Matrix::show()
{
unsigned i, j
cout <<"\n\n矩阵表示如下:"
for(i = 0i <xi++)
{
cout <<endl
for(j = 0j <yj++)
cout <<setw(7) <<setiosflags(ios::left) <<data[i][j]
}
}
//----------------------------------用Householder矩阵化为实对称矩阵
void Householder(Matrix &source)
{
unsigned i, lenth, k, m, n, flag
double temp_lie_x[Max_xy], temp_lie_y[Max_xy], temp[Max_xy]
double s
for(i = 0i <source.x - 2i++)
{
for(k = 0k <Max_xyk++)
{
//初始化为0
temp_lie_x[k] = 0
temp_lie_y[k] = 0
temp[k] = 0
}
for(lenth = 0lenth + i + 1 <source.xlenth++) //提取第data[i+1][i]到data[x-1][i]的数据存到temp_lie[Max_xy]
temp_lie_x[lenth] = source.data[lenth+i+1][i]
for(k = 0, s = 0k <lenthk++) //j为临时变量的个数.
s += temp_lie_x[k] * temp_lie_x[k]
s = sqrt(s)
temp_lie_y[0] = -s
for(k = 1k <lenthk++)
temp_lie_y[k] = 0
for(k = 0k <lenthk++)
temp[k] = temp_lie_x[k] - temp_lie_y[k]
for(k = 0, flag = 0k <lenthk++) //假如以上两个向量相等则退出,则跳出以进行下一次变换
if(temp[k] != 0)
{
flag = 1
break
}
if(flag == 0)
continue
Matrix part_h, I, x_y//定义Matrix变量
I.creat_unit(lenth)
x_y.x = lenth, x_y.y = 1
//对x_y赋值
for(k = 0k <lenthk++)
x_y.data[k][0] = temp[k]
//求x_y的转置
Matrix zhuanzhi_x_y( x_y ) //拷贝构造函数
zhuanzhi_x_y.transpose()
s = 2.0 / ( x_y.mode() * x_y.mode() )
x_y.shucheng(s)
temp_Matrix = x_y * zhuanzhi_x_y
part_h = I - x_y * zhuanzhi_x_y
Matrix H
H.creat_unit(source.x)
for(m = i + 1m <source.xm++)
for(n = i + 1n <source.yn++)
H.data[m][n] = part_h.data[m-i-1][n-i-1]//得到最后的Householder矩阵
source = source * H
source = H * source
}
for(i = 0i <source.xi++)
for(k = 0k <source.yk++)
if(fabs(source.data[i][k]) <1e-13)
source.data[i][k] = 0
}
//------------------------------------检查是否为实对称矩阵
void Matrix::check_shiduichen()
{
if(x != y)
{
cout <<"\n\n不是是对称矩阵(行列不相等)\n\n"
exit(0)
}
for(unsigned i = 0i <xi++)
for(unsigned j = 0j <yj++)
{
if(data[i][j] != data[j][i])
{
cout <<"\n\n不是实对称矩阵!(不对称!)\n\n"
exit(0)
}
}
}
ostream &operator <<(ostream &os, const Matrix &source)
{
unsigned i, j
for(i = 0i <source.xi++)
{
os <<"\n"
for(j = 0j <source.yj++)
os <<setw(10) <<setiosflags(ios::left) <<source.data[i][j] <<"\t"
}
os <<endl
return os
}
void Matrix::set_x(unsigned xx)//设置行数
{
x = xx
}
void Matrix::set_y(unsigned yy)//设置列数
{
y = yy
}
unsigned Matrix::get_x() //得到行数
{
return x
}
unsigned Matrix::get_y() //得到列数
{
return y
}
double Matrix::getdata(unsigned i, unsigned j)//定位得到数据
{
return data[i][j]
}
void Matrix::setdata(unsigned i, unsigned j, double source)
{
data[i][j] = source
if(x <i) x = i
if(y <j) y = j
}
void Matrix::creat_unit(unsigned lenth) //生成i行单位矩阵
{
init()
x = y = lenth
for(unsigned i = 0i <lenthi++)
data[i][i] = 1
}
void Matrix::shucheng(double changshu) //数乘运算
{
unsigned i, j
for(i = 0i <xi++)
for(j = 0j <yj++)
data[i][j] *= changshu
}
void Matrix::sturm()
{
double s, r, maxhang, a, b
unsigned i, j, m
for(i = 0, maxhang = 0i <xi++)
{
for( s = 0, j = 0j <yj++)
s += fabs(data[i][j])
if(s >maxhang)
maxhang = s
}
a = -maxhang
b = maxhang
m = sturm_s(a) - sturm_s(b)
for(i = 1i <= mi++)
{
a = -maxhang
b = maxhang
do
{
r = 0.5 * (a + b)
if(sturm_s(r) >= i)
a = r
else
b = r
}
while(fabs(a - b) >1e-11)
cout <<"\n特征值" <<i <<": " <<setiosflags(ios::left) <<setw(10) <<setprecision(10)
<<setiosflags(ios::fixed) <<0.5 * (a + b)
}
}
unsigned Matrix::sturm_s(double r)
{
double p[Max_xy+1]
int temp[Max_xy+1]
unsigned m, i
p[0] = 1
p[1] = data[0][0] - r
for(i = 2i <= xi++)
p[i] = (data[i-1][i-1] - r) * p[i-1] - data[i-2][i-1] * data[i-2][i-1] * p[i-2]
temp[0] = 1
for(i = 1i <= xi++)
if(p[i] >1e-14)
temp[i] = 1
else if(p[i] <-1e-14)
temp[i] = -1
else
temp[i] = temp[i-1]
for(i = 1, m = 0i <= xi++)
if(temp[i] + temp[i-1] != 0)
m++
return m
}
void main()
{
Matrix a
a.creat()
cout <<"输入的矩阵为:"
cout <<a
a.sturm()
cout <<endl
}
欢迎分享，转载请注明来源：内存溢出
微信扫一扫
支付宝扫一扫
评论列表（0条）