GNU Octave/Aproksymacja wielomianowa

Z Wikibooks, biblioteki wolnych podręczników.

[edytuj] Aproksymacja wielomianowa w przestrzeni Hilberta

Rozpatrzmy przestrzeń Hilberta funkcji całkowalnych z drugą potęgą L2(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