Analytical solution suffering from lack of precision

Question

I'm attempting to use python to solve the following analytical solution (eq. 4 from here) for a heat transport problem to verify the results calculated by a numerical model:

The small reproducible example below invokes the quad function in scipy.integrate to solve the integral. So far, so good. The result from the integration is then multiplied by what I'm loosely referring to as a prefix term to give the final "Delta T" (change in temperature) result - ultimately what I need.

Unfortunately when I multiply the integral result, a monotonically increasing value, by the prefix term, a monotonically decreasing value, I don't get the expected smooth solution. I'm wondering if python offers tricks for forcing greater precision in the multiplication operation (i.e., result = prefix * integral[0])?

The example script below ends with a 3-part plot highlighting what I think is the problem. That is, the left-most plot shows the monotonicity of the prefix term, the middle plot shows the monotonicity of the integral result, and the right most plot shows that multiplying the increasingly larger integrand values by ever-decreasing prefix values results in a non-smooth result - a result from a lack of precision possibly? Is there a fix for this?

x, y, and t values are related to space and time. I've chosen an arbitrarily small set of values for this small reproducible example. The larger T gets, the more "unsmooth" the results look.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.integrate import quad
import math


X = [val + 0.5 for val in np.arange(20)]
Y = [3.59]
days2seconds = 86400
times = [200, 400]  # days
v = 9.553639342343923e-06  # m/s (velocity)
D = 2.9059885319758944e-06  # m^2/s (diffusivity)
D_prime = 1.4150943396226415e-06  # m^2/s (diffusivity overburden)
h_prime = 0.8372827804107424  # unitless (heat capacity ratio)
H = 100.0  # m (aquifer thickness)
T0 = 80.0  # deg C (initial temperature)
T1 = 30.0  # deg C (injected temperature)


# Some functions for calculating the Barends analytical solution
def barends_eqn4(sigma, x, y, t):
    exponent = -1 * sigma ** 2 - ((x * v) / (4 * D * sigma)) ** 2
    term1 = math.exp(exponent)
    term2 = x ** 2 * h_prime * math.sqrt(D_prime) / (8 * D * H * sigma ** 2)
    term3 = y / (2 * math.sqrt(D_prime))
    term4 = t - x ** 2 / (4 * D * sigma ** 2)
    
    # put it all together
    eqn4_val = term1 * math.erfc((term2 + term3) * (term4) ** (-0.5))
    
    return eqn4_val


def calc_analytical_sln(times):
    times = [tm * days2seconds for tm in times]
    # 
    analytical_sln = np.zeros((len(times), len(X)))
    integral_rslt = np.zeros_like(analytical_sln)
    prefix_rslt = np.zeros_like(analytical_sln)
    for t_idx, t in enumerate(times):
        for k, y in enumerate(Y):  # 1 row
            for j, x in enumerate(X):  # columns 
                lower_lim = x / (2 * math.sqrt(D * t))
                integral = quad(barends_eqn4, lower_lim, np.inf, args=(x, y, t))
                integral_rslt[t_idx, j] = integral[0]
                
                # multiply the prefix by the solution to the integral
                prefix = (2 * (T1 - T0)) / (math.sqrt(np.pi)) * math.exp((x * v) / (2 * D))
                prefix_rslt[t_idx, j] = prefix
                result = prefix * integral[0]
                # store result for plotting
                analytical_sln[t_idx, j] = result + T0
    
    return analytical_sln, integral_rslt, prefix_rslt


analytical_answers, integral_rslt, prefix_rslt = calc_analytical_sln(times)

# Because the values in prefix_rslt are negative, need to flip the sign
# for a log-scaled plot
fig, (ax1, ax2, ax3) = plt.subplots(1, 3)
fig.set_size_inches(12, 3)
ax1.plot(X, -prefix_rslt[0], color='r', label="prefix term")
ax1.set_yscale("log")
ax1.legend(loc='lower right')

ax2.plot(X, integral_rslt[0], color='b', label="integral result")
ax2.set_yscale("log")
ax2.legend(loc='lower left')

ax3.plot(X, analytical_answers[0], 'k', label="analytical solution (product)")
ax3.yaxis.tick_right()
ax3.set_yscale("log")
ax3.legend(loc='lower right')

plt.show()

Answer 1

You have very large exponential terms - a large positive one outside the integral and a large negative one inside it. This is not very good numerically.

All you have to do is take the exponential part from outside the integral and move it inside. So, inside the integrand:

exponent = -sigma ** 2 - ((x * v) / (4 * D * sigma)) ** 2 + (x * v) / (2 * D)

or, equivalently and better,

exponent = -(sigma - x * v / ( 4 * D * sigma ) ) ** 2

Also remember to remove it from your "prefix" term:

prefix = (2 * (T1 - T0)) / (math.sqrt(np.pi))