
//利用变型QR方法计算实对称三对角矩阵全部特征值及特衡郑征向量
//n-矩阵的阶数
//b-长度为n的数组,返回时存放三对角阵的主对角线元素
//c-长度为n的数握饥组,返回时前n-1个元素存放次对角线咐皮颂元素
//q-长度为n*n的数组,若存放单位矩阵,则返回实对称三对角矩阵的特征向量组
// 若存放Householder变换矩阵,则返回实对称矩阵A的特征向量组
//a-长度为n*n的数组,存放n阶实对称矩阵
int ebstq(int n,double b[],double c[],double q[],double eps,int l)
{
int i,j,k,m,it,u,v
double d,f,h,g,p,r,e,s
c[n-1]=0.0
d=0.0
f=0.0
for (j=0j<=n-1j++)
{
it=0
h=eps*(fabs(b[j])+fabs(c[j]))
if (h>d)
{
d=h
}
m=j
while ((m<=n-1)&&(fabs(c[m])>d))
{
m=m+1
}
if (m!=j)
{
do
{
if (it==l)
{
printf("fail\n")
return(-1)
}
it=it+1
g=b[j]
p=(b[j+1]-g)/(2.0*c[j])
r=sqrt(p*p+1.0)
if (p>=0.0)
{
b[j]=c[j]/(p+r)
}
else
{
b[j]=c[j]/(p-r)
}
h=g-b[j]
for (i=j+1i<=n-1i++)
{
b[i]=b[i]-h
}
f=f+h
p=b[m]
e=1.0
s=0.0
for (i=m-1i>=ji--)
{
g=e*c[i]
h=e*p
if (fabs(p)>=fabs(c[i]))
{
e=c[i]/p
r=sqrt(e*e+1.0)
c[i+1]=s*p*r
s=e/r
e=1.0/r
}
else
{
e=p/c[i]
r=sqrt(e*e+1.0)
c[i+1]=s*c[i]*r
s=1.0/r
e=e/r
}
p=e*b[i]-s*g
b[i+1]=h+s*(e*g+s*b[i])
for (k=0k<=n-1k++)
{
u=k*n+i+1
v=u-1
h=q[u]
q[u]=s*q[v]+e*h
q[v]=e*q[v]-s*h
}
}
c[j]=s*p
b[j]=e*p
}
while (fabs(c[j])>d)
}
b[j]=b[j]+f
}
for (i=0i<=n-1i++)
{
k=ip=b[i]
if (i+1<=n-1)
{
j=i+1
while ((j<=n-1)&&(b[j]<=p))
{
k=j
p=b[j]
j=j+1
}
}
if (k!=i)
{
b[k]=b[i]
b[i]=p
for (j=0j<=n-1j++)
{
u=j*n+i
v=j*n+k
p=q[u]
q[u]=q[v]
q[v]=p
}
}
}
return(1)
}
http://zhidao.baidu.com/question/42993521.html?si=2
估计没人会免费给带宽虚你提供QR算法代码的,我自己写的VB.net的QR代码都费了不少巧脊劲。这些代码一般都是公司开发的,都商蠢燃品化了,你需要的话可以购买,或者那什么什么的,你懂的。
function l = rqrtz(A,M)%瑞利商位移的QR算法求矩阵全部特征值
%已知矩阵:A
%迭代步数:M
%求得的矩阵特征值:lA = 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)
end4.4 QR算 法QR算法也是一种迭代算法,是目前计算任意实的非奇异矩阵全部特征值问题的最有效的方法之一.该方法的基础是构造矩阵序列 ,并对它进行QR分解.由线性代数知识知道,若A为非奇异方阵,则A可迅洞以分解为正交矩阵Q与上三角形矩阵R的乘积,即A=QR,而且当R的对角线元素符号取定时,分解式是唯一的.若A为奇异方阵,则零为A的特征值.任取一数p不是A的特征值,则A-pI为非奇异方阵.只要求出A-pI的特征值,就很容易求出A的特征值,所以假设A为非奇异方阵,并不妨碍讨论的一般性.设A为非奇异方阵,令 ,对 进行QR分解,即把 分解为正交矩阵 与上三角形矩阵 的乘积= 做矩阵 继续对 进行QR分解 并定义 一般地,递推公式为QR算法就是利用矩阵的QR分解,按上述递推公式构造矩阵序列 .只要A为非奇异方阵,则由QR算法就完全确定 .这个矩阵序列 具有下列性质.性质1 所有 都相似,它们具有相同的特征值.证明 因为 若令 ,则 为正交阵,且有 因此 与A相似,它们具有相同的特征值.性质2 的QR分解式为其中 证明 用归纳法.显然当k=1时,有 假设 有分解式 于是因为 ,所以 因为 都是正交阵,所以 也是正交阵,同样 也是上三角形阵,从而 的QR分解式为 由前面的讨论知 .这说明QR算法的收敛性有正交矩阵序列 的性质决定.定理1 如果 收敛于非奇异矩阵 为上三角形矩阵,则 存在并且是上三角形矩阵.证明 因为 收敛,故下面极限存在 由于 为上三角形矩阵,所以 为上三角形矩阵.又亩者枯因为 所以 存在,并且是上三角形矩阵.定理2 (QR算法的收敛性)设A为n 阶实矩阵,且1) A的特征值满足: 2) ,其中 且设 有三角分解式 =LU(L为单位下三角阵,U为上三角阵),则由QR算法得到的矩阵序列 本质上收敛于上三角形矩阵.即 满足 当 当 的极限不一定存在证明 因为 ,矩阵 决定 的收敛性.又 我们利用 求 ,然后讨论 的收敛性.由定理条件 得 令 其中 的(i,j)元素 为 于是 由假设,当i>j时, 故设方阵X的QR分解式为 由由 知,对充分大的 非奇异,它应有唯一的QR分解式 ,并且 于是但上三角阵 的对角线元素不一定大于零.为此,引入对角矩阵以便保证( )的对角线元素都是正数,从而得到 的QR分解式 由 的QR分解式的唯一性得到 从而由于 ,所以 从而 其中 于是 因为 为上三角阵, 为对角阵,且嫌物元素为1或-1,所以 当当 的极限不一定存在例 用QR算法求矩阵 的特征值.A的特征值为-1,4,1+2i,1-2i.解 令 ,用施密特正交化过程将 分解为 将 与 逆序相乘,求出 用 代替A重复上面过程,计算11次得 由 不难看出,矩阵A的一个特征值是4,另一个特征值是-1,其他两个特征值是方程 的根.求得为
欢迎分享,转载请注明来源:内存溢出
微信扫一扫
支付宝扫一扫
评论列表(0条)