Source code for cubic_spline_integration

[docs]def x2_to_x3_integral(
        x1,
        x2,
        x3,
        x4,
        f1,
        f2,
        f3,
        f4,
):
    h = x3 - x2
    hm = x2 - x1
    hp = x4 - x3

    q1 = -hm / h
    q4 = 1 + hp / h

    avg_f = (
                               6 * q1 ** 3 * q4 ** 2 * (f2 + f3)
                               - 4 * q1 ** 3 * q4 * (2 * f2 + f3)
                               + 2 * q1 ** 3 * (f2 - f4)
                               - 6 * q1 ** 2 * q4 ** 3 * (f2 + f3)
                               + 3 * q1 ** 2 * q4 * (3 * f2 + f3)
                               + 3 * q1 ** 2 * (-f2 + f4)
                               + 4 * q1 * q4 ** 3 * (2 * f2 + f3)
                               - 3 * q1 * q4 ** 2 * (3 * f2 + f3)
                               + q1 * (f2 - f4)
                               + 2 * q4 ** 3 * (f1 - f2)
                               + 3 * q4 ** 2 * (-f1 + f2)
                               + q4 * (f1 - f2)
                       ) / (
                               12 * q1 * q4 * (q1 - 1) * (q1 - q4) * (q4 - 1)
                       )

    return avg_f * h




[docs]def f(x):
    return (
            -0.2 * x ** 3
            + x ** 2
            - x
            + 1
    )




[docs]a = 1


[docs]b = 3


from scipy import integrate


[docs]exact = integrate.quad(
    f,
    a,
    b,
)[0]

print(f"exact: {exact}")


[docs]x1 = a - 1


[docs]x2 = a


[docs]x3 = b


[docs]x4 = b + 1


[docs]f1 = f(x1)


[docs]f2 = f(x2)


[docs]f3 = f(x3)


[docs]f4 = f(x4)



[docs]approx = x2_to_x3_integral(
    x1,
    x2,
    x3,
    x4,
    f1,
    f2,
    f3,
    f4,
)

print(f"approx: {approx}")