GNU Octave/Aproksymacja wielomianowa

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

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

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

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

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

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

${\displaystyle 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ą ${\displaystyle \chi _{[0,1/2]}\!}$ :

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

Do obliczenia wektora ${\displaystyle {\overline {f}}}$  potrzebujemy funkcji ${\displaystyle 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 ${\displaystyle {\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) ${\displaystyle {\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.