GNU Octave/Aproksymacja wielomianowa

Aproksymacja wielomianowa w przestrzeni Hilberta

Rozpatrzmy przestrzeń Hilberta funkcji całkowalnych z drugą potęgą $L^{2}(0,1)$, z iloczynem skalarnym $<f,g>=\int _{0}^{1}{fg}$. Dana jest $f\in L^{2}(0,1)$. Znaleźć wielomian ${\mathcal {W}}=\sum _{i=0}^{N}{\alpha _{i}x^{i}}$ stopnia $N>1$ najlepiej przybliżający $f$, czyli taki, który minimalizuje normę $\|f-{\mathcal {W}}\|_{L^{2}}$.

W tym celu należy rozwiązać układ równań:

$G(1,x,x^{2},...,x^{N})\cdot {\overline {\alpha }}={\overline {f}},$

gdzie ${\overline {\alpha }}$ jest szukanym wektorem współczyników ${\mathcal {W}}$ macierz G jest macierzą Grama, a ${\overline {f}}$ jest wektorem:

$f_{i}=\int _{0}^{1}{f(x)x^{i}dx}$

W naszym przypadku macierz Grama G sprowadza się do macierzy Hilberta, gdyż

$G_{i,j}=<x^{i},x^{j}>=\int _{0}^{1}{x^{i+j}}={\frac {1}{i+j+1}}$

Algorytm ten można wykonać w środowisku Octave następująco. Zdefiniujmy funkcję charakterystyczną $\chi _{[0,1/2]}\!$:

function [y]=fcharakt(x)
   y=(x<=0.5);
endfunction;

Do obliczenia wektora ${\overline {f}}$ potrzebujemy funkcji $f\cdot {}x^{n}$:

function [y]=f_razy_xn(x)
   global n;
   global fun;
   z=feval(fun,x);
   y=(z).*(x^n);
endfunction;

Dla zadanego stopnia wielomianu N tworzymy macierz Grama:

N=15;
G=hilb(N+1);

Tworzymy wektor ${\overline {f}}$:

f=zeros(1,N+1);
global n;
global fun;
fun="fcharakt";
for (n=0:N)
   f(n+1)=quad("f_razy_xn", 0, 1);
endfor;

I wreszcie obliczamy współczynniki (w odwrotnej kolejności) ${\overline {\alpha }}=G^{-1}\cdot {}{\overline {f}}$:

alfa=G\f'

Uwaga. Ta metoda nie jest dobra dla dużych N, gdyż macierz Hilberta jest źle uwarunkowana.

Stopień wielomianu: 1

Stopień wielomianu: 3

Stopień wielomianu: 5

Stopień wielomianu: 10

Stopień wielomianu: 20