Как реализовать по-настоящему эффективные методы численного интегрирования с нуля? ⇐ Python

[Anonymous]

Это будет длинный пост.
Этот вопрос не касается домашнего задания и/или работы, я безработный, бросил школу в 2017 году, мне сейчас 26 лет, и я изучил программирование полностью самостоятельно.
Прежде чем вы скажете «просто используйте числовую библиотеку», меня не волнует практичность, я пишу программы просто для развлечения. А методы численного интегрирования произвольной точности я уже реализовал на Python с нуля, они полностью рабочие, я их все написал полностью сам, большинство из них используют чистую целочисленную арифметику, я нашел формулы в Википедии, упростил математику и написал весь код. Но они не сходятся достаточно быстро. Мне нужны более точные цифры за меньшее количество итераций.
Как же мне удалось повсюду использовать целочисленную арифметику? Это просто, я просто использую дроби. Основная идея заключается в следующем: сложение двух дробей таково: a/b + c/d = (a * d + b * c)/(b * d), затем мы можем просто использовать пару целых чисел для разделения числителя и знаменателя, мы определяем (a, b) + (c, d) = (a * d + b * c, b * d), таким образом, числа всегда имеют бесконечную точность, в отличие от просто двойных чисел IEEE-754. log10(2)*53=15.954589770191003 десятичных цифр.
Сначала я использовал формулы Ньютона-Котеса из Википедии
Вот формулы:


Самая простая формула, которую я реализовал, - это правило Буля 4-го порядка, я реализовал порядки 4, 5 и 6.
Я вывел точные формулы для расчета арктангенса, используя формулы Ньютона-Котеса для порядков 4 и 5. Сначала производная арктангенса равна 1 / (x^2 + 1), выразите x как a/b, то мы имеем 1 / (a^2/b^2 + 1) = 1 / ((a^2 + b^2) / b^2) = b^2 / (a^2 + b^2). Подставьте производную в формулу Ньютона-Котеса, выразите числа с помощью дробей, сократите дроби, сложите дроби, чтобы получить многочлены, а затем используйте метод Хорнера для вычисления многочленов, обеспечивая таким образом наименьшее количество умножений, и тогда мы получим следующее:

Код: Выделить всё

def arctan_Boole(x: int, y: int) -> tuple[int, int]:
y  tuple[int, int]:
return y * y, x * x + y * y

def add_fractions(num1: int, den1: int, num2: int, den2: int) -> tuple[int, int]:
return (num1 * den2 + num2 * den1, den1 * den2)

def Boole_integrate(
low_num: int, low_den: int, high_num: int, high_den: int, integrand: callable
) -> tuple[int, int]:
prod = low_num * high_den
num1, den1 = high_num * low_den - prod, low_den * high_den  tuple[int, int]:
return (num // (common := math.gcd(num, den)), den // common)

def integrate_by_parts(
low_num: int,
low_den: int,
high_num: int,
high_den: int,
level: int,
integrand: callable,
integrator: callable,
reduce: bool = False,
) -> tuple[int, int]:
start = low_num * high_den
step_num, den = high_num * low_den - low_num * high_den, low_den * high_den
points = [(low_num, low_den)]
for i in range(1, 1  tail), den  tuple[int, int]:
return integrate_by_parts(0, 1, x, y, level, arctan_derivative, integrator, reduce)

def arctan_Boole_parts(
x: int,
y: int,
level: int = 2,
reduce: bool = False,
) -> tuple[int, int]:
return arctan_parts(x, y, level, Boole_integrate, reduce)

def arctan_Newton_5_parts(
x: int,
y: int,
level: int = 2,
reduce: bool = False,
) -> tuple[int, int]:
return arctan_parts(x, y, level, Newton_Cotes_5_integrate, reduce)

def arctan_Newton_6_parts(
x: int,
y: int,
level: int = 2,
reduce: bool = False,
) -> tuple[int, int]:
return arctan_parts(x, y, level, Newton_Cotes_6_integrate, reduce)

Все работают. Попробуйте их, не верьте мне на слово.
Я также реализовал Гауссову квадратуру, используя веса из Википедии, я реализовал порядки 3, 4 и 5, веса указаны ниже:


Здесь есть квадратные корни. Теперь, если числа имеют вид a + b * sqrt(N), где N фиксировано, то сложение, вычитание, умножение и деление можно легко определить: (a + b * sqrt(N)) + (c + d * sqrt(N)) = (a + c + (b + d) * sqrt(N)) = (a + c, b + d), (a + b * sqrt(N)) - (c + d * sqrt(N)) = (a - c + (b -) d) * sqrt(N)) = (a - c, b - d), (a + b * sqrt(N)) * (c + d * sqrt(N)) = (a * c + N * b * d + (a * d + b * c) * sqrt(N)) = (a * c + N * b * d, a * d + b * c), (a + b * sqrt(N)) / (c + d * sqrt(N)) = ((a + b *) sqrt(N)) * (c - d * sqrt(N))) / ((c + d * sqrt(N)) * (c - d * sqrt(N))) = (a * c - N * b * d + (b * c - a * d) * sqrt(N)) / (c^2 - N * d^2) = ((a * c - N * b * d, b * c - a * d), (c^2 - N * d^2, 0)).
Но поскольку радикалы неоднородны, описанные выше трюки использовать нельзя, поэтому необходимо получать приближения. Я использовал gmpy2 для вычислений с плавающей запятой с произвольной точностью:

