❶ 如何用c語言編寫求對稱矩陣的特徵值和特徵向量的程序
//數值計算程序-特徵值和特徵向量
//////////////////////////////////////////////////////////////
//約化對稱矩陣為三對角對稱矩陣
//利用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=0; i<=n-1; i++)
{
for (j=0; j<=n-1; j++)
{
u=i*n+j; q[u]=a[u];
}
}
for (i=n-1; i>=1; i--)
{
h=0.0;
if (i>1)
{
for (k=0; k<=i-1; k++)
{
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=0; j<=i-1; j++)
{
q[j*n+i]=q[i*n+j]/h;
g=0.0;
for (k=0; k<=j; k++)
{
g=g+q[j*n+k]*q[i*n+k];
}
if (j+1<=i-1)
{
for (k=j+1; k<=i-1; k++)
{
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=0; j<=i-1; j++)
{
f=q[i*n+j];
g=c[j-1]-h2*f;
c[j-1]=g;
for (k=0; k<=j; k++)
{
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=0; i<=n-1; i++)
{
if ((b[i]!=0.0)&&(i-1>=0))
{
for (j=0; j<=i-1; j++)
{
g=0.0;
for (k=0; k<=i-1; k++)
{
g=g+q[i*n+k]*q[k*n+j];
}
for (k=0; k<=i-1; k++)
{
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=0; j<=i-1; j++)
{
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=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[i]=b[i]-h;
}
f=f+h;
p=b[m];
e=1.0;
s=0.0;
for (i=m-1; i>=j; i--)
{
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=0; k<=n-1; k++)
{
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=0; i<=n-1; i++)
{
k=i; p=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=0; j<=n-1; j++)
{
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=1; k<=n-2; k++)
{
d=0.0;
for (j=k; j<=n-1; j++)
{
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-1; j<=n-1; j++)
{
u=i*n+j;
v=k*n+j;
t=a[u];
a[u]=a[v];
a[v]=t;
}
for (j=0; j<=n-1; j++)
{
u=j*n+i;
v=j*n+k;
t=a[u];
a[u]=a[v];
a[v]=t;
}
}
for (i=k+1; i<=n-1; i++)
{
u=i*n+k-1;
t=a[u]/d;
a[u]=0.0;
for (j=k; j<=n-1; j++)
{
v=i*n+j;
a[v]=a[v]-t*a[k*n+j];
}
for (j=0; j<=n-1; j++)
{
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-1; it=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.0; v[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+2; j<=m-1; j++)
{
a[j*n+j-2]=0.0;
}
for (j=l+3; j<=m-1; j++)
{
a[j*n+j-3]=0.0;
}
for (k=l; k<=m-2; k++)
{
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=k; j<=m-1; j++)
{
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=l; i<=j; i++)
{
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=0; i<=n-1; i++)
{
v[i*n+i]=1.0;
for (j=0; j<=n-1; j++)
{
if (i!=j)
{
v[i*n+j]=0.0;
}
}
}
while (1==1)
{
fm=0.0;
for (i=0; i<=n-1; i++)
{
for (j=0; j<=n-1; j++)
{
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=0; j<=n-1; j++)
{
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=0; i<=n-1; i++)
{
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=0; i<=n-1; i++)
{
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);
}
❷ 四階矩陣,所有元素都是1,要怎麼算特徵值,求簡單點的方法
|A|=0,則它必有特徵值0,又因為r(A)=1,AX=0的解空間的維數是4-r(A)=3,從而0是A的三重特徵值,由於A的各行加起來都是4,則設X0=(1,1,1,1)^T,便有AX0=4X0,從而4也是A的特徵值,故A的全部特徵值0,0,0,4。
判斷矩陣可對角化的充要條件:
矩陣可對角化有兩個充要條件:
1、矩陣有n個不同的特徵向量。
2、特徵向量重根的重數等於基礎解系的個數。對於第二個充要條件,則需要出現二重以上的重特徵值可驗證(一重相當於沒有重根)。
(2)c語言求4階特徵值與特徵向量擴展閱讀:
求n階矩陣A的特徵值的基本方法:
根據定義可改寫為關系式,為單位矩陣(其形式為主對角線元素為λ-,其餘元素乘以-1)。要求向量具有非零解,即求齊次線性方程組有非零解的值。即要求行列式。
解此行列式獲得的值即為矩陣A的特徵值。將此值回代入原式求得相應的,即為輸入這個行列式的特徵向量。
❸ 如何用c語言寫求矩陣的特徵值和特徵向量
方法1:推導出det(aA-I)=0的解析式,這應該是個四次方程,因為只有4階,不是很困難的,寫出後就可以用方程求根的方法求解(如newton迭代法)
方法2:如果你是對角優勢陣,也就是對角線上的值的絕對值,比同行所有其他元素的絕對值的和還大,可以通過局部旋轉的方法把矩陣「能量」集中到對角線
這個是方法,你可以自己去寫一下試試~
❹ 特徵值和特徵向量怎麼求
對於特徵值λ和特徵向量a,得到Aa=aλ
於是把每個特徵值和特徵向量寫在一起
注意對於實對稱矩陣不同特徵值的特徵向量一定正交
得到矩陣P,再求出其逆矩陣P^(-1)
可以解得原矩陣A=PλP^(-1)
設A為n階矩陣,若存在常數λ及n維非零向量x,使得Ax=λx,則稱λ是矩陣A的特徵值,x是A屬於特徵值λ的特徵向量。
一個矩陣A的特徵值可以通過求解方程pA(λ) = 0來得到。 若A是一個n×n矩陣,則pA為n次多項式,因而A最多有n個特徵值。
反過來,代數基本定理說這個方程剛好有n個根,如果重根也計算在內的話。所有奇數次的多項式必有一個實數根,因此對於奇數n,每個實矩陣至少有一個實特徵值。在實矩陣的情形,對於偶數或奇數的n,非實數特徵值成共軛對出現。
(4)c語言求4階特徵值與特徵向量擴展閱讀
求矩陣的全部特徵值和特徵向量的方法如下:
第一步:計算的特徵多項式;
第二步:求出特徵方程的全部根,即為的全部特徵值;
第三步:對於的每一個特徵值,求出齊次線性方程組。
若是的屬於的特徵向量,則也是對應於的特徵向量,因而特徵向量不能由特徵值惟一確定.反之,不同特徵值對應的特徵向量不會相等,亦即一個特徵向量只能屬於一個特徵值。
在A變換的作用下,向量ξ僅僅在尺度上變為原來的λ倍。稱ξ是A 的一個特徵向量,λ是對應的特徵值(本徵值),是(實驗中)能測得出來的量,與之對應在量子力學理論中,很多量並不能得以測量,當然,其他理論領域也有這一現象。
❺ C語言求特徵向量
特徵值及特徵向量 //利用變型QR方法計算實對稱三對角矩陣全部特徵值及特徵向量 //n-矩陣的階數 //b- 特徵值及特徵向量 //利用變型QR方法計算實對稱三對角矩陣全部特徵值及特徵向量 //n-矩陣的階數 //b-
❻ (在線等!)求特徵值和特徵向量的步驟是
令|A-λE|=0,求出λ值。A是n階矩陣,Ax=λx,則x為特徵向量,λ為特徵值。
設矩陣為A,特徵向量是t,特徵值是x,At=x*t,移項得(A-x*I)t=0,
∵t不是零向量
∴A-x*I=0,(2-x)(1-x)(-x)-4(2-x)=0,化簡得(x-2)(x^2-x-4)=0,
∴矩陣有三個特徵值:2,(1±根號17)/2。把特徵值分別代入方程,設x=(a,b,c),可得到對於x=2,b=0,a+c=0,對應x=2的特徵向量為(-1,0,1)(未歸一化),其它x的一樣做。
求矩陣的全部特徵值和特徵向量:
1、計算的特徵多項式;
2、求出特徵方程的全部根,即為的全部特徵值;
3、對於的每一個特徵值,求出齊次線性方程組:的一個基礎解系,則的屬於特徵值的全部特徵向量是(其中是不全為零的任意實數)
[注]:若是的屬於的特徵向量,則也是對應於的特徵向量,因而特徵向量不能由特徵值惟一確定。反之,不同特徵值對應的特徵向量不會相等,亦即一個特徵向量只能屬於一個特徵值。
以上內容參考:網路-特徵值
❼ 如何用C語言求一般矩陣的特徵值和特徵向量
C語言並沒有封裝這類函數,只能自己實現。MATLAB倒是可以直接求。
自己實現的話可以用雅克比迭代法、高斯-賽戴爾迭代法等演算法
❽ 怎樣求特徵值和特徵向量
求特徵值的傳統方法是令特徵多項式| AE-A| = 0,求出A的特徵值,對於A的任一特徵值h,特徵方程( aE- A)X= 0的所有非零解X即為矩陣A的屬於特徵值N的特徵向量兩者的計算是分割的,一個是計算行列式,另一個是解齊次線性方程組,且計算量都較大。使用matlab可以方便的計算任何復雜的方陣的特徵值和特徵向量:
1、首先需要知道計算矩陣的特徵值和特徵向量要用eig函數,可以在命令行窗口中輸入help eig,查看一下eig函數的用法,如下圖所示:
注意事項:
特徵值和特徵向量的應用:
1、可以用在研究物理、化學領域的微分方程、連續的或離散的動力系統中。例如,在力學中,慣量的特徵向量定義了剛體的主軸。慣量是決定剛體圍繞質心轉動的關鍵數據;
2、數學生態學家用來預測原始森林遭到何種程度的砍伐,會造成貓頭鷹的種群滅亡;
3、著名的圖像處理中的PCA方法,選取特徵值最高的k個特徵向量來表示一個矩陣,從而達到降維分析+特徵顯示的方法,還有圖像壓縮的K-L變換。再比如很多人臉識別,數據流模式挖掘分析等方面。
❾ 如何用C語言求一般矩陣(非對稱矩陣)的特徵值和特徵向量
用C++或者VB編程很煩人的,matlab中命令:[a,b]=eig(A)就是求解矩陣A的特徵值和特徵值對應的向量,他們分別會構成一個由特徵值組成的對角矩陣b和一個由對應特徵值的特徵列向量組成的a矩陣。或者命令a=eig[A]就只有特徵值組成的對角矩陣a,別去想用C++和VB之類的,這些軟體用來求解矩陣和matlab相差太遠了。我之前也想過編程解決,人家一個命令就能解決的問題何不取巧呢?