如何用C語言編寫求對稱矩陣的特徵值和特徵向量的程式

2021-03-11 09:38:59 字數 5691 閱讀 4163

1樓:空麼

//數值計算程式-特徵值和特徵向量

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

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

if (h+1.0==1.0)

b[i]=0.0;

} else

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

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

f=0.0;

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

if (j+1<=i-1) }

c[j-1]=g/h;

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

} h2=f/(h+h);

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

b[i]=h;

} }b[0]=0.0;

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

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

} u=i*n+i;

b[i]=q[u];

q[u]=1.0;

if (i-1>=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)

m=j;

while ((m<=n-1)&&(fabs(c[m])>d))

if (m!=j)

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)

else

h=g-b[j];

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

f=f+h;

p=b[m];

e=1.0;

s=0.0;

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

else

p=e*b[i]-s*g;

b[i+1]=h+s*(e*g+s*b[i]);

for (k=0; k<=n-1; k++) }

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++) }

if (k!=i) }

} return(1);

} //約化實矩陣為赫申伯格(hessen berg)矩陣

//利用初等相似變換將n階實矩陣約化為上h矩陣

//a-長度為n*n的陣列,存放n階實矩陣,返回時存放上h矩陣

//n-矩陣的階數

void echbg(double a,int n) }

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

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

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

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

} }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)

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)

else if (l==m-2)

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

m=m-2;

it=0;

} else

it=it+1;

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

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

for (k=l; k<=m-2; k++) }

else

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

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

if (k!=l)

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++)

a[jj]=q;

a[ii]=p;

} j=k+3;

if (j>=m-1)

for (i=l; i<=j; i++)

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) }

} while (1==1) }

} if (fmjt)

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)

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

sn=omega/sqrt(2.0*sn);

**=sqrt(1.0-sn*sn);

fm=a[w];

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

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

a[u]=0.0;

a[t]=0.0;

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

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

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

return(1);}

c語言求矩陣的逆,C語言 求矩陣的逆

源程式如下 include include include include include include define max 100void inputstyle int 輸入函式 void input int int 輸入函式 long danx int int int sgnx int vo...

c語言中如何用多個檔案編寫程式C語言中,如何用多個檔案編寫程式

將一個函式寫在一個檔案裡,然後再在另一個檔案裡用 include 包含這個檔案。舉個例子 在檔案c1.c裡編一個函式 void printword 再建立一個檔案c2.c,檔案開頭寫上 include c1.c 就可以呼叫c1.c裡的函式printword 了 舉個簡單的例子!你會容易理解的。你寫一...

用c語言怎麼編寫輸入矩陣求其逆矩陣的程式

迴圈輸入矩陣元素,你想想求行列式的演算法,改一改就是求逆矩陣 通過 a e e a 1 這個初等變換來求逆矩陣。下面是實現gauss jordan法實矩陣求逆。include include include intbrinv double a,intn if d 1.0 1.0 if is k k ...