Ortogonalizacja
Mariusz:
#include <iostream>
#include <vector>
#include <gmpxx.h>
#include "Polynomial.h"
mpq
class InnerProductT(const Polynomial<mpq
class>& P, const Polynomial<mpq
class>& Q)
{
Polynomial<mpq
class> C = P * Q;
mpq
class s = 0;
mpq
class p = 1;
for(int k = 0; k <= C.degree(); k += 2)
{
s += p * C[k];
p *= (mpq
class)(k+1)/(k+2);
}
return s;
}
mpq
class InnerProductU(const Polynomial<mpq
class>& P, const Polynomial<mpq
class>& Q)
{
Polynomial<mpq
class> C = P * Q * Polynomial<mpq
class>({1, 0, −1});
mpq
class s = 0;
mpq
class p = 1;
for(int k = 0; k <= C.degree(); k += 2)
{
s += p * C[k];
p *= (mpq
class)(k+1)/(k+2);
}
return s;
}
mpq
class InnerProductP(const Polynomial<mpq
class>& P, const Polynomial<mpq
class>& Q)
{
Polynomial<mpq
class> C = P * Q;
mpq
class s = 0;
for(int k = 0;k <= C.degree();k += 2)
{
s += C[k]*(mpq
class)2/(k+1);
}
return s;
}
mpq
class InnerProductH(const Polynomial<mpq
class>& P, const Polynomial<mpq
class>& Q)
{
Polynomial<mpq
class> C = P * Q;
mpq
class s = 0;
mpq
class p = 1;
for(int k = 0; k <= C.degree(); k += 2)
{
s += p * C[k];
p *= (mpq
class)(k+1)/2;
}
return s;
}
mpq
class InnerProductL(const Polynomial<mpq
class>& P, const Polynomial<mpq
class>& Q)
{
Polynomial<mpq
class> C = P * Q;
mpq
class s = 0;
mpq
class p = 1;
for(int k = 0; k <= C.degree(); k++)
{
s += p * C[k];
p *= (mpq
class)(k+1);
}
return s;
}
void ModifiedGramSchmidt(int n, std::vector<Polynomial<mpq
class>>& Q,
mpq
class (*InnerProduct)(const Polynomial<mpq
class>&, const
Polynomial<mpq
class>&))
{
Polynomial<mpq
class> Qj, Qk;
mpq
class invnorm, w;
Q.resize(n+1);
for(int i = 0; i <= n; i++)
{
Q[i] = Polynomial<mpq
class>::monomial(i);
}
for(int j = 0; j <= n; j++)
{
Qj = Q[j];
invnorm = (mpq
class)1 / InnerProduct(Qj, Qj);
for(int k = j + 1; k <= n; k++)
{
Qk = Q[k];
w = invnorm * InnerProduct(Qj, Qk);
Q[k] −= Polynomial<mpq
class>(w * Q[j]);
}
}
}
void GramSchmidt(int n, std::vector<Polynomial<mpq
class>>& Q,
mpq
class (*InnerProduct)(const Polynomial<mpq
class>&, const
Polynomial<mpq
class>&))
{
mpq
class mu;
Polynomial<mpq
class> Qi, Qj;
Q.resize(n+1);
for(int i = 0; i <= n; i++)
{
Q[i] = Polynomial<mpq
class>::monomial(i);
}
for(int i = 0; i <= n; i++)
{
Qi = Q[i];
for(int j = 0; j <= i − 1; j++)
{
Qj = Q[j];
mu = InnerProduct(Qi, Qj) / InnerProduct(Qj, Qj);
Qi −= Polynomial<mpq
class>(mu * Qj);
}
Q[i] = Qi;
}
}
int main()
{
int n;
std::cout << "Podaj n=";
std::cin >> n;
// Pierwszy rodzaj wielomianów Czebyszewa (T)
std::vector<Polynomial<mpq
class>> QT(n+1);
GramSchmidt(n, QT, InnerProductT);
// Normalizacja − suma współczynników = 1
for(int i = 0; i <= n; i++)
{
mpq
class s = 0;
for(int j = 0; j <= i; j++)
{
s += QT[i][j];
}
QT[i] *= (mpq
class)1/s;
}
std::cout << "Wielomiany Czebyszewa pierwszego rodzaju (T):\n";
for(const auto& poly : QT)
{
std::cout << poly << std::endl;
}
// Drugi rodzaj wielomianów Czebyszewa (U)
std::vector<Polynomial<mpq
class>> QU(n+1);
GramSchmidt(n, QU, InnerProductU);
// Mnożenie przez potęgi 2
mpq
class p2 = 1;
std::cout << "\nWielomiany Czebyszewa drugiego rodzaju (U):\n";
for(int i = 0; i <= n; i++)
{
QU[i] *= p2;
std::cout << QU[i] << std::endl;
p2 *= 2;
}
std::vector<Polynomial<mpq
class>> QP(n+1);
GramSchmidt(n, QP, InnerProductP);
// Normalizacja − suma współczynników = 1
for(int i = 0; i <= n; i++)
{
mpq
class s = 0;
for(int j = 0; j <= i; j++)
{
s += QP[i][j];
}
QP[i] *= (mpq
class)1/s;
}
std::cout<<std::endl;
std::cout << "Wielomiany Legendre'a (P):\n";
for(const auto& poly : QP)
{
std::cout << poly << std::endl;
}
std::vector<Polynomial<mpq
class>> QH(n+1);
GramSchmidt(n, QH, InnerProductH);
// Mnożenie przez potęgi 2
p2 = 1;
std::cout<<std::endl;
std::cout << "\nWielomiany Hermite'a (H):\n";
for(int i = 0; i <= n; i++)
{
QH[i] *= p2;
std::cout << QH[i] << std::endl;
p2 *= 2;
}
std::cout<<std::endl;
std::vector<Polynomial<mpq
class>> QL(n+1);
GramSchmidt(n, QL, InnerProductL);
mpq
class tmp;
std::cout<<std::endl;
std::cout << "\nWielomiany Laguerre'a (L):\n";
for(int i = 0; i <= n; i++)
{
tmp = 1/QL[i][0];
QL[i] *= tmp;
std::cout << QL[i] << std::endl;
}
system("PAUSE");
return 0;
}
Aby programik zadziałał musicie mieć gmp
(Dla Windows wystarczą pliki dll oraz pliki nagłówkowe)
oraz musicie napisać własną szablonową klasę wielomianów
Użyłem tutaj (zmodyfikowanej) ortogonalizacji Grama−Schmidta
Jeśli chodzi o iloczyny skalarne to:
Dla wielomianów T Czebyszowa iloczyn skalarny zdefiniowany jest następująco
| | p(x)q(x) | |
<p(x),q(x)> = ∫−11 |
| dx |
| | √1−x2 | |
Podczas pisania funkcji obliczającej iloczyn skalarny skorzystałem ze
| | xn | |
wzoru rekurencyjnego na całkę ∫−11 |
| dx |
| | √1−x2 | |
Aby otrzymać wielomiany T Czebyszowa trzeba jeszcze skorzystać z warunku T
n(1) = 1
Dla wielomianów U Czebyszowa iloczyn skalarny zdefiniowany jest następująco
<p(x),q(x)> = ∫
−11p(x)q(x)
√1−x2dx
Podczas pisania funkcji obliczającej iloczyn skalarny skorzystałem z tego że
| | 1−x2 | |
√1−x2 = |
| (stąd ten dodatkowy czynnik (1−x2)) |
| | √1−x2 | |
oraz z wcześniej wyprowadzonego wzoru redukcyjnego
Aby otrzymać wielomiany U Czebyszowa trzeba jeszcze skorzystać z warunku że
współczynnik wiodący wynosi 2
n
Dla wielomianów Legendre'a iloczyn skalarny zdefiniowany jest następująco
iloczyn skalarny zdefiniowany jest następująco
<p(x),q(x)> = ∫
−11p(x)q(x)dx
i łatwo napisać bez wzoru redukcyjnego
Aby otrzymać wielomiany Legendre'a trzeba jeszcze skorzystać z warunku T
n(1) = 1
Dla wielomianów Hermite'a iloczyn skalarny zdefiniowany jest następująco
<p(x),q(x)> = ∫
−∞∞p(x)q(x)e
−x2dx
Podczas pisania funkcji obliczającej iloczyn skalarny można skorzystać
ze wzoru rekurencyjnego na całkę ∫
−∞∞x
ne
−x2dx
Aby otrzymać wielomiany Hermite'a trzeba jeszcze skorzystać z warunku że
współczynnik wiodący wynosi 2
n
Do tej pory w funkcjach obliczających iloczyn skalarny zaczynałem pętle od zera
i w każdej iteracji zwiększałem licznik o dwa bo
funkcja podcałkowa dla nieparzystego n była nieparzysta a
przedział całkowania był symetryczny względem zera
Dla wielomianów Laguerre'a iloczyn skalarny zdefiniowany jest następująco
<p(x),q(x)> = ∫
0∞p(x)q(x)e
−xdx
Podczas pisania funkcji obliczającej iloczyn skalarny można skorzystać
ze wzoru rekurencyjnego na całkę ∫
0∞x
ne
−xdx
Aby otrzymać wielomiany Laguerre'a trzeba uwzględnić jeszcze to że L
n(0) = 1
Tutaj przy obliczaniu iloczynu skalarnego zacząłem pętle od zera
i w każdej iteracji zwiększałem licznik o jeden
Podczas obliczania iloczynu skalarnego najbardziej kosztowny jest iloczyn wielomianów
(Mnożenie znane ze szkoły wymaga O(n
2) operacji choć można ten koszt zmniejszyć)
W ortogonalizacji Grama−Schmidta mamy dwie zagnieżdżone pętle napisane w ten sposób
że operacje o koszcie O(1) dałyby w sumie złożoność O(n
2)
ale w tych zagnieżdżonych pętlach mamy mnożenie wielomianów
co daje nam złożoność czasową O(n
4)
Jeżeli chodzi o złożoność pamięciową to mamy tutaj n wielomianów
z których każdy ma swój wektor współczynników co daje O(n
2)
Gdy próbowałem się bawić i obliczać te wielomiany w Paincie to już dla n = 6
zajmowało to dość dużo czasu