Jako że nie ma działu dla metod numerycznych umieszczam zadanie tutaj
Próbowałem napisać samodzielnie program realizujący tę metodę
Co mam do tej pory
Sprowadzenie macierzy do postaci Hessenberga ale tylko za pomocą eliminacji Gaussa
bo była ona dobrze opisana w jednej z książek do metod numerycznych
Wolałbym jednak użyć do tej redukcji odbić Householdera bo eliminacja Gaussa bywa niestabilna numerycznie
Czym jest postać Hessenberga?
Macierz A jest w postaci Hessenberga jeżeli można ją zapisać w następujący sposób
A = U+T , gdzie U to macierz trójkątna górna a T macierz trójdiagonalna
Rozkład QR zrealizowałem za pomocą mnożenia przez macierze obrotów
Rozpisałem sobie mnożenie przez macierze obrotów i okazało się że lewostronne mnożenie modyfikuje tylko dwa wiersze
a prawostronne mnożenie modyfikuje tylko dwie kolumny
(Mnożymy lewostronnie macierz A przez macierze obrotów aby dostać macierz R
Mnożymy prawostronnie macierz I przez te same macierze obrotów aby dostać macierz Q)
Mnożenie macierzy niby łatwe - zastosowanie wzorku - trzy pętle ale aby dobrze działało musiałem macierz skopiować
Na razie algorytm wykonuje tylko z góry określoną liczbę iteracji
Podstawowa metoda QR wygląda następująco
A_{k} = Q_{k}R_{k}
A_{k+1} = R_{k}Q_{k}
Nietrudno pokazać że metoda ta tworzy nam ciąg macierzy podobnych do macierzy A
Jeżeli metoda ta jest zbieżna do postaci Schura to stosunkowo łatwo odczytać wartości własne
Macierz w postaci Schura jest macierzą trójkątną (dla macierzy o elementach zespolonych)
oraz taką macierzą w której na przekątnej występują co najwyżej bloki 2x2 , poniżej tej blokowej przekątnej występują zera
a powyżej dowolna liczba rzeczywista
(te bloki 2x2 występują tylko jeżeli wartości własne z nich policzone są zespolone a cała macierz ma tylko elementy rzeczywiste)
Problem w tym że ta metoda w podstawowej wersji nie zawsze jest zbieżna do postaci Schura
a nawet jeśli jest zbieżna to dość często jest zbieżna bardzo wolno
No i teraz jakie wybrać przesunięcie i jak je zrealizować?
Jak zrealizować deflację
Co ze zbieżnością dla wielokrotnych wartości własnych
Jaki warunek stopu , inny niż maksymalna liczba iteracji ?
A i jeszcze jedno chciałbym uniknąć arytmetyki zespolonej
using System;
using System.IO;
using System.Globalization;
namespace NamespaceName
{
public class ClassName
{
const int maxiter = 1000;
const double eps = 1e-12;
public static void PrintMatrix(int m,int n,double[,] A)
{
NumberFormatInfo nfi = new NumberFormatInfo();
nfi.NumberDecimalDigits = 12;
for(int i = 0; i < m; i++)
{
for(int j = 0; j < n; j++)
Console.Write("{0} ,",A[i,j].ToString("N",nfi));
Console.WriteLine();
}
Console.WriteLine();
}
public static void PrintMatrix(StreamWriter sw,int m,int n,double[,] A)
{
NumberFormatInfo nfi = new NumberFormatInfo();
nfi.NumberDecimalDigits = 12;
for(int i = 0; i < m; i++)
{
for(int j = 0; j < n; j++)
sw.Write("{0} ",A[i,j].ToString("N",nfi));
sw.WriteLine();
}
sw.WriteLine();
}
public static double Hypot(double a,double b)
{
double absa,absb;
absa = Math.Abs(a);
absb = Math.Abs(b);
if(absa > absb) return absa * Math.Sqrt(1 + ((double) absb / absa) * ((double) absb / absa));
else return (absb == 0 ? 0 : absb * Math.Sqrt(1 + ((double) absa / absb) * ((double) absa / absb)));
}
public static void ToHessenberg(int n,double[,] A)
{
for(int i = 0; i < n - 2;i++)
{
int k = i + 1;
double maxA = Math.Abs(A[k,i]);
for(int j = i + 1; j < n; j++)
if(Math.Abs(A[j,i]) > maxA)
{
k = j;
maxA = Math.Abs(A[k,i]);
}
if(maxA > 0)
{
for(int j = 0; j < n; j++)
{
double temp = A[k,j];
A[k,j] = A[i + 1,j];
A[i + 1,j] = temp;
}
for(int j = 0; j < n;j++)
{
double temp = A[j,k];
A[j,k] = A[j,i + 1];
A[j,i + 1] = temp;
}
for(int j = i + 2; j < n; j++)
{
double m = (double) A[j,i] / A[i + 1,i];
for(int l = 0; l < n; l++)
A[j,l] -= m * A[i + 1,l];
for(int l = 0;l < n;l++)
A[l,i + 1] += m * A[l,j];
}
}
}
}
public static void QR_Givens(int m,int n,double[,] A,double[,] Q)
{
for(int i = 0; i < m; i++)
for(int j = 0;j < m;j++)
Q[i,j] = (i == j ? 1 : 0);
int min = (m < n ? m : n);
for(int i = 0; i < min; i++)
for(int j = i + 1; j < m; j++)
if(A[j,i] != 0)
{
double r = Hypot(A[i,i],A[j,i]);
double c = (double) A[i,i] / r;
double s = (double) A[j,i] / r;
for(int k = 0; k < n; k++)
{
double 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(int k = 0;k < m;k++)
{
double temp = Q[k,i];
Q[k,i] = c * Q[k,i] + s * Q[k,j];
Q[k,j] = -s * temp + c * Q[k,j];
}
}
}
public static 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];
}
public static void MultiplyMatrix(int m,int n,int p,double[,] A,double[,] B,double[,] C)
{
for(int i = 0; i < m; i++)
for(int k = 0; k < p; k++)
{
C[i,k] = 0;
for(int j =0; j < n; j++)
C[i,k] += A[i,j] * B[j,k];
}
}
public static void Main(string[] args)
{
Console.Clear();
char esc;
double[,] A,Q,R;
int n;
string path;
Console.WriteLine("Podaj sciezke do pliku do ktorego chcesz zapisac wynik");
path = Console.ReadLine();
using(StreamWriter sw = new StreamWriter(path,true))
{
do
{
Console.WriteLine("Podaj stopień macierzy ktorej wartosci wlasne chcesz policzyc");
int.TryParse(Console.ReadLine(),out n);
A = new double[n,n];
R = new double[n,n];
Q = new double[n,n];
for(int i = 0;i < n;i++)
for(int j = 0;j < n;j++)
{
Console.Write("A[{0},{1}]=",i,j);
double.TryParse(Console.ReadLine(),out A[i,j]);
}
Console.WriteLine("Wczytana macierz A");
sw.WriteLine("Wczytana macierz A");
PrintMatrix(n,n,A);
PrintMatrix(sw,n,n,A);
ToHessenberg(n,A);
Console.WriteLine("Macierz A po sprowadzeniu do postaci Hessenberga");
PrintMatrix(n,n,A);
sw.WriteLine("Macierz A po sprowadzeniu do postaci Hessenberga");
PrintMatrix(sw,n,n,A);
for(int i = 0;i < maxiter; i++)
{
QR_Givens(n,n,A,Q);
CopyMatrix(n,n,A,R);
MultiplyMatrix(n,n,n,R,Q,A);
}
Console.WriteLine("Macierz A po iteracyjnym sprowadzeniu do postaci Schura");
sw.WriteLine("Macierz A po iteracyjnym sprowadzeniu do postaci Schura");
PrintMatrix(n,n,A);
PrintMatrix(sw,n,n,A);
Console.WriteLine("Przyblizone wartosci wlasne macierzy A");
sw.WriteLine("Przyblizone wartosci wlasne macierzy A");
int k = 0;
NumberFormatInfo nfi = new NumberFormatInfo();
nfi.NumberDecimalDigits = 12;
while(k < n)
{
if(k + 1 < n && Math.Abs(A[k+1,k]) > eps)
{
double p = 0.5 * (A[k,k] + A[k + 1,k + 1]);
double q = A[k,k] * A[k + 1,k + 1] - A[k,k + 1] * A[k + 1,k];
double d = Math.Sqrt(Math.Abs(q - p * p));
Console.WriteLine("x[{0}]={1} - {2}i",k,p.ToString("N",nfi),d.ToString("N",nfi));
Console.WriteLine("x[{0}]={1} + {2}i",k+1,p.ToString("N",nfi),d.ToString("N",nfi));
sw.WriteLine("x[{0}]={1} - {2}i",k,p.ToString("N",nfi),d.ToString("N",nfi));
sw.WriteLine("x[{0}]={1} + {2}i",k+1,p.ToString("N",nfi),d.ToString("N",nfi));
k += 2;
}
else
{
Console.WriteLine("x[{0}]={1}",k,A[k,k].ToString("N",nfi));
sw.WriteLine("x[{0}]={1}",k,A[k,k].ToString("N",nfi));
k++;
}
}
esc = (char) Console.ReadKey().Key;
}
while(esc != (char) ConsoleKey.Escape);
}
}
}
}