matlab用QR方法怎么求特征值,把程序写出来，谢谢

function l = rqrtz(A,M)

%QR算法求矩阵全部特征值

%已知矩阵: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

//数值计算程序-特征值和特征向量

//////////////////////////////////////////////////////////////

//约化对称矩阵为三对角对称矩阵

//利用Householder变换将n阶实对称矩阵约化为对称三对角矩阵

//a-长度为n*n的数组，存放n阶实对称矩阵

//n-矩阵的阶数

//q-长度为n*n的数组，返回时存放Householder变换矩阵

//b-长度为n的数组，返回时存放三对角阵的主对角线元素

//c-长度为n的数组，返回时前n-1个元素存放次对角线元素

void eastrq(double a[],int n,double q[],double b[],double c[])

//////////////////////////////////////////////////////////////

//求实对称三对角对称矩阵的全部特征值及特征向量

//利用变型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)

//////////////////////////////////////////////////////////////

//约化实矩阵为赫申伯格(Hessen berg)矩阵

//利用初等相似变换将n阶实矩阵约化为上H矩阵

//a-长度为n*n的数组，存放n阶实矩阵,返回时存放上H矩阵

//n-矩阵的阶数

void echbg(double a[],int n)

//////////////////////////////////////////////////////////////

//求赫申伯格(Hessen berg)矩阵的全部特征值

//返回值小于0表示超过迭代jt次仍未达到精度要求

//返回值大于0表示正常返回

//利用带原点位移的双重步QR方法求上H矩阵的全部特征值

//a-长度为n*n的数组，存放上H矩阵

//n-矩阵的阶数

//u-长度为n的数组，返回n个特征值的实部

//v-长度为n的数组，返回n个特征值的虚部

//eps-控制精度要求

//jt-整型变量，控制最大迭代次数

int edqr(double a[],int n,double u[],double v[],double eps,int jt)

//////////////////////////////////////////////////////////////

//求实对称矩阵的特征值及特征向量的雅格比法

//利用雅格比(Jacobi)方法求实对称矩阵的全部特征值及特征向量

//返回值小于0表示超过迭代jt次仍未达到精度要求

//返回值大于0表示正常返回

//a-长度为n*n的数组，存放实对称矩阵，返回时对角线存放n个特征值

//n-矩阵的阶数

//u-长度为n*n的数组，返回特征向量(按列存储)

//eps-控制精度要求

//jt-整型变量，控制最大迭代次数

int eejcb(double a[],int n,double v[],double eps,int jt)

//////////////////////////////////////////////////////////////

选自<<徐世良数值计算程序集(C)>>

每个程序都加上了适当地注释，陆陆续续干了几个月才整理出来的啊。

今天都给贴出来了

#include "stdio.h"

#include "math.h"

//约化对称矩阵为三对角对称矩阵

//利用Householder变换将n阶实对称矩阵约化为对称三对角矩阵

//a-长度为n*n的数组，存放n阶实对称矩阵

//n-矩阵的阶数

//q-长度为n*n的数组，返回时存放Householder变换矩阵

//b-长度为n的数组，返回时存放三对角阵的主对角线元素

//c-长度为n的数组，返回时前n-1个元素存放次对角线元素

void eastrq(double a[],int n,double q[],double b[],double c[])

{

int i,j,k,u,v

double h,f,g,h2

for (i=0i<=n-1i++)

{

for (j=0j<=n-1j++)

{

u=i*n+jq[u]=a[u]

}

}

for (i=n-1i>=1i--)

{

h=0.0

if (i>1)

{

for (k=0k<=i-1k++)

{

u=i*n+k

h=h+q[u]*q[u]

}

}

if (h+1.0==1.0)

{

c[i-1]=0.0

if (i==1)

{

c[i-1]=q[i*n+i-1]

}

b[i]=0.0

}

else

{

c[i-1]=sqrt(h)

u=i*n+i-1

if (q[u]>0.0)

{

c[i-1]=-c[i-1]

}

h=h-q[u]*c[i-1]

q[u]=q[u]-c[i-1]

f=0.0

for (j=0j<=i-1j++)

{

q[j*n+i]=q[i*n+j]/h

g=0.0

for (k=0k<=jk++)

{

g=g+q[j*n+k]*q[i*n+k]

}

if (j+1<=i-1)

{

for (k=j+1k<=i-1k++)

{

g=g+q[k*n+j]*q[i*n+k]

}

}

c[j-1]=g/h

f=f+g*q[j*n+i]

}

h2=f/(h+h)

for (j=0j<=i-1j++)

{

f=q[i*n+j]

g=c[j-1]-h2*f

c[j-1]=g

for (k=0k<=jk++)

{

u=j*n+k

q[u]=q[u]-f*c[k-1]-g*q[i*n+k]

}

}

b[i]=h

}

}

b[0]=0.0

for (i=0i<=n-1i++)

{

if ((b[i]!=0.0)&&(i-1>=0))

{

for (j=0j<=i-1j++)

{

g=0.0

for (k=0k<=i-1k++)

{

g=g+q[i*n+k]*q[k*n+j]

}

for (k=0k<=i-1k++)

{

u=k*n+j

q[u]=q[u]-g*q[k*n+i]

}

}

}

u=i*n+i

b[i]=q[u]

q[u]=1.0

if (i-1>=0)

{

for (j=0j<=i-1j++)

{

q[i*n+j]=0.0

q[j*n+i]=0.0

}

}

}

return

//求实对称三对角对称矩阵的全部特征值及特征向量

//利用变型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)

}

//约化实矩阵为赫申伯格(Hessen berg)矩阵

//利用初等相似变换将n阶实矩阵约化为上H矩阵

//a-长度为n*n的数组，存放n阶实矩阵,返回时存放上H矩阵

//n-矩阵的阶数

void echbg(double a[],int n)

{ int i,j,k,u,v

double d,t

for (k=1k<=n-2k++)

{

d=0.0

for (j=kj<=n-1j++)

{

u=j*n+k-1

t=a[u]

if (fabs(t)>fabs(d))

{

d=t

i=j

}

}

if (fabs(d)+1.0!=1.0)

{

if (i!=k)

{

for (j=k-1j<=n-1j++)

{

u=i*n+j

v=k*n+j

t=a[u]

a[u]=a[v]

a[v]=t

}

for (j=0j<=n-1j++)

{

u=j*n+i

v=j*n+k

t=a[u]

a[u]=a[v]

a[v]=t

}

}

for (i=k+1i<=n-1i++)

{

u=i*n+k-1

t=a[u]/d

a[u]=0.0

for (j=kj<=n-1j++)

{

v=i*n+j

a[v]=a[v]-t*a[k*n+j]

}

for (j=0j<=n-1j++)

{

v=j*n+k

a[v]=a[v]+t*a[j*n+i]

}

}

}

}

return

}

//求赫申伯格(Hessen berg)矩阵的全部特征值

//利用带原点位移的双重步QR方法求上H矩阵的全部特征值

//返回值小于0表示超过迭代jt次仍未达到精度要求

//返回值大于0表示正常返回

//a-长度为n*n的数组，存放上H矩阵

//n-矩阵的阶数

//u-长度为n的数组，返回n个特征值的实部

//v-长度为n的数组，返回n个特征值的虚部

//eps-控制精度要求

//jt-整型变量，控制最大迭代次数

int edqr(double a[],int n,double u[],double v[],double eps,int jt)

{

int m,it,i,j,k,l,ii,jj,kk,ll

double b,c,w,g,xy,p,q,r,x,s,e,f,z,y

it=0

m=n

while (m!=0)

{

l=m-1

while ((l>0)&&(fabs(a[l*n+l-1])>eps*(fabs(a[(l-1)*n+l-1])+fabs(a[l*n+l]))))

{

l=l-1

}

ii=(m-1)*n+m-1

jj=(m-1)*n+m-2

kk=(m-2)*n+m-1

ll=(m-2)*n+m-2

if (l==m-1)

{

u[m-1]=a[(m-1)*n+m-1]

v[m-1]=0.0

m=m-1it=0

}

else if (l==m-2)

{

b=-(a[ii]+a[ll])

c=a[ii]*a[ll]-a[jj]*a[kk]

w=b*b-4.0*c

y=sqrt(fabs(w))

if (w>0.0)

{

xy=1.0

if (b<0.0)

{

xy=-1.0

}

u[m-1]=(-b-xy*y)/2.0

u[m-2]=c/u[m-1]

v[m-1]=0.0v[m-2]=0.0

}

else

{

u[m-1]=-b/2.0

u[m-2]=u[m-1]

v[m-1]=y/2.0

v[m-2]=-v[m-1]

}

m=m-2

it=0

}

else

{

if (it>=jt)

{

printf("fail\n")

return(-1)

}

it=it+1

for (j=l+2j<=m-1j++)

{

a[j*n+j-2]=0.0

}

for (j=l+3j<=m-1j++)

{

a[j*n+j-3]=0.0

}

for (k=lk<=m-2k++)

{

if (k!=l)

{

p=a[k*n+k-1]

q=a[(k+1)*n+k-1]

r=0.0

if (k!=m-2)

{

r=a[(k+2)*n+k-1]

}

}

else

{

x=a[ii]+a[ll]

y=a[ll]*a[ii]-a[kk]*a[jj]

ii=l*n+l

jj=l*n+l+1

kk=(l+1)*n+l

ll=(l+1)*n+l+1

p=a[ii]*(a[ii]-x)+a[jj]*a[kk]+y

q=a[kk]*(a[ii]+a[ll]-x)

r=a[kk]*a[(l+2)*n+l+1]

}

if ((fabs(p)+fabs(q)+fabs(r))!=0.0)

{

xy=1.0

if (p<0.0)

{

xy=-1.0

}

s=xy*sqrt(p*p+q*q+r*r)

if (k!=l)

{

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<=m-1j++)

{

ii=k*n+j

jj=(k+1)*n+j

p=x*a[ii]+e*a[jj]

q=e*a[ii]+y*a[jj]

r=f*a[ii]+g*a[jj]

if (k!=m-2)

{

kk=(k+2)*n+j

p=p+f*a[kk]

q=q+g*a[kk]

r=r+z*a[kk]

a[kk]=r

}

a[jj]=q

a[ii]=p

}

j=k+3

if (j>=m-1)

{

j=m-1

}

for (i=li<=ji++)

{

ii=i*n+k

jj=i*n+k+1

p=x*a[ii]+e*a[jj]

q=e*a[ii]+y*a[jj]

r=f*a[ii]+g*a[jj]

if (k!=m-2)

{

kk=i*n+k+2

p=p+f*a[kk]

q=q+g*a[kk]

r=r+z*a[kk]

a[kk]=r

}

a[jj]=q

a[ii]=p

}

}

}

}

}

return(1)

}

