Gauss-Legendre su intervalli -x -> infinito: algoritmo adattivo per trasformare efficientemente pesi e nodi

Question

Ok so che questo è stato chiesto prima con un esempio limitato per il ridimensionamento [-1, 1] intervalli [a, b]Different intervals for Gauss-Legendre quadrature in numpy MA nessuno ha pubblicato come generalizzare questo per [-a, Infinity] (come si fa di seguito, ma non (ancora) veloce). Anche questo mostra come chiamare una funzione complessa (in ogni caso il prezzo delle opzioni quantitative) con diverse implementazioni. Esiste il codice di riferimento quad, seguito da leggauss, con collegamenti a esempi di codice su come implementare un algoritmo adattivo. Ho lavorato sulla maggior parte delle difficoltà collegate adaptive algorithm - attualmente stampa la somma dell'integrale diviso per mostrare che funziona correttamente. Qui troverai le funzioni per convertire un intervallo da [-1, 1] a [0, 1] a [a, Infinity] (grazie @AlexisClarembeau). Per utilizzare l'algoritmo adattivo, ho dovuto creare un'altra funzione per convertire da [-1, 1] a [a, b] che viene reimmesso nella funzione [a, Infinity].

import numpy as np 
from scipy.stats import norm, lognorm 
from scipy.integrate import quad 

a = 0 
degrees = 50 

flag=-1.0000 
F = 1.2075 
K = 0.1251 
vol = 0.43 
T2 = 0.0411 
T1 = 0.0047 

def integrand(x, flag, F, K, vol, T2, T1): 
    d1 = (np.log(x/(x+K)) + 0.5 * (vol**2) * (T2-T1))/(vol * np.sqrt(T2 - T1)) 
    d2 = d1 - vol*np.sqrt(T2 - T1) 
    mu = np.log(F) - 0.5 *vol **2 * T1 
    sigma = vol * np.sqrt(T1) 
    return lognorm.pdf(x, mu, sigma) * (flag * x*norm.cdf(flag * d1) - flag * (x+K)*norm.cdf(flag * d2)) 

def transform_integral_0_1_to_Infinity(x, a): 
    return integrand(a+(x/(1-x)), flag, F, K, vol, T2, T1) *(1/(1-x)**2); 

def transform_integral_negative1_1_to_0_1(x, a): 
    return 0.5 * transform_integral_0_1_to_Infinity((x+1)/2, a) 

def transform_integral_negative1_1_to_a_b(x, w, a, b): 
    return np.sum(w*(0.5 * transform_integral_0_1_to_Infinity(((x+1)/2*(b-a)+a), a))) 

def adaptive_integration(x, w, a=-1, b=1, lastsplit=False, precision=1e-10): 
    #split the integral in half assuming [-1, 1] range 
    midpoint = (a+b)/2 
    interval1 = transform_integral_negative1_1_to_a_b(x, w, a, midpoint) 
    interval2 = transform_integral_negative1_1_to_a_b(x, w, midpoint, b) 
    return interval1+interval2 #just shows this is correct for splitting the interval 

def integrate(x, w, a): 
    return np.sum(w*transform_integral_negative1_1_to_0_1(x, a)) 

x, w = np.polynomial.legendre.leggauss(degrees) 
quadresult = quad(integrand, a, np.Inf, args=(flag, F, K, vol, T2, T1), epsabs=1e-1000)[0] 
GL = integrate(x, w, a) 
print("Adaptive Sum Result:") 
print(adaptive_integration(x, w)) 
print("GL result"); 
print(GL) 
print("QUAD result") 
print(quadresult) 

ancora bisogno di aumentare la velocità e la precisione con minori dimensioni non posso regolare manualmente l'intervallo degrees per -a ottenere convergenza. Per illustrare il motivo per cui questo è un problema, inserisci questi valori: a=-20, F=50, quindi esegui. È possibile aumentare degrees=1000 e vedere che non vi è alcun vantaggio per questo algoritmo di Gauss-Legendre se non viene applicato in modo intelligente. Il mio requisito per la velocità è di arrivare a 0.0004s per loop, mentre l'ultimo algoritmo I Cythonized ha richiesto circa 0,75 secondi, motivo per cui sto cercando di utilizzare un algoritmo di bassa precisione e alta precisione con Gauss-Legendre. Con Cython e multi-filettatura questo requisito da un'implementazione Python completamente ottimizzata è approssimativamente 0.007s per loop (un cappio guidato, di routine non vectorized inefficiente potrebbe essere 0.1s per loop, con degrees=20, cioè %timeit adaptive_integration(x,w).

Una possibile soluzione che ho implementato a metà è qui http://online.sfsu.edu/meredith/Numerical_Analysis/improper_integrals alle pagine 5/6, adaptive integration mentre l'intervallo a-b (in questo caso, ho scritto la funzione transform_integral_negative1_1_to_a_b) dove l'intervallo è diviso in 2 (@0.5), la funzione viene quindi valutata su questi 1/2 intervalli e la somma dei due 0->0.5 + 0.5->1 viene confrontata con i risultati della funzione per l'intero intervallo 0->1. Se la precisione non è entro la tolleranza, l'intervallo è ulteriormente s suddiviso in 0.25 e 0.75, la funzione viene nuovamente valutata per ciascun sottointervallo e confrontata con le somme dell'intervallo 1/2 precedenti @0.5. Se 1 lato è entro la tolleranza (ad esempio abs(0->0.5 - (0->0.25 + 0.25->0.5)) < precision), ma l'altro lato non lo è, la divisione si arresta sul lato entro la tolleranza, ma continua dall'altro lato finché non viene raggiunto precision. A questo punto i risultati per ciascuna porzione dell'intervallo vengono sommati per ottenere l'integrale completo con maggiore precisione.

