Metoda QR przybliżonego obliczania wartości własnych macierzy

Otrzymałeś(aś) rozwiązanie do zamieszczonego zadania? - podziękuj autorowi rozwiązania! Kliknij
Robakks
Czasem tu bywam
Czasem tu bywam
Posty: 149
Rejestracja: 30 wrz 2012, 20:36
Podziękowania: 2 razy
Otrzymane podziękowania: 13 razy
Płeć:

Metoda QR przybliżonego obliczania wartości własnych macierzy

Post autor: Robakks »

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

Redukcję do postaci Hessenberga zrealizowałem za pomocą eliminacji Gaussa bo była ona dość dobrze
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
ODPOWIEDZ