Kod: Zaznacz cały
// Example program
#include <iostream>
#include <iomanip>
#include <cmath>
#include <ctime>
const int maxiter = 1000;
const double eps = 1e-12;
double abs(double x)
{
return (x < 0? -x:x);
}
void elimHess(int n,double **A)
{
int i,j,l,k;
double m,maxA,temp;
for(i = 0;i < n - 2;i++)
{
k = i + 1;
maxA = abs(A[k][i]);
for(j = i + 1;j < n;j++)
if(abs(A[j][i]) > maxA)
{
k = j;
maxA = abs(A[k][i]);
}
if(maxA > 0)
{
for(j = 0;j < n;j++)
{
temp = A[k][j];
A[k][j] = A[i + 1][j];
A[i + 1][j] = temp;
}
for(j = 0;j < n;j++)
{
temp = A[j][k];
A[j][k] = A[j][i + 1];
A[j][i + 1] = temp;
}
for(j = i + 2;j < n;j++)
{
m = (double)A[j][i]/A[i + 1][i];
for(l = 0; l < n;l++)
A[j][l] -= m * A[i + 1][l];
for(l = 0;l < n;l++)
A[l][i+1] += m * A[l][j];
}
}
}
}
void QRdecomp(int n,int m,double **A,double **Q)
{
double r,c,s,temp;
int i,j,k,min;
for(i = 0;i < m;i++)
for(j = 0;j < m;j++)
Q[i][j] = (i == j ? 1 : 0);
min = (m < n ? m : n);
for(i = 0;i < min;i++)
for(j = i + 1;j < m; j++)
if(A[j][i] != 0)
{
r = hypot(A[i][i],A[j][i]);
c = (double)A[i][i]/r;
s = (double)A[j][i]/r;
for(k = 0;k < n; k++)
{
temp = A[i][k];
A[i][k] = c * A[i][k] + s * A[j][k];
A[j][k] = -s * temp + c * A[j][k];
}
for(k = 0 ;k < m;k++)
{
temp = Q[k][i];
Q[k][i] = c * Q[k][i] + s * Q[k][j];
Q[k][j] = -s * temp + c * Q[k][j];
}
}
}
void multiplyMatrix(int n,int m ,int p,double **A,double **B,double **C)
{
int i,j,k;
for(i=0;i<n;i++)
for(k = 0;k < p;k++)
{
C[i][k] = 0;
for(j = 0;j < n;j++)
C[i][k] += A[i][j] * B[j][k];
}
}
void copyMatrix(int m,int n, double **A,double **B)
{
for(int i=0; i< m;i++)
for(int j=0;j < n;j++)
B[i][j] = A[i][j];
}
void printMatrix(int m,int n,double**A)
{
for(int i=0;i<n;i++)
{
for(int j=0;j<n;j++)
if(abs(A[i][j]) > eps)
std::cout<<std::setprecision(12)<<A[i][j]<<" ";
else
std::cout<<"0 ";
std::cout<<"\n";
}
std::cout<<"\n";
}
int main()
{
int m,n,k;
double **A;
double **Q;
double **R;
double p,q,d;
std::cout<<"Podaj stopien macierzy \n";
std::cin>>n;
std::cout<<"Podaj gorny zakres przedzialu dla elementow macierzy \n";
std::cin>>m;
A = new double*[n];
for(int i=0;i<n;i++)
A[i] = new double[n];
Q = new double*[n];
for(int i=0;i<n;i++)
Q[i] = new double[n];
R = new double*[n];
for(int i=0;i<n;i++)
R[i] = new double[n];
srand(time(NULL));
for(int i=0;i<n;i++)
for(int j=0;j<n;j++)
A[i][j] = rand()%m;
std::cout<<"Wylosowana macierz \n";
printMatrix(n,n,A);
elimHess(n,A);
std::cout<<"Macierz po sprowadzeniu do postaci Hessenberga \n";
printMatrix(n,n,A);
for(int i=0;i<maxiter;i++)
{
QRdecomp(n,n,A,Q);
copyMatrix(n,n,A,R);
multiplyMatrix(n,n,n,R,Q,A);
}
std::cout<<"Macierz po iteracyjnym sprowadzeniu do postaci Schura \n";
printMatrix(n,n,A);
std::cout<<"Przyblizone wartosci wlasne macierzy \n";
k = 0;
while(k < n)
{
if(k + 1 < n && abs(A[k+1][k]) > eps)
{
p = 0.5 * (A[k][k] + A[k + 1][k + 1]);
q = A[k][k] * A[k + 1][k + 1] - A[k][k + 1] * A[k+1][k];
d = sqrt(abs(q - p * p));
std::cout<<"x["<<k<<"]="<<std::setprecision(12)<<p<<"-"<<std::setprecision(12)<<d<<"*i \n";
std::cout<<"x["<<k+1<<"]="<<std::setprecision(12)<<p<<"+"<<std::setprecision(12)<<d<<"*i \n";
k += 2;
}
else
{
std::cout<<"x["<<k<<"]="<<std::setprecision(12)<<A[k][k]<<"\n";
k++;
}
}
for(int i=0;i<n;i++)
delete[] A[i];
delete[] A;
for(int i=0;i<n;i++)
delete[] Q[i];
delete[] Q;
for(int i=0;i<n;i++)
delete[] R[i];
delete[] R;
return 0;
}
opisana w książce Fortuna, Macukow, Wąsowski Metody numeryczne
ale eliminacja Gaussa może być niestabilna numerycznie więc lepiej by było zastosować odbicia Householdera
Jak mógłby wyglądać pseudokod redukcji do postaci Hessenberga z wykorzystaniem odbić Householdera ?
Pozostałe pytania to
Jakie wybrać przesunięcie
(Macierz ma elementy rzeczywiste ale wartości własne mogą być zespolone
obliczane z bloków czy jak kto woli klatek 2x2 występujących na przekątnej)
Jak zrealizować deflację i czy można by ją zrealizować w miejscu (bez allokacji dodatkowej pamięci)
Zbieżność dla wielokrotnych wartości własnych jest bardzo wolna
Jakieś pomysły na przyśpieszenie zbieżności ?
Jaki dać warunek stopu , inny niż maksymalna liczba iteracji