Skip to content
Snippets Groups Projects
Select Git revision
  • fda0cd082dcac7057795000a1917e348619088ce
  • main default protected
2 results

test_quad_transformed.py

Blame
  • nguyed99's avatar
    nguyed99 authored
    fda0cd08
    History
    Code owners
    Assign users and groups as approvers for specific file changes. Learn more.
    test_quad_transformed.py 1.11 KiB
    import numpy as np
    import quadrature
    import polynomial_new as poly
    
    """ Settings: """
    # Parameters of the quadrature rule on reference interval:
    x = np.array([-1.0, 0.0, 1.0])
    w = np.array([1.0/3, 4.0/3, 1.0/3])
    
    # Define linear transformation of [-1, 1] to [1, 5]:
    a = 1.0
    b = 5.0
    T = lambda x: 2 * x + 3
    dT = 2.0
    # Number of tests:
    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: """
    S = quadrature.LinearTransformed(x, w)
    
    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, T, dT)
        print("Test %d: diff = "%ii, np.abs(I - ISimp))