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>=\underset{0}{\overset{1}{\int }}fg$. Dana jest $f\in L^2(0,1)$. Znaleźć wielomian $𝒲=\sum _{i=0}^N{\alpha }_i{x}^i$ stopnia $N>1$ najlepiej przybliżający $f$, czyli taki, który minimalizuje normę $\Vert f-𝒲{\Vert }_{L^2}$.

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

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

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

$f_i=\underset{0}{\overset{1}{\int }}f(x)x^{i}dx$

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

$G_{i,j}=<x^i,x^j>=\underset{0}{\overset{1}{\int }}{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 $\bar{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 $\bar{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) $\bar{\alpha }=G^{-1}\cdot \bar{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