Автор работы: Пользователь скрыл имя, 19 Января 2012 в 12:42, реферат
Цель курсовой работы: Ознакомиться и изучить проблему собственных значений исследование эффективных методов ее решения проблемы собственных значений.
Задачами курсовой работы являются: Определение эффективности работы алгоритмов решения проблемы собственных значений.
Введение
1 Классификация методов решения проблем собственных значений.
2 Эффективные методы решения проблем собственных значений.
2.1 Метод Якоби для действительных симметрических матриц.
2.2 Преобразование Хаусхолдера для приведения симметричной матрице к трехдиагональной форме.
2.3 QL-алгоритм с неявными сдвигами.
3 Программная реализация методов решения проблем собственных значений.
3.1 Метод Якоби для действительных симметрических матриц.
3.2 Преобразование Хаусхолдера для приведения симметричной матрице к трехдиагональной форме.
3.3 QL-алгоритм с неявными сдвигами.
4. Проведение численных экспериментов и анализ эффективности.
4.1 Метод Якоби для действительных симметрических матриц.
4.2 Преобразование Хаусхолдера для приведения симметричной матрице к трехдиагональной форме.
4.3 QL-алгоритм с неявными сдвигами.
Заключение
Список использованных источников
include "f2c.h"
static integer c__35 = 35;
int aej2d_c(integer *n, doublereal *a, doublereal *ev,
doublereal *rab1, integer *ierr)
{
static doublereal sys051 = 1.1107652e-16;
/* System generated locals */
integer a_dim1, a_offset, i__1, i__2;
doublereal d__1, d__2, d__3;
/* Builtin functions */
double sqrt(doublereal), d_sign(doublereal *, doublereal *);
/* Local variables */
static integer b, c__, d__, e, f, g, h__;
static doublereal i__, j, k, l, m, o, p;
extern /* Subroutine */ int utae10_c(integer *, integer *, integer *);
#define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
/* Parameter adjustments */
--rab1;
--ev;
a_dim1 = *n;
a_offset = a_dim1 + 1;
a -= a_offset;
/* Function Body */
*ierr = 0;
i__1 = *n;
for (b = 1; b <= i__1; ++b) {
if (b == 1) {
goto L3;
}
rab1[b] = a_ref(b, 1) * a_ref(b - 1, 3);
if ((d__1 = rab1[b]) < 0.) {
goto L5;
} else if (d__1 == 0) {
goto L1;
} else {
goto L2;
}
L1:
if (a_ref(b, 1) == 0. && a_ref(b - 1, 3) == 0.) {
goto L2;
}
*ierr = -(*n * 3 + b);
L2:
rab1[b] = sqrt(rab1[b]);
L3:
ev[b] = a_ref(b, 2);
/* L4: */
}
goto L6;
L5:
*ierr = *n + b;
L6:
if (*ierr != 0) {
goto L20;
}
*ierr = 0;
if (*n == 1) {
goto L20;
}
i__1 = *n;
for (c__ = 2; c__ <= i__1; ++c__) {
/* L7: */
rab1[c__ - 1] = rab1[c__];
}
rab1[*n] = 0.;
i__1 = *n;
for (e = 1; e <= i__1; ++e) {
d__ = 0;
L8:
i__2 = *n;
for (f = e; f <= i__2; ++f) {
if (f == *n) {
goto L10;
}
if ((d__1 = rab1[f], abs(d__1)) <= sys051 * ((d__2 = ev[f], abs(
d__2)) + (d__3 = ev[f + 1], abs(d__3)))) {
goto L10;
}
/* L9: */
}
L10:
m = ev[e];
if (f == e) {
goto L14;
}
if (d__ == 30) {
goto L19;
}
++d__;
l = (ev[e + 1] - m) / (rab1[e] * 2.);
if (abs(l) > 1.) {
/* Computing 2nd power */
d__1 = 1. / l;
o = abs(l) * sqrt(d__1 * d__1 + 1.);
}
if (abs(l) > 1.) {
/* Computing 2nd power */
d__1 = 1 / l;
o = abs(l) * sqrt(d__1 * d__1 + 1.);
}
if (abs(l) <= 1.) {
o = sqrt(l * l + 1.);
}
l = ev[f] - m + rab1[e] / (l + d_sign(&o, &l));
p = 1.;
j = 1.;
m = 0.;
h__ = f - e;
i__2 = h__;
for (g = 1; g <= i__2; ++g) {
c__ = f - g;
k = p * rab1[c__];
i__ = j * rab1[c__];
if (abs(k) < abs(l)) {
goto L11;
}
j = l / k;
o = sqrt(j * j + 1.);
rab1[c__ + 1] = k * o;
p = 1. / o;
j *= p;
goto L12;
L11:
p = k / l;
o = sqrt(p * p + 1.);
rab1[c__ + 1] = l * o;
j = 1. / o;
p *= j;
L12:
l = ev[c__ + 1] - m;
o = (ev[c__] - l) * p + j * 2. * i__;
m = p * o;
ev[c__ + 1] = l + m;
l = j * o - i__;
/* L13: */
}
ev[e] -= m;
rab1[e] = l;
rab1[f] = 0.;
goto L8;
L14:
if (e == 1) {
goto L16;
}
i__2 = e;
for (g = 2; g <= i__2; ++g) {
c__ = e + 2 - g;
if (m >= ev[c__ - 1]) {
goto L17;
}
ev[c__] = ev[c__ - 1];
/* L15: */
}
L16:
c__ = 1;
L17:
ev[c__] = m;
/* L18: */
}
goto L20;
L19:
*ierr = e;
L20:
if (*ierr != 0) {
utae10_c(ierr, n, &c__35);
}
return 0;
} /* aej2d_c */
#undef
a_ref
}
2.Преобразование Хаусхолдера для приведения симметричной матрице к трехдиагональной форме.
#include <math.h>
void tred2(float **a, int n, float *d, float *e) {
int l,k,j,i;
float scale,hh,h,g,f;
for(i=n;i>=2;i--) {
l=i-1; h=scale=0.;
for(k=1;k<=l;k++) scale += fabs(a[i][k]);
if(scale==0.) e[i]=a[i][l];
else {
for(k=1;k<=l;k++) {
a[i][k]/=scale; h += a[i][k]*a[i][k];
}
f=a[i][l];
g=(f>=0.?-sqrt(h):sqrt(h));
e[i]=scale*g; h -= f*g;
a[i][l]=f-g;
f=0.;
for(j=1;j<=l;j++) { a[j][i]=a[i][j]/h;
g=0.;
for(k=1;k<=j;k++) g += a[j][k]*a[i][k];
for(k=j+1;k<=l;k++) g += a[k][j]*a[i][k];
e[j]=g/h;
f += e[j]*a[i][j];
}
hh=f/(h+h);
for(j=1;j<=l;j++) {
/* Сформировать q и поместить на место
p (в e) */
f=a[i][j]; e[j]=g=e[j]-hh*f;
for(k=1;k<=j;k++) a[j][k] -= (f*e[k]+g*a[i][k]);
}
}
}
else e[i]=a[i][l];
d[i]=h;
}
d[1]=0.;
e[1]=0.;
for(i=1;i<=n;i++) {
l=i-1; if(d[i]!=0.) {
for(j=1;j<=l;j++) {
g=0.;
for(k=1;k<=l;k++) g += a[i][k]*a[k][j];
for(k=1;k<=l;k++) a[k][j] -= g*a[k][i];
}
} d[i]=a[i][i];
a[i][i]=0.;
for(j=1;j<=l;j++) a[j][i]=a[i][j]=0.;
}
}
#include "f2c.h"
static integer c__34 = 34;
int aee2r_c(integer *n, real *a, real *ev, real *rab1,
integer *ierr)
{
static real sys051 = 1.192093e-7f;
integer a_dim1, a_offset, i__1, i__2;
real r__1, r__2, r__3;
/* Builtin functions */
double sqrt(doublereal), r_sign(real *, real *);
/* Local variables */
static integer b, c__, d__, e, f, g;
static real h__;
static integer i__;
static real j, k, l, m, o, p;
extern /* Subroutine */ int utae10_c(integer *, integer *, integer *);
#define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
/* Parameter adjustments */
--rab1;
--ev;
a_dim1 = *n;
a_offset = a_dim1 + 1;
a -= a_offset;
/* Function Body */
i__1 = *n;
for (i__ = 1; i__ <= i__1; ++i__) {
ev[i__] = a_ref(i__, 2);
rab1[i__] = a_ref(i__, 1);
/* L1: */
}
*ierr = 0;
if (*n == 1) {
goto L15;
}
i__1 = *n;
for (b = 2; b <= i__1; ++b) {
/* L2: */
rab1[b - 1] = rab1[b];
}
rab1[*n] = 0.f;
i__1 = *n;
for (d__ = 1; d__ <= i__1; ++d__) {
c__ = 0;
L3:
i__2 = *n;