Numerical and symbolic integral computations using NumPy/SciPy unexpectedly disagree

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.

I took your integral and solved it using sympy and got a different answer than the one you got from Math Exchange (I also expanded K using the definition you show in your code).

import sympy as smp

z, A, B, C, beta = smp.symbols("z A B C \\beta", real=True)
K = C*smp.sqrt(A**2-beta**2)/beta

pre = A*smp.sqrt((A-beta)*(A+beta)*(C**2+K**2))/abs(K)

integrand = 1/(A**2 + 2*A*B*smp.exp(C*z) + B**2*smp.exp(2*C*z) - beta**2)
res = smp.integrate(integrand, z)

combined = (pre*res).expand().simplify()

With that I get the following result:

I've input this into your code. I've also fixed your numerical integral. You had a K**2 instead of abs(K). I split up the numerical integration (since there is no reason to integrate a constant).

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):
        Abeta2 = A**2 - beta**2
        term1 = A*np.sqrt(Abeta2/beta**2)
        term2 = 2*C*beta*z*Abeta2 + (A-beta)*Abeta2*np.log((A + B*np.exp(C*z) + beta)/B)
        term3 = (A+beta)*Abeta2*np.log((A + B*np.exp(C*z) - beta)/B)
        term4 = np.abs(A*beta/np.sqrt(Abeta2))
        numer = term1*(term2 - term3)*term4
        denom = 2*C*beta*Abeta2**2
        return numer/denom

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

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

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

Now the numerical and analytical integrals match.

Approximate Expected Value: 499.7
Numerical Result: 999.1911207875748
Analytical Result: 999.1911207875743