Wielomiany Czebyszowa ortogonalizacja
Mariusz:
using System;
using System.Globalization;
using System.Numerics;
using ExtendedNumerics;
namespace NamespaceName
{
public class ClassName
{
public static void Main(string[] args)
{
int n;
BigRational[,] T,U;
BigRational sum;
int error;
char esc;
do
{
Console.WriteLine("Podaj stopień macierzy");
int.TryParse(Console.ReadLine(),out n);
T = new BigRational[n+1,n+1];
U = new BigRational[n,n];
ModifiedGramSchmidt(n,T);
for(int i=0;i<=n;i++)
{
sum = 0;
for(int j=0; j<=i;j++)
sum += T[j,i];
for(int j=0; j<=i;j++)
T[j,i] /= (BigRational)sum;
}
for(int j=0;j<U.GetLength(1);j++)
for(int i=0;i<U.GetLength(0);i++)
{
U[i,j] = (i+1)*T[i+1,j+1];
U[i,j] /= (j+1);
}
Console.WriteLine("Macierz współczynników wielomianów Czebyszowa I rodzaju");
PrintMatrix(T);
Console.WriteLine("Macierz współczynników wielomianów Czebyszowa II rodzaju");
PrintMatrix(U);
Inverse(T,out error);
if(error == 0)
{
Console.WriteLine("Macierz odwrotna do macierzy współczynników wielomianów Czebyszowa I
rodzaju");
PrintMatrix(T);
}
esc = (char)Console.ReadKey().Key;
}
while(esc != (char)ConsoleKey.Escape);
}
public static void PrintMatrix(BigRational[,] A)
{
NumberFormatInfo nfi = new NumberFormatInfo();
nfi.NumberDecimalDigits = 12;
nfi.NumberGroupSeparator = "";
for(int i=0;i<A.GetLength(0);i++)
{
for(int j=0;j<A.GetLength(1);j++)
{
Console.Write("{0} , ",A[i,j].GetImproperFraction());
}
Console.WriteLine();
}
Console.WriteLine();
}
public static BigRational InnerProduct(int m, int n,BigRational[] P,BigRational[] Q)
{
BigRational sum = BigRational.Zero;
BigRational prod = BigRational.One;
BigRational[] C = new BigRational[n + m + 1];
for(int i = 0;i < C.Length;i++)
C[i] = BigRational.Zero;
for(int i = 0;i <= m;i++)
{
for(int j = 0; j <= n;j++)
{
C[i + j] += P[i]*Q[j];
}
}
for(int k = 0;k <= m + n;k += 2)
{
sum += C[k] * prod;
prod *= (BigRational)(k + 1)/(k + 2);
}
return sum;
}
public static BigRational[] GetColumn(BigRational[,] A,int j)
{
BigRational[] temp = new BigRational[A.GetLength(1)];
for(int k = 0;k < A.GetLength(1);k++)
{
temp[k] = A[k,j];
}
return temp;
}
public static void Inverse(BigRational[,] A,out int error)
{
if(A.GetLength(0) != A.GetLength(1))
{
error = 2;
return;
}
error = 0;
int n = A.GetLength(0);
int k;
BigRational maxA,d,e;
int[] M = new int[n];
for(int i = 0;i < n;i++)
{
maxA = 0;
k = i;
for(int j = i;j < n;j++)
{
d = A[j,i];
if(BigRational.Abs(maxA) < BigRational.Abs(d))
{
maxA = d;
k = j;
}
}
if(maxA == BigRational.Zero)
{
error = 1;
return;
}
M[i] = k;
A[k,i] = 1;
for(int j=0;j < n;j++)
{
d = (BigRational)A[k,j]/maxA;
A[k,j] = A[i,j];
A[i,j] = d;
}
for(int j = 0;j<n;j++)
if(j != i)
{
d = A[j,i];
A[j,i] = 0;
for(int l = 0;l < n;l++)
{
e = d * A[i,l];
A[j,l] −= e;
}
}
}
for(int i = n−1;i >= 0;i−−)
{
k = M[i];
if(k != i)
for(int j = 0;j < n;j++)
{
d = A[j,i];
A[j,i] = A[j,k];
A[j,k] = d;
}
}
}
public static void ModifiedGramSchmidt(int n, BigRational[,] Q)
{
/*
for j = 1:n
v
j = x
j;
end for
for j = 1:n
q
j = v
j / ||v
j||
2
for k = j + 1:n
v
k = v
k − (q
jTv
k)*q
j
end for
end for
*/
BigRational[] Qj = new BigRational[n+1];
BigRational[] Qk = new BigRational[n+1];
for(int i=0;i<=n;i++)
for(int j=0;j<=n;j++)
Q[i,j] = (i==j?1:0);
for(int j = 0;j <= n;j++)
{
Qj = GetColumn(Q,j);
BigRational invnorm = BigRational.One/(InnerProduct(j,j,Qj,Qj));
for(int k = j+1;k <= n;k++)
{
Qk = GetColumn(Q,k);
BigRational w = invnorm*InnerProduct(j,k,Qj,Qk);
for(int i = k;i >= 0;i−−)
{
Q[i,k] −= w*Q[i,j] ;
}
}
}
}
}
}
Tutaj zmodyfikowałem nieco iloczyn skalarny aby uniknąć
konieczności korzystania z liczb niewymiernych
Użyłem tutaj zmodyfikowanej metody ortogonalizacji Grama−Schmidta
Ciekawe jak ją zastąpić np metodą odbić Householdera
Dość kosztowny czasowo jest też iloczyn skalarny
czy nie dałoby się go jakoś przyśpieszyć ?
Skompilowałem ze źródeł znalezionych na githubie klasę BigRational do dll