Ci sono probabilmente modi più veloci e migliori per affrontare questo problema. Non mi interessa finché è veloce e preciso. Ecco la migliore descrizione delle routine di integrazione che ho trovato per riferimento http://orion.math.iastate.edu/keinert/computation_notes/chapter5.pdf Il premio è 100 punti bounty + 15pts per l'accettazione della risposta. Grazie per aver contribuito a rendere questo codice VELOCE e ACCURATO!

EDIT:

Qui è la mia modifica al codice adaptive_integration - se qualcuno può fare questo lavoro velocemente posso accettare una risposta e di aggiudicazione di taglie. Questo codice Mathematica a pagina 7 http://online.sfsu.edu/meredith/Numerical_Analysis/improper_integrals esegue la routine che ho tentato. Ha funzionato su una routine che non converge bene, vedi le variabili sottostanti.In questo momento i miei errori di codice fuori: RecursionError: maximum recursion depth exceeded in comparison su alcuni input, o se lo degrees è impostato troppo alto, o non si avvicina al risultato quad quando funziona, quindi qualcosa è apparentemente sbagliato qui.

def adaptive_integration(x, w, a, b, integralA2B, remainingIterations, firstIteration, precision=1e-9): 
    #split the integral in half assuming [-1, 1] range 
    if remainingIterations == 0: 
     print('Adaptive integration failed on the interval',a,'->',b) 
    if np.isnan(integralA2B): return np.nan 

    midpoint = (a+b)/2 
    interval1 = transform_integral_negative1_1_to_a_b(x, w, a, midpoint) 
    interval2 = transform_integral_negative1_1_to_a_b(x, w, midpoint, b) 
    if np.abs(integralA2B - (interval1 + interval2)) < precision : 
     return(interval1 + interval2)  
    else: 
     return adaptive_integration(x, w, a, midpoint, interval1, (remainingIterations-1), False) + adaptive_integration(x, w, midpoint, b, interval2, (remainingIterations-1), False) 

#This example doesn't converge to Quad 

# non-converging interval inputs 
a = 0 # AND a = -250 
degrees = 10 

flag= 1 
F = 50 
K = 0.1251 
vol = 0.43 
T2 = 0.0411 
T1 = 0.0047 

print(adaptive_integration(x, w, -1, 1, GL, 500, False)) 

L'uscita con degrees=100 (dopo aver calcolato GL con degrees=10000 per una migliore stima iniziale, in caso contrario, l'algoritmo è sempre d'accordo con la sua accuratezza, apparentemente e non richiama il percorso di adattamento che non riesce ogni volta):

GL result: 
60.065205169286379 
Adaptive Sum Result: 
RecursionError: maximum recursion depth exceeded in comparison 
QUAD result: 
68.338 

Answer 1

Penso che il codice fa il lavoro:

 
    import numpy as np 
    import math 

    deg = 10 
    x, w = np.polynomial.legendre.leggauss(deg) 

    def function(x): 
     # the function to integrate 
     return math.exp(-x) 

    def function2(x, a): 
     return function(a+x/(1-x))/((1-x)**2); 

    def anotherOne(x, a): 
     return 0.5 * function2(x/2 + 1/2, a) 

    def integrate(deg, a): 
     sum = 0 
     x, w = np.polynomial.legendre.leggauss(deg) 
     for i in range(deg): 
      print("sum({}) += {} * {} (eval in {})".format(sum, w[i], anotherOne(x[i], a), x[i])) 
      sum += w[i]*anotherOne(x[i], a) 
     return sum; 

    print("result"); 
    print(integrate(10, 1)) 

Esso combina l'equazione di integrare da una a inf e l'equazione di modificare i limiti di un integrale.

Spero che risolve il problema (funziona per exp (-x) almeno) :)

Se si desidera un calcolo in linea, il programma fa la somma di:

Si tratta di una combinazione di:

E:

E:

Answer 2

in "Programmazione numerica: Una guida pratica per gli scienziati e gli ingegneri che usano Python e C/C++" da Tito A. Beu, è possibile trovare i metodi negli esempi di codice integral.py e specfunc.py qui: http://phys.ubbcluj.ro/~tbeu/INP/libraries.html si chiama la funzione 012.che chiama Laguerre dall'altro file .py e restituisce lo (x,w) corretto tra a e infinity. Ecco come impostare questa funzione (nota appena sopra deg=80 è molto lento, sto solo che vi mostra come applicarlo per la modifica delle linee sopra):

x, w = np.array(xGaussLag(a,deg)) 
gauss = sum(w * integrand(x, flag, F, K, vol, T2, T1)) 

diventa piuttosto stretta convergenza su deg=80 (più veloce), ma io appena messo il eps=1e-13 in xGaussLag e spinto il deg=150 con questi risultati, comunque più veloce di quad del 33%:

la soluzione Quadpack: 0,149221620346 con l'errore: 1.49870924498e-12 Gauss-Legendre soluzione: 0,149238273747 Differenza tra Quadpack e Gauss-Legendre : 1.66534003601e-05

In Cython questo è 6 volte più veloce di Python dritto BTW ancora troppo lento, quindi ho intenzione di provare il pacchetto "FastGL" con la risposta di @Alexis per ora, semplicemente postando come penso che questo sarà essere utile per altri utenti SO in futuro.