#include #include #define Nmax 11 void GJACOBI(int nf,double a[Nmax][Nmax],double b[Nmax][Nmax],double ev[Nmax],double vec[Nmax][Nmax]); /*--------------------------------------------------------------------*/ int main() { /* nf:固有方程式の次数(配列宣言はnf+1) */ /* ev[k]:固有値(大きい順に出力) */ /* vec[i=0〜n-1][k]:固有ベクトル */ int i,j,k; int nf=3; double a[Nmax][Nmax]; double b[Nmax][Nmax]; double ev[Nmax]; double vec[Nmax][Nmax]; /*検算用*/ double ar[Nmax][Nmax]; double br[Nmax][Nmax]; double c[Nmax][Nmax]; double d[Nmax][Nmax]; //Case-1 a[1][1]=4;a[1][2]=2;a[1][3]=3; a[2][1]=2;a[2][2]=6;a[2][3]=1; a[3][1]=3;a[3][2]=1;a[3][3]=5; b[1][1]=1;b[1][2]=0;b[1][3]=0; b[2][1]=0;b[2][2]=1;b[2][3]=0; b[3][1]=0;b[3][2]=0;b[3][3]=1; for(i=1;i<=nf;i++){ for(j=1;j<=nf;j++){ ar[i][j]=a[i][j]; br[i][j]=b[i][j]; } } GJACOBI(nf,a,b,ev,vec); printf("***Eigenvalue***\n"); printf("%16.6e %16.6e %16.6e\n",ev[1],ev[2],ev[3]); printf("***Eigenvector***\n"); for(i=1;i<=nf;i++){ printf("%16.6e %16.6e %16.6e\n",vec[i][1],vec[i][2],vec[i][3]); } for(i=1;i<=nf;i++){ for(j=1;j<=nf;j++){ c[i][j]=0.0; d[i][j]=0.0; for(k=1;k<=nf;k++){ c[i][j]=c[i][j]+ar[i][k]*vec[k][j]; d[i][j]=d[i][j]+br[i][k]*vec[k][j]*ev[j]; } } } printf("***検算結果 [A]{x}***\n"); for(i=1;i<=nf;i++){ printf("%16.6e %16.6e %16.6e\n",c[i][1],c[i][2],c[i][3]); } printf("***検算結果 λ[b]{x}***\n"); for(i=1;i<=nf;i++){ printf("%16.6e %16.6e %16.6e\n",d[i][1],d[i][2],d[i][3]); } /* Case-1出力結果 */ /* ***Eigenvalue*** */ /* 9.000000e+000 4.732051e+000 1.267949e+000 */ /* ***Eigenvector*** */ /* 9.227825e-001 -2.633891e-001 1.060632e+000 */ /* 9.227825e-001 9.829816e-001 -2.841954e-001 */ /* 9.227825e-001 -7.195925e-001 -7.764364e-001 */ printf("リターンキーを押して下さい"); getchar(); return 0; } /*--------------------------------------------------------------------*/ void GJACOBI(int nf,double a[Nmax][Nmax],double b[Nmax][Nmax],double ev[Nmax],double vec[Nmax][Nmax]) { /**************************************/ /* 一般化JACOBI法による固有値計算 */ /* [パソコン構造力学,宮澤健二著] */ /* 固有値ev[k],固有ベクトルvec[i][k] */ /**************************************/ int i,j,k,kmax,ip,iq,im; double fmk,fmg,epsk,epsg; double fact=1.0E-20; double amax,dam,de1,de2,de3,dd,ss,x,gg,ww; fmk=0.0; fmg=0.0; for(i=1;i<=nf;i++){ for(j=1;j<=nf;j++){ if(fabs(a[i][j])>fmk)fmk=fabs(a[i][j]); if(fabs(b[i][j])>fmg)fmg=fabs(b[i][j]); } } epsk=fmk*fact; epsg=fmg*fact; for(i=1;i<=nf;i++){ for(j=1;j<=nf;j++){ vec[i][j]=0.0; } vec[i][i]=1.0; } kmax=5*nf*nf; for(k=1;k<=kmax;k++){ ip=1; iq=2; amax=a[1][2]*a[1][2]/(a[1][1]*a[1][1]+a[2][2]*a[2][2]); for(i=1;i<=nf-1;i++){ for(j=i+1;j<=nf;j++){ dam=a[i][j]*a[i][j]/(a[i][i]*a[i][i]+a[j][j]*a[j][j]); if(amax=amax){ im=i;amax=ev[i]; } } ww=ev[k];ev[k]=ev[im];ev[im]=ww; for(i=1;i<=nf;i++){ ww=vec[i][k]; vec[i][k]=vec[i][im]; vec[i][im]=ww; } } } /*--------------------------------------------------------------------*/