newtoncotes.py

# calculate coefficients of Newton-Cotes integral formula
#  for \int_a^b f(x)dx
#  using sampling points x0(=a), x1, ..., xn-1, xn(=b)

import sympy as sym

n = 7;

x = sym.Symbol('x');
a = sym.Symbol('a');
b = sym.Symbol('b');

h = (b - a) / n;

lx = [0] * (n+1);

for i in range(0,n+1):
	lx[i] = a + i * h;

for i in range(0,n+1):
	tmp = 1
	for j in range(0,n+1):
		if j != i:
			tmp *= (x - lx[j]);
			tmp /= (lx[i] - lx[j]);
	result = sym.simplify(sym.integrate(tmp, (x, a, b)));
	print("(b - a) * " + str(sym.simplify(result / (b - a))) + " * f(x" + str(i) + ")");