对称三对角矩阵特征值的二分法

数学软件实验任务书

课程名称数学软件实验班级数0901

实验课题对称三对角矩阵特征值的二分法

实验目的熟悉对称三对角矩阵特征值的二分法

运用Matlab/C/C++/Java/Maple/Mathematica等其中实验要求

一种语言完成

对称三对角矩阵特征值的二分法

实验内容

成绩教师

实验1对称三对角矩阵特征值的二分法

1 实验原理

对于对称三对角矩阵：

11122

2111

n n n n n c b b c b C b c b b c ----⎛⎫ ⎪ ⎪ ⎪= ⎪ ⎪ ⎪⎝⎭ 设0i b ≠，1,2,,1i n =- ，记特征矩阵I C λ-的左上角的k 阶子式为()k p λ，设0()1p λ=，利用行列式的展开式，可得()k p λ的递推公式：

0112112()1() 1,2,,()()()

k k k k k p p a k n p a p b p λλλλλλ---=⎧⎪=-⎪=⎨⎪⎪=--⎩ ()det()n p I C λλ=-为C 的特征多项式，n 个零点为矩阵C 的n 个特征值 2 实验程序

程序：

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=0; j<=n-1; j++)

{

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+1; i<=n-1; i++) {

b=b-h;

}

f=f+h;

p=b[m];

e=1.0;

s=0.0;

for (i=m-1; i>=j; i--)

{

g=e*c;

h=e*p;

if (fabs(p)>=fabs(c)) {

e=c/p;

r=sqrt(e*e+1.0);

c[i+1]=s*p*r;

s=e/r;

e=1.0/r;

}

else

{

e=p/c;

r=sqrt(e*e+1.0);

c[i+1]=s*c*r;

s=1.0/r;

e=e/r;

}

p=e*b-s*g;

b[i+1]=h+s*(e*g+s*b); for (k=0; k<=n-1; k++) {

u=k*n+i+1;

v=u-1;

h=q;

q=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=0; i<=n-1; i++)

{

k=i; p=b;

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;

b=p;

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

{

u=j*n+i;

v=j*n+k;

p=q; q=q[v];

q[v]=p;

}

}

}

return(1);

}