Код: Выделить всё

import gmpy2
import math
from gmpy2 import get_context, mpfr, sqrt

def sqrt_frac(a: int, b: int) -> mpfr:
return sqrt(mpfr(a) / b)

Number = int | float | mpfr

class Gaussian_Integrator:
def __init__(self) -> None:
get_context().precision = 53
self.calculate_constants()

def set_precision(self, decimal_places: int) -> None:
assert isinstance(decimal_places, int) and decimal_places
get_context().precision = int(math.ceil(decimal_places / math.log10(2)))
self.calculate_constants()

def calculate_constants(self) -> None:
root1 = sqrt_frac(3, 5)
frac1 = mpfr(5) / 9
root2 = sqrt_frac(6, 5) * 2 / 7
frac2 = mpfr(3) / 7
root3 = sqrt(frac2 - root2)
root4 = sqrt(frac2 + root2)
root5 = sqrt(30)
frac3 = (18 + root5) / 36
frac4 = (18 - root5) / 36
root6 = sqrt_frac(10, 7) * 2
root7 = sqrt(5 - root6) / 3
root8 = sqrt(5 + root6) / 3
root9 = sqrt(70) * 13
frac5 = (322 + root9) / 900
frac6 = (322 - root9) / 900
self.data = [
[(-root1, frac1), (0, mpfr(8) / 9), (root1, frac1)],
[(-root4, frac4), (-root3, frac3), (root3, frac3), (root4, frac4)],
[
(-root8, frac6),
(-root7, frac5),
(0, mpfr(128) / 225),
(root7, frac5),
(root8, frac6),
],
]

def integrate(
self, low: Number, high: Number, integrand: callable, level: int
) -> mpfr:
assert (
isinstance(low, Number)
and isinstance(high, Number)
and low < high
and isinstance(level, int)
and 0  None:
self.integrator = Gaussian_Integrator()

def set_precision(self, decimal_places: int) -> None:
self.integrator.set_precision(decimal_places)

def calculate(self, x: Number, level: int = 2) -> mpfr:
return self.integrator.integrate(0, x, Arctan_Gaussian.integrand, level)

def calculate_parts(self, x: Number, bisect_level: int, level: int = 2) -> mpfr:
return self.integrator.integrate_by_parts(
0, x, Arctan_Gaussian.integrand, bisect_level, level
)

ATAN_GAUSS = Arctan_Gaussian()
ATAN_GAUSS.set_precision(100)

И снова мой код работает отлично.
Что дальше? Я реализовал метод Ромберга:


Он выглядит довольно сложным, но на самом деле он довольно прост:


Итак, код:

Код: Выделить всё

def Romberg_integration(
low_num: int,
low_den: int,
high_num: int,
high_den: int,
level: int,
integrand: callable,
) -> tuple[int, int]:
step_num = high_num * low_den - low_num * high_den
step_den = low_den * high_den 

Подробнее здесь: [url]https://stackoverflow.com/questions/79885495/how-to-implement-truly-efficient-numerical-integration-methods-from-scratch[/url]