import numpy as np import quadrature import polynomial_new as poly """ Settings: """ # Create an instance of the quadrature rule: x = np.array([-1.0, 0.0, 1.0]) w = np.array([1.0/3, 4.0/3, 1.0/3]) S = quadrature.QuadratureRule(x, w) # Boundary values of interval: a = -1.0 b = 1.0 # Number of test: ntest = 10 """ Function to integrate polynomial: """ def int_poly(P, a, b): # Coefficients of the primitive function: pc = P.coeff / np.arange(1, P.deg+2, dtype=float) # Evaluate powers at both boundary points: bx = np.array([b**i for i in range(1, P.deg+2)]) ax = np.array([a**i for i in range(1, P.deg+2)]) return (pc @ bx - pc @ ax) """ Main Test: """ for ii in range(ntest): # Generate random third order polynomial: #c = np.array([1.0, 0.0, 0.0, 1.0]) c = np.random.randn(4) print("Test %d: coeffs: "%ii, c) P = poly.polynomial(c) # Compute integral: I = int_poly(P, a, b) # Approximation with Simpson rule: ISimp = S.integrate(P) print("Test %d: diff = "%ii, np.abs(I - ISimp))