//求实对称矩阵的特征值及特征向量的雅格比法

//利用雅格比(Jacobi)方法求实对称矩阵的全部特征值及特征向量

//返回值小于0表示超过迭代jt次仍未达到精度要求

//返回值大于0表示正常返回

//a-长度为n*n的数组，存放实对称矩阵，返回时对角线存放n个特征值

//n-矩阵的阶数

//u-长度为n*n的数组，返回特征向量(按列存储)

//eps-控制精度要求

//jt-整型变量，控制最大迭代次数

int eejcb(double a[],int n,double v[],double eps,int jt)

{

int i,j,p,q,u,w,t,s,l

double fm,cn,sn,omega,x,y,d

l=1

for (i=0i<=n-1i++)

{

v[i*n+i]=1.0

for (j=0j<=n-1j++)

{

if (i!=j)

{

v[i*n+j]=0.0

}

}

}

while (1==1)

{

fm=0.0

for (i=0i<=n-1i++)

{

for (j=0j<=n-1j++)

{

d=fabs(a[i*n+j])

if ((i!=j)&&(d>fm))

{

fm=d

p=i

q=j

}

}

}

if (fm<eps)

{

return(1)

}

if (l>jt)

{

return(-1)

}

l=l+1

u=p*n+q

w=p*n+p

t=q*n+p

s=q*n+q

x=-a[u]

y=(a[s]-a[w])/2.0

omega=x/sqrt(x*x+y*y)

if (y<0.0)

{

omega=-omega

}

sn=1.0+sqrt(1.0-omega*omega)

sn=omega/sqrt(2.0*sn)

cn=sqrt(1.0-sn*sn)

fm=a[w]

a[w]=fm*cn*cn+a[s]*sn*sn+a[u]*omega

a[s]=fm*sn*sn+a[s]*cn*cn-a[u]*omega

a[u]=0.0

a[t]=0.0

for (j=0j<=n-1j++)

{

if ((j!=p)&&(j!=q))

{

u=p*n+j

w=q*n+j

fm=a[u]

a[u]=fm*cn+a[w]*sn

a[w]=-fm*sn+a[w]*cn

}

}

for (i=0i<=n-1i++)

{

if ((i!=p)&&(i!=q))

{

u=i*n+p

w=i*n+q

fm=a[u]

a[u]=fm*cn+a[w]*sn

a[w]=-fm*sn+a[w]*cn

}

}

for (i=0i<=n-1i++)

{

u=i*n+p

w=i*n+q

fm=v[u]

v[u]=fm*cn+v[w]*sn

v[w]=-fm*sn+v[w]*cn

}

}

return(1)

}

#include <stdio.h>

void main(){

int i,j,k=1,n,s=0

scanf("%d",&n)

for(i=1ni++){

j=n%10n/=10

if(j%2&&i%2||j%2==0&&i%2==0)

s+=kk*=2

}

printf("%d",s)

}

//运行示例：

---