Метод квадратурного численного интегрирования Tanh-sinh сходится к неправильному значению

Для чего это стоит, Scipy имеет функции численного интегрирования

Например,

from scipy import integrate
check = integrate.quad(lambda x: 1 / math.sqrt(1 - x ** 2), -1, 1)
print 'Scipy quad integral = ', check

дает результат

Квадратный интеграл Scipy = (3.141592653589591, 6.200897573194197e-10)

где второе число в кортеже — абсолютная ошибка.

Тем не менее, я смог заставить вашу программу работать с некоторой настройкой (хотя это всего лишь начальная попытка):

1) Установите размер шага h равным 0,0002 (примерно 1/2^12), как предлагается эта статья

Но обратите внимание: в документе на самом деле предлагается итеративно изменять размер шага — при фиксированном размере шага вы достигнете точки, в которой sinh или ch становятся слишком большими для достаточно больших значений kh. Вероятно, было бы лучше попытаться реализовать подход, основанный на этом документе.

Но придерживаясь поставленного вопроса,

2) Убедитесь, что вы установили достаточное количество итераций для действительной сходимости интегрирования, т. е. достаточное количество итераций, чтобы math.fabs(w_k * func(x_k)) ‹ 1.0e-9

Благодаря этим настройкам мне удалось добиться сходимости интегрирования к правильному значению числа пи до 4 значащих цифр, используя > 30000 итераций.

Например, при 31111 итерациях вычисленное значение числа пи составило 3,14159256208.

Пример кода с изменениями (обратите внимание, я заменил sum на thesum, sum — это имя встроенной функции Python):

import math

def func(x):
    # Function to be integrated, with singular points set = 0
    if x == 1 or x == -1 :
        return 0
    else:
        return 1 / math.sqrt(1 - x ** 2)

# Input number of evaluations
N = input("Please enter number of evaluations \n")
if N % 2 == 0:
    print "The number of evaluations must be odd"
else:
    print "N =", N  

# Set step size
#h = 2.0 / (N - 1)
h=0.0002 #(1/2^12)
print "h =", h

# k ranges from -(N-1)/2 to +(N-1)/2
k = -1 * ((N - 1) / 2.0)
k_max  = ((N - 1) / 2.0)
thesum = 0

# Loop across integration interval
actual_iter =0
while k < k_max + 1:

    # Compute abscissa
    x_k = math.tanh(math.pi * 0.5 * math.sinh(k * h))

    # Compute weight
    numerator = 0.5 * h * math.pi * math.cosh(k * h)
    dcosh  = math.cosh(0.5 * math.pi * math.sinh(k * h))
    denominator = dcosh*dcosh
    #denominator = math.pow(math.cosh(0.5 * math.pi * math.sinh(k * h)),2)
    w_k =  numerator / denominator

    thesum += w_k * func(x_k)
    myepsilon = math.fabs(w_k * func(x_k))
    if actual_iter%2000 ==0 and actual_iter > k_max/2:
        print "Iteration = %d , myepsilon = %g"%(actual_iter,myepsilon)


    k += 1
    actual_iter += 1

print 'Actual iterations = ',actual_iter
print "Integral =", thesum

Используя библиотеку мультиточности mpmath:

from mpmath import *

mp.dps = 100

h = mpf(2**-12);

def weights(k):
    num = mpf(0.5)*h*pi*cosh(k*h)
    den = cosh(mpf(0.5)*pi*sinh(k*h))**2
    return (num/den)

def abscissas(k):
    return tanh(mpf(0.5)*pi*sinh(k*h))

def f(x):
    return 1/sqrt(1 - mpf(x)**2)

N = 20000

result = 0
for k in range(-N, N+1):
    result = result + weights(k)*f(abscissas(k))

print result - pi

дает за ошибку

-3.751800610920472803216259350430460844457732874052618682441090144344372471319795201134275503228835472e-45

Я думаю, что часть проблемы может быть связана с диапазоном и размером шага. Я изменил код, чтобы вы могли указать диапазон и размер шага отдельно, и переписал некоторые математические операции. Кажется, он дает правильные ответы. Попробуйте, например, 5 и 0,1 в качестве входных данных.

Особая проблема заключается в вычислении math.cosh(0.5 * math.pi * math.sinh(k * h)) по мере того, как k * h становится большим math.sinh(k * h) растет экспоненциально, и вычисление math.cosh может быть трудным. импортировать математику

