Numerical and symbolic integral computations using NumPy/SciPy unexpectedly disagree

Question

Problem

I'm attempting to write some Python code to compute the value of a definite integral that describes a well-understood physical system. I know enough about this physical system so as to be able to approximate the result for some specific values of the integrand's constants and some sets of bounds, and can thereby run a sort of sanity check on any code I write.

For a variety of reasons (one of which is execution time), I'd strongly prefer to compute the result of the integration by evaluating the symbolic antiderivative of the integrand, as opposed to performing some numerical computation. I've determined the form of this symbolic antiderivative, and have written some Python code to compare the result of the integration using this symbolic computation to the result obtained via numerical integration and also the approximate expected result. Here's a minimum viable example of that code, which compares these results for a set of specific values of the constants and bounds:

import numpy as np
from scipy.integrate import quad

# Example bound values
z1 = -500.7
z0 = -1000.4

# Example constant values
A = 1.78
B = 0.454
C = 0.0132

theta = np.pi/4
beta = (A - B * np.exp(C * z0)) * np.cos(theta)
K = C * np.sqrt(A**2 - beta**2)/beta

# Symbolic elementary antiderivative
def integral_result(z):
        return (A / (B * C * np.abs(K))) * np.sqrt((A - beta) * (A + beta) * (C**2 + K**2)) * np.arccosh((B * np.exp(C * z) + A) / beta)

# The exact integrand
def integrand(z):
        return A*np.sqrt(((A - beta)*(A + beta)*(C**2 + K**2))/((K**2*(A**2 + 2*A*B*np.exp(C*z) + B**2*np.exp(2*C*z) - beta**2))))

result, error = quad(integrand, z0, z1)

print("Approximate Expected Value:", (z1 - z0))
print("Numerical Result:", result)
print("Analytical Result:", (integral_result(z1) - integral_result(z0)))

This code produces the following output:

Approximate Expected Value: 499.7
Numerical Result: 1257.7634049781054
Analytical Result: 0.2566028102053224

The numerical result has low estimated error, and is roughly as close to the approximate expected value as I had anticipated (they will differ, perhaps by a lot, but should at least be of the same or neighboring orders of magnitude). The result determined using the symbolic elementary antiderivative disagrees with the result determined using numerical methods, and differs from the approximate expected value by several orders of magnitude.

What I've Tried

I have run this by several members of the math and physics faculty at my institution, and I am confident that I have not erred in my determination of the symbolic elementary antiderivative, and also that it is valid for all values of the constant and bounds for which I am concerned. The result obtained via numerical methods also agrees well with expectation and also with several pieces of old code designed to solve a similar problem. It is clear, then, that my issue is not a mathematical one, but a coding issue.

I've tried all of the obvious things, including:

1. Term-by-Term Breakdown: Breaking the computation down into each of its terms, printing their values, and verifying that they make sense (e.g. checking the argument and output of np.arccosh).
2. Parameter Variation Testing: Comparing the results for a range of A,B,C and theta values.
3. Floating Point Precision and Round-Off Errors: I've investigated whether the issue could be related to floating-point precision or round-off errors. Given the complex nature of some of the functions, particularly those involving logarithms and square roots, I checked if small numerical inaccuracies in intermediate steps might be leading to larger discrepancies in the final result.
4. Mathematical Analysis: Various things to ensure that I haven't erred in my mathematics, which I won't repeat here as they're most probably entirely irrelevant to what is almost surely a coding issue.

I've talked this over with a number of different people, and I unfortunately still don't have the faintest idea why these results disagree. I suspect that there is some floating point precision issue, something I've somehow missed in the documentation of some function, or similar that hasn't occurred to me.

Any insights or suggestions would be greatly appreciated!

---

Addendum

The symbolic elementary antiderivative is given as follows (and is due not to me but to the very helpful folks over at Math StackExchange):

I'm using Python 3.9.6, NumPy 1.26.0, and SciPy 1.11.2, if that is an any way useful.

Answer 1

Jared has provided a great solution to the question asked.

However, the linked question (after edits) on Maths Stack Exchange (https://math.stackexchange.com/questions/4968907/computing-the-antiderivative-of-frac1-sqrta22-a-b-ec-zb2-e2cz-b/4969007#4969007) shows the integral to be slightly different. It is the square root of what was asked here. Moreover, the Maths Stack Exchange answer is wrong for either integral.

For completeness I'll suggest a different solution and check it by code for the (slightly different) integral on Maths Stack Exchange.

Following the initial substitution of Maths Stack Exchange, but NOT its subsequent implementation, let

whilst

Hence, the integral becomes (different from Maths Stack Exchange):

By analogy by what you would do for trigonometric, rather than hyperbolic functions, of this form (https://en.wikipedia.org/wiki/Tangent_half-angle_substitution), substitute

Then,

The integral above becomes

Combining this with your original gives

This can be simplified since

leads to

Final solution (minus a constant of integration):

Code:

import numpy as np
from scipy.integrate import quad

# Example bound values
z1 = -500.7
z0 = -1000.4

# Example constant values
A = 1.78
B = 0.454
C = 0.0132

theta = np.pi / 4
beta = ( A - B * np.exp( C * z0 ) ) * np.cos( theta )
K = C * np.sqrt( A ** 2 - beta ** 2 ) / beta

def integral_result( z ):
    t = np.sqrt( ( A + B * np.exp( C * z ) - beta ) / ( A + B * np.exp( C * z ) + beta ) )
    alpha = np.sqrt( ( A - beta ) / ( A + beta ) )
    return A * np.sqrt( ( C ** 2 + K ** 2 ) / ( C ** 2 * K ** 2 ) ) * np.log( np.abs( ( t - alpha ) / ( t + alpha ) ) )

def integrand(z):
    return 1 / np.sqrt( A ** 2 + 2 * A * B *np.exp( C * z ) + B ** 2 * np.exp( 2 * C * z ) - beta ** 2 )

result, _ = quad( integrand, z0, z1 )
r = A * np.sqrt( ( A ** 2 - beta ** 2 ) * ( C ** 2 + K ** 2 ) ) / np.abs( K )
result *= r

print( "Numerical Result:", result )
print( "Analytical Result:", ( integral_result( z1 ) - integral_result( z0 ) ) )

Output:

Numerical Result: 1257.7634049781054
Analytical Result: 1257.7634050261993