Lausekkeella f(c_0, ..., c_n)=\int_{r_1}^{r_2} (g(x)-\sum_{k=1}^n c_k x^{n-k})^2\ dx
on yksikäsitteinen minimi, joka voidaan löytää etsimällä sellaiset polynomikertoimet c_k, että \frac{\partial I}{\partial c_i}=0 kaikilla i = 0, .., n.
Olkoon g(x)=e^x ja n = 2, 3, 4. Normaaliyhtälöt \frac{\partial f}{\partial c_i}=0 johtavat ehtoihin\sum_{k=1}^n c_k\mathop{\Bigg/}_{r_1}^{r_2} \frac{r^{2n-k-j+1}}{2n-k-j+1}=\int_{r_1}^{r_2} g(x)x^{n-j}\ dx,
josta kertoimet c_k voidaan ratkaista. Esimerkkikoodissa oikean puolen integraali määritetään numeerisesti metodin numerical_integral avulla.
reset()
from numpy import arange, polyval, zeros, linalg
x = var("x")
f(x) = e^x
interval = [-1, 1]
nmax = 4
data = zeros((nmax, nmax))
r1 = interval[0]
r2 = interval[1]
for n in range(2, nmax+1):
a = zeros((n, n))
b = zeros(n)
c = zeros(n)
for j in range(1, n+1):
for k in range(1, n+1):
a[j-1, k-1] = (r2^(2*n-k-j+1) - r1^(2*n-k-j+1))/(2*n-k-j+1)
b[j-1] = numerical_integral(f*x^(n-j), r1, r2)[0]
c = linalg.solve(a,b)
h = (r2-r1)/40
xs = arange(r1, r2+h, h)
y1 = [f(xi) for xi in xs]
y2 = polyval(c, xs)
erhe = abs(y2-y1)
maxer = max(erhe)
# Use trapezoidal rule to compute error
int1 = h*(sum(erhe) - 0.5*(erhe[0] + erhe[-1]))
int2 = h*(sum(erhe^2) - 0.5*(erhe[0]^2 + erhe[-1]^2))
# Plots
eplot = plot(f(x), (x, r1, r2), color="black")
polyplot = plot(polyval(c, x), (x, r1, r2), color="red", figsize=3)
epoly = eplot + polyplot
erheplot = plot(abs(f(x)-polyval(c, x)), (x, r1, r2), figsize=3)
# HTML code
html("<hr>$n=%s:$"%n)
html.table([["$%s$"%latex(matrix(a).n(digits=4)), "$%s$"%latex(vector(b).column().n(digits=4)), "$%s$"%latex(vector(c).column().n(digits=4))]], header=["$a$", "$b$", "$c$"])
html("$\\text{Abs. error = } %s\qquad\qquad\\text{L2 error = } %s$"%(maxer, int2))
html.table([["$\\text{Approksimaatio (n = %s)}$"%n,"$\\text{Virhe (n = %s)}$"%n], [epoly, erheplot]], header=True)
|