def func(x):
#    return 1   # very simple test function
    # Function to be integrated, with singular points set = 0
    if x == 1 or x == -1 :
        return 0
    else:
        return 1 / math.sqrt(1 - x ** 2)

# Input number of evaluations
N = input("Please enter max value for range \n")
    print "N =", N
h = input("Please the step size \n")
print "h =", h

k = -N
k_max = N
sum = 0
count = 0
print "k ", k , " " , k_max

# Loop across integration interval
while k < k_max + 1:

    # Compute abscissa
    v = k
    u = math.pi * 0.5 * math.sinh(v)
    x_k = math.tanh(u)
    #print u
    # Compute weight 
    numerator = 0.5 * math.pi * math.cosh(v)
    csh = math.cosh(u)
    denominator = csh*csh
    w_k =  numerator / denominator
    print k, v, x_k , w_k
    sum += w_k * func(x_k)
    count += 1
    k += h      # note changed
res = sum * h
print "Integral =", sum * h

Вы должны понимать, что +1 и -1 являются особыми точками вашего интеграла, f(x)-->+infinity как x-->+1,-1. Таким образом, вы можете использовать свою любимую квадратурную формулу вдали от граничных точек, но вам нужно разработать специальную квадратуру, основанную на локальном расширении f(x) в окрестности их.

Эскиз подхода:

1. Выберите несколько epsilon<<1.

2. Разобьем интеграл I на гладкую и сингулярную части:

   • I_smooth is the integral inside [-1+epsilon, 1-epsilon]
   • I_singular — это интегралы от [-1, -1+epsilon] и [1-epsilon, 1].

3. Примените стандартное квадратурное правило внутри интервала [-1+epsilon, 1-epsilon], чтобы получить I_smooth.

4. Выполните локальное расширение вокруг особых точек (например, x=1):

   f(x) = 1/sqrt(1-x) * (a0 + a1*(1-x) + a2*(1-x)^2 + ...)
   
        = f0(x-1) + f1(x-1) + f2(x-1) + ..
   

   что является просто разложением Тейлора о x=1 от f(x)*sqrt(1-x), предварительно умноженного на 1/sqrt(1-x). (К сожалению, вам придется заняться математикой и вычислить разложение Тейлора, если только у вас нет Mathematica или вы не найдете его где-нибудь в виде таблицы.)

5. Каждое отдельное слагаемое fn(x-1) = an*(1-x)^n/sqrt(1-x) можно точно проинтегрировать (это просто степенная функция). Пусть Fn будет интегралом от fn от 1-epsilon до 1. Приблизительно I_singular = F0 + F1 + F2 + ... до нужного вам порядка.

6. Окончательно:

   I = I_smooth + I_singular  
   

Примечание: чтобы повысить точность, вы не должны делать epsilon слишком маленьким, потому что разрушение интеграла делает задачу численно плохо обусловленной, а скорее увеличиваете порядок разложения Тейлора.

Когда дело доходит до квадратуры тангенса и синха, возникает множество ловушек, одна из которых заключается в том, что подынтегральная функция должна оцениваться очень близко к границам интервала, на расстояниях, меньших машинной точности, например, 1.0 - 1.0e-20 в оригинальный пример. Когда эта точка оценивается, она округляется до 1.0, в которой f имеет сингулярность, и может случиться что угодно. Вот почему вам сначала нужно преобразовать функцию так, чтобы сингулярности находились в 0.

В случае 1 / sqrt(1 - x**2) это 1 / numpy.sqrt(-x**2 + 2*x) как для левой, так и для правой особенности. С помощью tanh_sinh (мой проект) можно получить

import numpy
import tanh_sinh

# def f(x):
#    return 1 / numpy.sqrt(1 - x ** 2)

val, error_estimate = tanh_sinh.integrate_lr(
      [lambda x: 1 / numpy.sqrt(-x**2 + 2*x)],  # = 1 / sqrt(1 - (x-1)**2)
      [lambda x: 1 / numpy.sqrt(-x**2 + 2*x)],  # = 1 / sqrt(1 - (-(x-1))**2)
      2,  # length of the interval
      1.0e-10
      )
print(val, val - numpy.pi)

3.1415926533203944 -2.693987255497632e-10