Skip to content

Instantly share code, notes, and snippets.

@tux21b
Created December 12, 2012 17:06
Show Gist options
  • Save tux21b/4269610 to your computer and use it in GitHub Desktop.
Save tux21b/4269610 to your computer and use it in GitHub Desktop.
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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment