Документ взят из кэша поисковой машины. Адрес оригинального документа : http://num-anal.srcc.msu.ru/meth_mat/prac_alg/arfil/symc1.txt Дата изменения: Tue Dec 17 12:59:03 2002 Дата индексирования: Mon Oct 1 23:20:03 2012 Кодировка: Windows-1251

#include
#include

double fdet ( float d, float c, int n );
int fau ( float d, float c, int n, double ua[] );
void fld ( float d, float c, int n, double ld[] );
void fv ( int n, double vld[] );

void main()
{
int i, i1, j, j1, n = 5 ;
float d = 2., c = 1., det;
double ld[ 10 ], ua[100], vld[100];

printf("\n НАЧАЛЬНЫЕ ДАННЫЕ \n");
printf(" d = %5.2f c = %5.2f n = %5d \n", d, c, n );

det = fdet ( d, c, n );
printf("\n Детерминант \n" );
printf(" det = %5.2f ", det );

fau ( d, c, n, ua );
printf("\n Обратная матрица \n" );
for ( i=0; i < n; i++)
{
for ( j=0; j < n; j++)
{
i1 = i + 1;
j1 = j + 1;
printf("ua(%1d,%1d)=%2.4f ", i1, j1, ua[j + i*n] );
}
printf("\n " );
}

fld (d, c, n, ld );
printf("\n Собственные значения \n" );
for ( i=0; i < n; i++)
{
i1 = i + 1;
printf(" ld(%1d) = %2.2f ", i1, ld[i] );
}

fv ( n, vld );
printf("\n Собственные векторы \n" );
for ( i=0; i < n; i++)
{
for ( j=0; j < n; j++)
{
i1 = i + 1;
j1 = j + 1;
printf(" vld(%1d,%1d) = %5.2f ", i1, j1, vld[j + i*n] );
}
printf("\n " );
}
} // end_of_main

double fdet ( float d, float c, int n )
// Вычисление определителя
// симметричной якобиевой матрицы
// с одинаковыми диагональными элементами.
//
// Параметры программы:
// d - диагональный элемент
// c - внедиагональный элемент
// n - порядок матрицы
// функция возвращает значение определителя
//
{
int i;
double det, eps = 1.e-10, fi, p, r, r1, r2, r3, s1, s2;
if ( fabs (d*d - 4*c*c) <= eps ) // совпадающие корни
// характеристического многочлена
{
for ( i=0, r = 1.; i < n; i++)
r *= d/2.;
det = r*( 1.+ n );
return det;
}
if ( d*d > 4*c*c) // различные вещественные корни
// характеристического многочлена
{
r1 = sqrt( d*d - 4.*c*c );
r2 = ( d + r1 );
r3 = ( d - r1 );
for ( i=0, s1 = 1., s2 = 1., p = 1.; i < n + 1; i++)
{
s1 *= r2;
s2 *= r3;
p *= 2.;
}
det = ( s1 - s2 )/( p*r1);
return det;
}
if ( d*d < 4*c*c) // мнимые корни
// характеристического многочлена
for ( i=0, r = 1.; i < n; i++)
r *= c;
r1 = d/( 2 * c );
fi = acos( r1 );
det = r * sin( ( n + 1.)*fi)/sin(fi);
return det;
} // end_of fdet

