Численное интегрирование на разностных схемах#
Численное интегрирование на разностных схемах#
Опишем методику численного расчета интеграла \(\int_a^b f(x) d x\) с помощью разностных схем в старших порядках. В отличие от квадратур, тут нам будут нужны производные.
Теорема
Пусть функция \(f(x)\) аналитична в каждой точке отрезка \([a ; b]\). Тогда её интеграл можно вычислять квадратурой на разностных схемах старшего порядка по формуле:
\[
\begin{gathered}
\int_a^b f(x) d x=h \sum_{j=0}^{J-1} \sum_{k=-m}^m W_{m-k}^m y_{j+k}+O\left(h^{2 m+2}\right), \quad y_j = f\left(a + (j + 0.5) \cdot h\right)
\end{gathered}
\]
где коэффициенты не зависят от функции и задаются рекурсивными формулами:
\[
W_k^m=\sum_{n=0}^m \frac{A_{k, 2 n}^{m-n}}{2^{2 n}(2 n+1) !}
\]
\[
A_{k n}^m=\sum_{l=0}^m(-1)^{k-m} D_n^{(l)} C_{n+2 l}^{k-m+l}
\]
\[
D_n^{(j+1)}=\sum_{k=0}^{n+2 j} \sum_{l=0}^j(-1)^k D_n^{(l)} C_{n+2 l}^{k-j+l} \frac{(n+2 j-2 k)^{n+2 j+2}}{2^{n+2 j+2}(n+2 j+2) !}, \quad D_n^{(0)} \equiv 1
\]
Примечание
Обратите внимание, что квадратура вылезает за сетку на отрезке \([a ; b]\)! Это специфика этого метода. Также он использует нестандартную сетку.
👉Доказательство
Пусть функция \(f(x)\) аналитична в каждой точке отрезка \([a ; b]\). Тогда она аналитична хотя бы немного за пределами отрезка \([a ; b]\), т. е. на некотором интервале \(\left(a^{\prime} ; b^{\prime}\right)\supset[a ; b]\). Это условие нам понадобится для того чтобы использовать только симметричные разностные схемы даже на границах отрезка.
Разобьем отрезок \([a ; b]\) на \(J\) одинаковых кусочков и на каждом кусочке разложим \(f(x)\) в ряд Тейлора в средней точке каждого кусочка, после чего вычислим интеграл на каждом кусочке. Имеем,
\[\begin{split}
\begin{gathered}
\int_a^b f(x) d x=\sum_{j=0}^{J-1} \int_{x_j-\Delta x}^{x_j+\Delta x} f(x) d x= \\
=\sum_{j=0}^{J-1} \sum_{n=0}^m \frac{f^{(2 n)}\left(x_j\right)}{(2 n) !} \int_{x_j-\Delta x}^{x_j+\Delta x}\left(x-x_j\right)^{2 n} d x+O\left(\Delta x^{2 m+2}\right)= \\
=\sum_{j=0}^{J-1} \sum_{n=0}^m \frac{f^{(2 n)}\left(x_j\right)}{(2 n+1) !} 2 \Delta x^{2 n+1}+O\left(\Delta x^{2 m+2}\right) .
\end{gathered}
\end{split}\]
Где \(\Delta x=(b-a) / 2 J=\cfrac{h}{2}\)- половина длины кусочка, \(J-\) число кусочков, \(x_j=a+(2 j+1) \Delta x-\) средняя точка \(j\) го кусочка, и учтено, что \(\sum_{j=0}^{J-1} O(\Delta x)=O(1)\) и что интеграл в симметричных пределах от нечетных степеней \(\left(x-x_j\right)^{2 n+1}\) равен нулю.
Итого получили метод численного расчёта интеграла произвольного порядка сходимости \(2m+2\), в котором нам нужны не только значения функции в узлах сетки, но и значения всех производных до \(2m\) включительно в этих же узлах с нужной точностью.
В следующем году в семинаре 4.1 будет выведена общая формула для симметричной разностной схемы произвольного порядка аппроксимации. Пока её вывод можно посмотреть в файле numer.pdf в доп. материалах. Сама формула:
\[
f^{(n)}\left(x_j\right)=\frac{1}{2^n h^n} \sum_{k=0}^{n+2 m} A_{k n}^m \cdot y_{j+n+2 m-2 k}+O\left(h^{2 m+2}\right)
\]
\[
A_{k n}^m=\sum_{l=0}^m(-1)^{k-m} D_n^{(l)} C_{n+2 l}^{k-m+l}
\]
\[
D_n^{(j+1)}=\sum_{k=0}^{n+2 j} \sum_{l=0}^j(-1)^k D_n^{(l)} C_{n+2 l}^{k-j+l} \frac{(n+2 j-2 k)^{n+2 j+2}}{2^{n+2 j+2}(n+2 j+2) !}, \quad D_n^{(0)} \equiv 1
\]
Отметим, что это разностная схема не сходится с ростом \(m\)! при фиксированном шаге. Поэтому рекомендуется брать \(m\) меньше, чем точность исходных данных. Для машинной точности (аналитические функции, 16 значящих знаков) хорошее значение \(m=7\).
Нам нужны только чётные производные и даже немаксимального порядка аппроксимации (т.к. для каждой из них производные домножаются на \(\Delta x = \frac{h}{2}\) в каких-то степенях). Внесём нужные модификации в формулу:
\[
f^{(2 n)}\left(x_j\right)=\frac{1}{(2 \Delta x)^{2 n}} \sum_{k=0}^{2 m} A_{k, 2 n}^{m-n} \cdot y_{j+m-k}+O\left(\Delta x^{2 m-2 n+2}\right)
\]
Подставляя эту формулу в оценку интеграла, имеем:
\[
\int_a^b f(x) d x=2 \Delta x \sum_{j=0}^{J-1} \sum_{n=0}^m \sum_{k=0}^{2 m} \frac{A_{k, 2 n}^{m-n}}{2^{2 n}(2 n+1) !} \cdot y_{j+m-k}+O\left(\Delta x^{2 m+2}\right)
\]
Вводя коэфициенты
\[
W_k^m=\sum_{n=0}^m \frac{A_{k, 2 n}^{m-n}}{2^{2 n}(2 n+1) !}
\]
получаем квадратуру, основанную на разностных схемах на \(2m\) точках (заметье, что она вылезает за отрезок интегрирования!):
\[\begin{split}
\begin{gathered}
\int_a^b f(x) d x&=2 \Delta x \sum_{j=0}^{J-1} \sum_{k=-m}^m W_{m-k}^m y_{j+k}+O\left(\Delta x^{2 m+2}\right)= \\
&=h \sum_{j=0}^{J-1} \sum_{k=-m}^m W_{m-k}^m y_{j+k}+O\left(h^{2 m+2}\right)
\end{gathered}
\end{split}\]
Важно заметить, что квадратура для интеграла «вылезает» за пределы области интегрирования - отрезка \([a ; b]\). Это специфика применяемых здесь разностных схем в старших порядках. Но зато интеграл вычисляется очень быстро, \(J+2 m\) обращений к функции дают ошибку \(O\left(\Delta x^{2 m+2}\right)\). Так, при \(m=7\) вычисление сетки значений \(\exp \left(x_j\right)\) для интерграла \(\int_a^b e^x d x\) машинная точность \(10^{-16}\) достигается при длине кусочка \(x_{j+1}-x_j=0.25\), тогда как обычный метод Симпсона (с ошибкой \(O\left(\Delta x^4\right)\) ) для достижения машинной точности требует шаг \(\Delta x=0.0001\), что на три порядка меньше!
Также несложно показать, что с ростом \(m\) данная схема расходится, поэтому важно его опдбирать исходя из погрешностей исходных данных.
Более того, для данной квадратуры естественно можно применять алгоритм Рунге для оценки практической погрешности.
Что делать, если вылезать за отрезок нельзя?#
Отметим, что если выход сетки \(\left\{x_j\right\}\) за пределы области интегрирования недопустим (неаналитичность или бесконечные пределы), можно, как один из вариантов, вместо интеграла \(\int_a^b f(x) d x\) вычислять интеграл \(\int_0^1 f(x(t)) \dot{x}(t) d t\) с заменами, например,
\[x(t)=a+(b-a) \cdot \frac{1-\cos (\pi t)}{2}\]
или, при \(b=+\infty\),
\[
x(t)=\frac{2 a}{(1+\cos (\pi t))}
\]
Аналитичность функции \(f(x(t)) \dot{x}(t)\) в каждой точке отрезка \(t \in[0 ; 1]\) должна гарантироваться заменой в обязательном порядке. В противном случае нужно использовать другую замену.
Примеры#
Show code cell source
from decimal import Decimal, getcontext
import math
from functools import lru_cache
###############################################################################
# 1) Вспомогательные функции для разностных квадратур #
###############################################################################
getcontext().prec = 100 # повышенная точность для Decimal
@lru_cache(None)
def factorial_dec(n: int) -> Decimal:
if n < 0:
raise ValueError("Negative factorial not defined")
if n == 0 or n == 1:
return Decimal(1)
return Decimal(n) * factorial_dec(n - 1)
@lru_cache(None)
def comb_dec(n: int, k: int) -> Decimal:
if k < 0 or k > n:
return Decimal(0)
return factorial_dec(n) / (factorial_dec(k) * factorial_dec(n - k))
@lru_cache(None)
def D(n: int, j: int) -> Decimal:
"""
Вычисляет D_n^(j) по заданной рекуррентной формуле (см. теорию).
"""
if j == 0:
return Decimal(1)
# "прогреваем" кэш:
_ = D(n, j-1)
return _compute_all_D_up_to(n, j)[j]
@lru_cache(None)
def _compute_all_D_up_to(n: int, J: int):
results = [Decimal(0)]*(J+1)
results[0] = Decimal(1)
for j in range(J):
s = Decimal(0)
for k in range(n + 2*j + 1):
for l in range(j+1):
val = (
Decimal((-1)**k)
* results[l]
* comb_dec(n+2*l, k - (j - l))
* Decimal((n + 2*j - 2*k)**(n + 2*j + 2))
/ (Decimal(2)**(n + 2*j + 2) * factorial_dec(n + 2*j + 2))
)
s += val
# Финальная формула берёт только s, без + results[j].
results[j+1] = s
return results
@lru_cache(None)
def A(k: int, n: int, m: int) -> Decimal:
"""
A_{k,n}^m = sum_{l=0}^m [ (-1)^(k-m) * D_n^(l) * C_{n+2l}^{k - m + l} ].
"""
s = Decimal(0)
sign = Decimal((-1)**(k - m))
for l in range(m+1):
s += (
sign
* D(n, l)
* comb_dec(n+2*l, k - m + l)
)
return s
@lru_cache(None)
def W_list(m: int):
"""
Коэффициенты W_k^m (k=0..2m) в Decimal.
W_k^m = sum_{n=0}^m [ A_{k,2n}^{m-n} / (2^(2n)*(2n+1)!) ].
"""
result = []
for k in range(2*m + 1):
s = Decimal(0)
for n_ in range(m+1):
a_val = A(k, 2*n_, m - n_)
denom = (Decimal(2)**(2*n_) * factorial_dec(2*n_+1))
s += a_val / denom
result.append(s)
return result
def build_function_values_diff_scheme(f, a: float, b: float, J: int, m: int):
"""
Центральная сетка: x_j = a + (j+0.5)*h, j=-m..(J-1+m).
Возвращает (y_values, h, calls).
"""
h = (b - a)/J
values = []
for j_real in range(-m, (J - 1) + m + 1):
x_val = a + (j_real + 0.5)*h
values.append(f(x_val))
calls = len(values)
return values, h, calls
def integrate_by_diff_scheme(f, a: float, b: float, J: int, m: int):
"""
Разностная (дифференсная) квадратура порядка 2m:
∫[a..b] f(x) dx ≈ h Σ_j Σ_k [ W_{m-k} * f( x_{j+k} ) ].
Возвращает (approx, calls).
"""
y_values, h, calls = build_function_values_diff_scheme(f, a, b, J, m)
W_dec = W_list(m)
W = [float(wd) for wd in W_dec]
offset = m
total_sum = 0.0
for j in range(J):
loc_sum = 0.0
for k in range(-m, m+1):
idx = (j + k) + offset
loc_sum += W[m - k] * y_values[idx]
total_sum += loc_sum
return (h * total_sum, calls)
###############################################################################
# 2) Метод Симпсона (классический) с заданным числом вызовов f #
###############################################################################
def build_function_values_simpson(f, a: float, b: float, n: int):
h = (b - a)/n
y = [f(a + i*h) for i in range(n+1)]
return y, h
def simpson_by_calls(f, a: float, b: float, calls: int):
"""
Метод Симпсона «под число вызовов calls».
Если (calls-1) нечётно, уменьшаем n на 1.
Если n < 2, берём n=2.
Возвращает (approx, real_calls).
"""
n = calls - 1
if n % 2 == 1:
n -= 1
if n < 2:
n = 2
y, h = build_function_values_simpson(f, a, b, n)
real_calls = n + 1
# Классическая формула Симпсона
s_odd = 0.0
s_even = 0.0
for i in range(1, n, 2):
s_odd += y[i]
for i in range(2, n, 2):
s_even += y[i]
simpson_approx = (h/3.0)*(y[0] + y[n] + 4.0*s_odd + 2.0*s_even)
return (simpson_approx, real_calls)
###############################################################################
# 3) Основная часть #
###############################################################################
def main():
# ------------------------------------------------------
# Пример 1: f1(x)= sin(3x)/(1 + x^2) на [0,1].
# Нет простого закрытого выражения; оставляем без аналитической формулы.
# ------------------------------------------------------
def f1(x: float) -> float:
return math.sin(3*x)/(1.0 + x*x)
# "Имитация" точного значения с помощью большого числа разбиений Симпсона.
# (Если хотите исключить численную оценку, уберите пример 1.)
def get_almost_exact_for_f1():
big_calls = 10_000_001
val, _ = simpson_by_calls(f1, 0.0, 1.0, big_calls)
return val
exact_val_1 = get_almost_exact_for_f1()
a1, b1 = 0.0, 1.0
print("\n=========== Пример 1 ===========")
print("Интеграл: I1 = ∫[0..1] sin(3x)/(1+x^2) dx")
print("(псевдо-«точное» значение, большое разбиение Симпсона)\n")
for m in [3, 4, 5, 6, 7]:
print(f"--- m={m} ---")
for J in [2, 4, 8, 16]:
diff_approx, diff_calls = integrate_by_diff_scheme(f1, a1, b1, J, m)
diff_err = abs(diff_approx - exact_val_1)
simpson_approx, simp_real_calls = simpson_by_calls(f1, a1, b1, diff_calls)
simpson_err = abs(simpson_approx - exact_val_1)
print(f" J={J}, calls={diff_calls:2d} | "
f"DiffScheme={diff_approx:.8f}, err={diff_err:.2e} | "
f"Simpson={simpson_approx:.8f}, err={simpson_err:.2e} "
f"(n={simp_real_calls - 1} subintrvls)")
print()
# ------------------------------------------------------
# Пример 2: f2(x)= e^(-x^2) на [-1,1].
# Аналитически: ∫ e^(-x^2) dx = нет в элементарных, но
# ∫_{-1}^1 e^(-x^2) dx = √π * erf(1).
# ------------------------------------------------------
def f2(x: float) -> float:
return math.exp(-x*x)
a2, b2 = -1.0, 1.0
# Точное аналитическое значение:
exact_val_2 = math.sqrt(math.pi) * math.erf(1)
print("\n=========== Пример 2 ===========")
print("Интеграл: I2 = ∫[-1..1] e^(-x^2) dx = √π * erf(1)\n")
for m in [3, 4, 5, 6, 7]:
print(f"--- m={m} ---")
for J in [2, 4, 8, 16]:
diff_approx, diff_calls = integrate_by_diff_scheme(f2, a2, b2, J, m)
diff_err = abs(diff_approx - exact_val_2)
simpson_approx, simp_real_calls = simpson_by_calls(f2, a2, b2, diff_calls)
simpson_err = abs(simpson_approx - exact_val_2)
print(f" J={J}, calls={diff_calls:2d} | "
f"DiffScheme={diff_approx:.8f}, err={diff_err:.2e} | "
f"Simpson={simpson_approx:.8f}, err={simpson_err:.2e} "
f"(n={simp_real_calls - 1} subintrvls)")
print()
# ------------------------------------------------------
# Пример 3: f3(x) = sin^2(100 x) на [-1,1].
# Точно: ∫_{-1}^1 sin^2(100x) dx = 1 - sin(200)/200.
# ------------------------------------------------------
def f3(x: float) -> float:
return math.sin(100*x)**2
a3, b3 = -1.0, 1.0
# Аналитическое значение:
exact_val_3 = 1.0 - math.sin(200.0)/200.0
print("\n=========== Пример 3 ===========")
print("Интеграл: I3 = ∫[-1..1] sin^2(100x) dx = 1 - sin(200)/200\n")
for m in [3, 4, 5, 6, 7]:
print(f"--- m={m} ---")
for J in [2, 4, 8, 16]:
diff_approx, diff_calls = integrate_by_diff_scheme(f3, a3, b3, J, m)
diff_err = abs(diff_approx - exact_val_3)
simpson_approx, simp_real_calls = simpson_by_calls(f3, a3, b3, diff_calls)
simpson_err = abs(simpson_approx - exact_val_3)
print(f" J={J}, calls={diff_calls:2d} | "
f"DiffScheme={diff_approx:.8f}, err={diff_err:.2e} | "
f"Simpson={simpson_approx:.8f}, err={simpson_err:.2e} "
f"(n={simp_real_calls - 1} subintrvls)")
print()
if __name__ == "__main__":
main()
=========== Пример 1 ===========
Интеграл: I1 = ∫[0..1] sin(3x)/(1+x^2) dx
(псевдо-«точное» значение, большое разбиение Симпсона)
--- m=3 ---
J=2, calls= 8 | DiffScheme=0.52003163, err=2.80e-03 | Simpson=0.51749322, err=2.58e-04 (n=6 subintrvls)
J=4, calls=10 | DiffScheme=0.51726196, err=2.64e-05 | Simpson=0.51731432, err=7.87e-05 (n=8 subintrvls)
J=8, calls=14 | DiffScheme=0.51723573, err=1.46e-07 | Simpson=0.51725080, err=1.52e-05 (n=12 subintrvls)
J=16, calls=22 | DiffScheme=0.51723559, err=6.32e-10 | Simpson=0.51723754, err=1.95e-06 (n=20 subintrvls)
--- m=4 ---
J=2, calls=10 | DiffScheme=0.51930554, err=2.07e-03 | Simpson=0.51731432, err=7.87e-05 (n=8 subintrvls)
J=4, calls=12 | DiffScheme=0.51724824, err=1.27e-05 | Simpson=0.51726737, err=3.18e-05 (n=10 subintrvls)
J=8, calls=16 | DiffScheme=0.51723561, err=2.63e-08 | Simpson=0.51724376, err=8.17e-06 (n=14 subintrvls)
J=16, calls=24 | DiffScheme=0.51723559, err=3.22e-11 | Simpson=0.51723692, err=1.33e-06 (n=22 subintrvls)
--- m=5 ---
J=2, calls=12 | DiffScheme=0.51896190, err=1.73e-03 | Simpson=0.51726737, err=3.18e-05 (n=10 subintrvls)
J=4, calls=14 | DiffScheme=0.51724341, err=7.82e-06 | Simpson=0.51725080, err=1.52e-05 (n=12 subintrvls)
J=8, calls=18 | DiffScheme=0.51723559, err=6.69e-09 | Simpson=0.51724036, err=4.78e-06 (n=16 subintrvls)
J=16, calls=26 | DiffScheme=0.51723559, err=2.47e-12 | Simpson=0.51723653, err=9.38e-07 (n=24 subintrvls)
--- m=6 ---
J=2, calls=14 | DiffScheme=0.51877832, err=1.54e-03 | Simpson=0.51725080, err=1.52e-05 (n=12 subintrvls)
J=4, calls=16 | DiffScheme=0.51724102, err=5.43e-06 | Simpson=0.51724376, err=8.17e-06 (n=14 subintrvls)
J=8, calls=20 | DiffScheme=0.51723559, err=2.15e-09 | Simpson=0.51723856, err=2.98e-06 (n=18 subintrvls)
J=16, calls=28 | DiffScheme=0.51723559, err=2.63e-13 | Simpson=0.51723627, err=6.81e-07 (n=26 subintrvls)
--- m=7 ---
J=2, calls=16 | DiffScheme=0.51867099, err=1.44e-03 | Simpson=0.51724376, err=8.17e-06 (n=14 subintrvls)
J=4, calls=18 | DiffScheme=0.51723962, err=4.03e-06 | Simpson=0.51724036, err=4.78e-06 (n=16 subintrvls)
J=8, calls=22 | DiffScheme=0.51723559, err=8.48e-10 | Simpson=0.51723754, err=1.95e-06 (n=20 subintrvls)
J=16, calls=30 | DiffScheme=0.51723559, err=3.97e-14 | Simpson=0.51723609, err=5.06e-07 (n=28 subintrvls)
=========== Пример 2 ===========
Интеграл: I2 = ∫[-1..1] e^(-x^2) dx = √π * erf(1)
--- m=3 ---
J=2, calls= 8 | DiffScheme=1.49190419, err=1.74e-03 | Simpson=1.49383992, err=1.92e-04 (n=6 subintrvls)
J=4, calls=10 | DiffScheme=1.49360711, err=4.12e-05 | Simpson=1.49371076, err=6.25e-05 (n=8 subintrvls)
J=8, calls=14 | DiffScheme=1.49364800, err=2.64e-07 | Simpson=1.49366078, err=1.25e-05 (n=12 subintrvls)
J=16, calls=22 | DiffScheme=1.49364826, err=1.17e-09 | Simpson=1.49364990, err=1.63e-06 (n=20 subintrvls)
--- m=4 ---
J=2, calls=10 | DiffScheme=1.49094566, err=2.70e-03 | Simpson=1.49371076, err=6.25e-05 (n=8 subintrvls)
J=4, calls=12 | DiffScheme=1.49361774, err=3.05e-05 | Simpson=1.49367411, err=2.58e-05 (n=10 subintrvls)
J=8, calls=16 | DiffScheme=1.49364820, err=6.19e-08 | Simpson=1.49365504, err=6.77e-06 (n=14 subintrvls)
J=16, calls=24 | DiffScheme=1.49364827, err=7.29e-11 | Simpson=1.49364938, err=1.11e-06 (n=22 subintrvls)
--- m=5 ---
J=2, calls=12 | DiffScheme=1.49045179, err=3.20e-03 | Simpson=1.49367411, err=2.58e-05 (n=10 subintrvls)
J=4, calls=14 | DiffScheme=1.49362749, err=2.08e-05 | Simpson=1.49366078, err=1.25e-05 (n=12 subintrvls)
J=8, calls=18 | DiffScheme=1.49364825, err=1.42e-08 | Simpson=1.49365224, err=3.98e-06 (n=16 subintrvls)
J=16, calls=26 | DiffScheme=1.49364827, err=4.50e-12 | Simpson=1.49364905, err=7.87e-07 (n=24 subintrvls)
--- m=6 ---
J=2, calls=14 | DiffScheme=1.49016512, err=3.48e-03 | Simpson=1.49366078, err=1.25e-05 (n=12 subintrvls)
J=4, calls=16 | DiffScheme=1.49363399, err=1.43e-05 | Simpson=1.49365504, err=6.77e-06 (n=14 subintrvls)
J=8, calls=20 | DiffScheme=1.49364826, err=3.40e-09 | Simpson=1.49365075, err=2.48e-06 (n=18 subintrvls)
J=16, calls=28 | DiffScheme=1.49364827, err=2.92e-13 | Simpson=1.49364884, err=5.72e-07 (n=26 subintrvls)
--- m=7 ---
J=2, calls=16 | DiffScheme=1.48998437, err=3.66e-03 | Simpson=1.49365504, err=6.77e-06 (n=14 subintrvls)
J=4, calls=18 | DiffScheme=1.49363818, err=1.01e-05 | Simpson=1.49365224, err=3.98e-06 (n=16 subintrvls)
J=8, calls=22 | DiffScheme=1.49364826, err=8.50e-10 | Simpson=1.49364990, err=1.63e-06 (n=20 subintrvls)
J=16, calls=30 | DiffScheme=1.49364827, err=1.95e-14 | Simpson=1.49364869, err=4.25e-07 (n=28 subintrvls)
=========== Пример 3 ===========
Интеграл: I3 = ∫[-1..1] sin^2(100x) dx = 1 - sin(200)/200
--- m=3 ---
J=2, calls= 8 | DiffScheme=0.17756175, err=8.27e-01 | Simpson=0.81310534, err=1.91e-01 (n=6 subintrvls)
J=4, calls=10 | DiffScheme=0.17763152, err=8.27e-01 | Simpson=0.17760894, err=8.27e-01 (n=8 subintrvls)
J=8, calls=14 | DiffScheme=0.17763182, err=8.27e-01 | Simpson=1.08564473, err=8.13e-02 (n=12 subintrvls)
J=16, calls=22 | DiffScheme=0.17763182, err=8.27e-01 | Simpson=1.07678341, err=7.24e-02 (n=20 subintrvls)
--- m=4 ---
J=2, calls=10 | DiffScheme=0.17761897, err=8.27e-01 | Simpson=0.17760894, err=8.27e-01 (n=8 subintrvls)
J=4, calls=12 | DiffScheme=0.17763181, err=8.27e-01 | Simpson=1.10415952, err=9.98e-02 (n=10 subintrvls)
J=8, calls=16 | DiffScheme=0.17763182, err=8.27e-01 | Simpson=0.85172200, err=1.53e-01 (n=14 subintrvls)
J=16, calls=24 | DiffScheme=0.17763182, err=8.27e-01 | Simpson=0.88096640, err=1.23e-01 (n=22 subintrvls)
--- m=5 ---
J=2, calls=12 | DiffScheme=0.17762932, err=8.27e-01 | Simpson=1.10415952, err=9.98e-02 (n=10 subintrvls)
J=4, calls=14 | DiffScheme=0.17763182, err=8.27e-01 | Simpson=1.08564473, err=8.13e-02 (n=12 subintrvls)
J=8, calls=18 | DiffScheme=0.17763182, err=8.27e-01 | Simpson=0.17763040, err=8.27e-01 (n=16 subintrvls)
J=16, calls=26 | DiffScheme=0.17763182, err=8.27e-01 | Simpson=0.95775150, err=4.66e-02 (n=24 subintrvls)
--- m=6 ---
J=2, calls=14 | DiffScheme=0.17763131, err=8.27e-01 | Simpson=1.08564473, err=8.13e-02 (n=12 subintrvls)
J=4, calls=16 | DiffScheme=0.17763182, err=8.27e-01 | Simpson=0.85172200, err=1.53e-01 (n=14 subintrvls)
J=8, calls=20 | DiffScheme=0.17763182, err=8.27e-01 | Simpson=0.85501868, err=1.49e-01 (n=18 subintrvls)
J=16, calls=28 | DiffScheme=0.17763182, err=8.27e-01 | Simpson=1.07412498, err=6.98e-02 (n=26 subintrvls)
--- m=7 ---
J=2, calls=16 | DiffScheme=0.17763172, err=8.27e-01 | Simpson=0.85172200, err=1.53e-01 (n=14 subintrvls)
J=4, calls=18 | DiffScheme=0.17763182, err=8.27e-01 | Simpson=0.17763040, err=8.27e-01 (n=16 subintrvls)
J=8, calls=22 | DiffScheme=0.17763182, err=8.27e-01 | Simpson=1.07678341, err=7.24e-02 (n=20 subintrvls)
J=16, calls=30 | DiffScheme=0.17763182, err=8.27e-01 | Simpson=1.03893703, err=3.46e-02 (n=28 subintrvls)