Аналитическое решение, страдающее от недостатка точности

Я пытаюсь использовать Python для решения следующего аналитического решения (уравнение 4 отсюда) для задачи переноса тепла, чтобы проверить результаты, рассчитанные с помощью числовой модели:

Небольшой воспроизводимый пример ниже вызывает функцию quad в scipy.integrate для решения интеграла. Все идет нормально. Затем результат интегрирования умножается на то, что я условно называю префиксным термином, чтобы получить окончательный результат «Дельта Т» (изменение температуры) - в конечном итоге то, что мне нужно.

К сожалению, когда я умножаю интегральный результат, монотонно возрастающее значение, на префиксный член, монотонно убывающее значение, я не получаю ожидаемого гладкого решения. Мне интересно, предлагает ли Python способы повышения точности операции умножения (т. е. result = prefix * integral[0])?

Приведенный ниже пример сценария заканчивается сюжетом из трех частей, показывающим, в чем, по моему мнению, заключается проблема. То есть самый левый график показывает монотонность префиксного термина, средний график показывает монотонность интегрального результата, а самый правый график показывает, что умножение все больших значений подынтегрального выражения на постоянно уменьшающиеся значения префикса приводит к неудовлетворительному результату. -гладкий результат - возможно, результат недостаточной точности? Это можно как-то исправить?

Значения x, y и t связаны с пространством и временем. Для этого небольшого воспроизводимого примера я выбрал произвольно небольшой набор значений. Чем больше T, тем более «негладкими» выглядят результаты.

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()