int fau ( float d, float c, int n, double ua[] )
// Вычисление обратной матрицы для
// симметричной якобиевой матрицы
// с одинаковыми диагональными элементами.
//
// Параметры программы:
// d - диагональный элемент
// c - внедиагональный элемент
// n - порядок матрицы
// ua - массив элементов обратной матрицы
// Функция возвращает 1,если матрица вырождена и
// 0 в противном случае.
{
int i, j, k, l ;
double c1, c2, q1, q2, eps = 1.e-10, fi, r, r1, r2, r3, s1, s2, xjk;
double sfi1, sfi2, sfi3;
if ( fabs( d*d - 4*c*c) <= eps ) // совпадающие корни
// характеристического многочлена
{
q1 = -d/(2.*c);
for ( j = 1 ; j < n + 1; j++ )
for ( k = 1 ; k < n + 1; k++ )
{
if ( k + 1 >= j )
{
for ( l = 0, r = 1.; l < k - j + 1; l++ )
r *= q1;
}
if ( k + 1 < j )
{
for ( l = 0, r = 1.; l < j - k - 1; l++ )
r /= q1;
}
r1 = c*( j*q1*q1 - j + n + 1 );
xjk = r/r1;
if ( k < j )
xjk *= ( j - n - 1 )*k;
else
xjk *= ( k - n - 1 )*j;
ua[k-1+n*(j-1)] = xjk;
}
return 0;
}
if ( d*d > 4*c*c ) // различные вещественные корни
// характеристического многочлена
{
r1 = sqrt( d*d - 4.*c*c );
q1 = ( -d + r1 )/(2.*c);
q2 = ( -d - r1 )/(2.*c);

for ( l = 0, r1 = 1., r2 = 1.; l < n + 1 ; l++ )
{
r1 *= q1;
r2 *= q2;
}

c1 = -1./( c*( q1 - q2 ) * ( r1 - r2 ));

for ( j = 1 ; j < n + 1; j++ )
{
for ( l = 0, r1 = 1., r2 = 1.; l < n + 1 - j ; l++ )
{
r1 *= q1;
r2 *= q2;
}
c2 = c1*( r1 - r2 );

for ( k = 1 ; k < n + 1; k++ )
{
for ( l = 0, r1 = 1., r2 = 1.; l < k ; l++ )
{
r1 *= q1;
r2 *= q2;
}
xjk = c2*( r1 - r2 );
ua[ k - 1 + n*(j - 1) ] = xjk;
if ( k > j )
{
for ( l = 0, r1 = 1., r2 = 1.; l < k - j ; l++ )
{
r1 *= q1;
r2 *= q2;
}
xjk += ( r1 - r2 )/( c * ( q1 - q2 ) );
ua[ k - 1 + n*(j - 1) ] = xjk;
}
}
}
return 0;
}

if ( d*d < 4*c*c ) // мнимые корни
// характеристического многочлена
{
r1 = -d/ ( 2 * c );
fi = acos( r1 );
for ( j = 1 ; j < n + 1; j++ )
{
for ( k = 1 ; k < n + 1; k++ )
{
r1 = sin( fi );
r2 = sin( (n + 1)*fi );
r1 *= c*r2;
if ( fabs(r1) <= eps ) return 1;
r2 = sin(( n + 1 - j )*fi);
r3 = sin( k*fi );
xjk = -r2*r3/r1;
if ( k > j )
{
r1 = sin( (k - j )*fi );
r2 = c*sin( fi );
if ( fabs(r2) <= eps ) return 1;
xjk += r1/r2;
}
ua[k - 1 + n*(j - 1)] = xjk;
}
}
}
return 0;
}

void fld ( float d, float c, int n, double ld[] )
// Вычисление собственных значений для
// симметричной якобиевой матрицы
// с одинаковыми диагональными элементами.
//
// Параметры программы:
// d - диагональный элемент
// c - внедиагональный элемент
// n - порядок матрицы
// ld - массив собственных значений
//
{
int i;
float k ;
double a, pi, pin1;
pi = atan(1.)*4.;
pin1=pi/(n+1);
if(c>0.)
{
for ( k=1; k <= n; k++)
{
ld[k-1] = d - 2*c*cos(pin1*k);
} return;
}
for ( k=n; k >= 1; k--)
{
ld[n-k] = d - 2*c*cos(pin1*k);
}return;
}

void fv ( int n, double vld[] )
// Вычисление собственных векторов
// симметричной якобиевой матрицы
// с одинаковыми диагональными элементами.
//
// Параметры программы:
// n - порядок матрицы
// vld - массив собственных векторов
{
int j, k;
double a, fik, fik_j, a1, pi;
{
pi = atan(1.)*4.;
for(k=1; k<=n; k++)
{
vld[(k-1)*n+0]=1.;
fik = pi*k/(n + 1.);
a1 = sin(fik);
for(j=2; j<=n; j++)
{
fik_j = fik*j;
a = sin(fik_j);
vld[(k-1)*n+(j-1)] = pow(-1,j+1)*a/a1;
}
}
}
}