Beispiel 22

def dreitermrek(s, n):
    if n == 0:
        return 1
    r1 = dreitermrek(s, n-1)
    a = s(x*r1, r1) / s(r1, r1)
    if n == 1:
        return x - a
    r2 = dreitermrek(s, n-2)
    b = s(r1, r1) / s(r2, r2)
    return ((x - a) * r1 - b * r2).expand()

def gauss(n, w, f, a, b):
    s = lambda u, v: integrate(w(x) * u * v, x, a, b)
    p = dreitermrek(s, n+1)
    print 'p:', p
    xi = [i[0] for i in RR[x](p).roots()]
    wi = [integrate(w(x) * prod([(x - xi[j]) / (xi[i] - xi[j]) \
        for j in range(n+1) if i != j]), x, a, b).n() \
        for i in range(n+1)]
    for i in range(n+1):
        print 'x[%d] = %s, w[%d] = %s' % (i, xi[i], i, wi[i])
    result = sum([wi[i] * f(xi[i]) for i in range(n+1)])
    print 'gauss(%s) = %s' % (f(x)*w(x), result)
    return result

exact = -0.239811742000565

f(x) = ln(x) * sin(x)
w(x) = 1
g = gauss(3, w, f, 0, 1)
print 'error:', abs(g - exact).n()
print 

f(x) = sin(x)
w(x) = ln(x)
g = gauss(3, w, f, 0, 1)
print 'error:', abs(g - exact).n()
print