I wrote my trapezoidal integration function to approximate the area under the curve for vehicle crash analysis. The function will take in two numpy vectors: (1) time vector and (2) acceleration or force. Thus, when I integrate, I’ll get velocity (if I put in acceleration vector) or total impulse exerted on the vehicle (if I put in the force).

def numerical_integration(x: np.array, fx: np.array, initial_value: float = 0.0):
    # error checking:
    if len(x) != len(fx):
        sys.exit("ERROR: x and y vectors are not the same length.")

    # array to keep track of integration value at each intervals:
    areas = np.zeros(len(fx), dtype = float)
    # we already have initial velocity so set that equal to initial area
    areas[0] = initial_value # total_area is used to keep track of total area under the curve
    total_area = initial_value
    
    # we'll interate for each value in x-data while keep the running total.
    for pos in range(1, len(x)):
        # 1: calculate the small (incremental) trapezoid area and add it to the total area
        # formula for trapezoid area is: width / 2 * (height1 + height2)
        total_area += 0.5 * (x[pos] - x[pos - 1]) * (fx[pos] + fx[pos - 1])
        # 2: assign the newly updated integration value to the areas array:
        areas[pos] = total_area
    
    # return the newly integrated data points:
    return areas, total_area

The problem is when I compare my function to the numpy’s trapezoidal function, np.trapz(), I get different results – very similar but noticeably off. If i run the following code in my main():

######################## START: IMPULSE AND MOMENTUM ANALYSIS ########################
    # TODO 1: Determine the momentum of the vehicle prior to impact, p1 , for each case
    p0_1 = mass * v0_1
    p0_2 = mass * v0_2
    
    # TODO 2: Determine the momentum of the vehicle after impact, p2 , for each case
    # plot_data(times_1, mass * velocities_1, times_2, mass * velocities_2, "MOMENTUM VS TIME", 'Time ($s$)', 'Momentum ($N-s$)', 'Case 1', 'Case 2')
    pf_1 = mass * velocities_1[-1]
    pf_2 = mass * velocities_2[-1]

    # TODO 3: Determine the total impulse imparted on the vehicle as a result of the collision for both cases 
            # – you should use the numerical integration algorithm discussed in class.
    # total impulse is the area under the curve of Force VS Time
    _ , impulse_1 = numerical_integration(times_1, forces_1, forces_1[0]) # FIXME
    _ , impulse_2 = numerical_integration(times_2, forces_2, forces_2[0]) # FIXME

    print('-' * 55)
    print(f"Momentum before impact: CASE 1: {p0_1}, CASE 2: {p0_2}")
    print(f"Momentum after  impact: CASE 1: {pf_1}, CASE 2: {pf_2}")
    print(f"nDelta P (not integration):nCase 1: {pf_1-p0_1}, Case 2: {pf_2-p0_2}")
    print(f"nTotal Impulse Case 1: {impulse_1}nTotal Impulse Case 2: {impulse_2}")
    print(f"nNumpy's Integration:nCase 1: {np.trapz(forces_1,times_1)}, Case 2: {np.trapz(forces_2, times_2)}")
    print('-' * 55)
######################## END: IMPULSE AND MOMENTUM ANALYSIS ########################

I get following output:

Momentum before impact: CASE 1: 0.7200000000000001, CASE 2: 0.6400000000000001
Momentum after  impact: CASE 1: -0.37721155161600284, CASE 2: -0.050379520229791856

Delta P (not integration):
Case 1: -1.0972115516160028, Case 2: -0.690379520229792

Total Impulse Case 1: -1.097211551616002
Total Impulse Case 2: -0.6903794111388608

Numpy's Integration:
Case 1: -1.0972115516160028, Case 2: -0.6903795202297909

Why such discrepancy between my trapezoidal and numpy’s trapezoidal integration? And why is numpy’s output closer to the actual value than mine? What am i doing wrong in my numerical integration function?