SciPy में अनंत समाकलन: quad की qagie मैपिंग को समझें

इस लेख में SciPy के quad द्वारा अनंत समाकलन की QUADPACK qagie मैपिंग समझें: X=A+(1−T)/T रूपांतरण, [0,1] पर इंटीग्रेशन, स्थिरता व points उपयोग के सुझाव.

अनंत अंतरालों पर फलनों का समाकलन SciPy में quad कॉल करने पर देखने में सरल लगता है, लेकिन अंदर एक सटीक तंत्र काम करता है। QUADPACK की qagie रूटीन में उपयोग होने वाले निर्देशांक-परिवर्तन को समझने से संख्यात्मक स्थिरता, अभिसरण, और सीमित डोमेन पर points जैसी सुविधाओं के साथ काम करने की समझ बेहतर होती है।

समस्या की रूपरेखा

यह न्यूनतम उदाहरण SciPy के साथ 0 से अनंत तक मानक नॉर्मल PDF का समाकलन करता है। उद्देश्य यह देखना है कि quad अनंत सीमाओं को कैसे संभालता है।

import scipy
import numpy as np
def pdf_normal(z):
    return np.exp(-(z**2) / 2) / np.sqrt(2 * np.pi)
val, err, dbg = scipy.integrate.quad(pdf_normal, 0.0, np.inf, full_output=True)
print(val)

SciPy के दस्तावेज़ बताते हैं कि अनंत सीमाएँ QUADPACK की qagie के माध्यम से निर्देशांक-परिवर्तन से संभाली जाती हैं। उनके शब्दों में:

qagie

अनंत अंतरालों पर समाकलन को संभालती है। अनंत दायरे को एक सीमित अंतराल पर प्रतिचित्रित किया जाता है और इसके बाद QAGS की वही रणनीति लागू की जाती है।

दिलचस्प बात यह है कि [0, 1] के सीमित अंतराल पर सटीक मैपिंग कैसी होती है, और उसका इंटेग्रैंड के आकार तथा संख्यात्मक बर्ताव पर क्या असर पड़ता है।

अंदरूनी तौर पर वास्तव में क्या होता है

qagie द्वारा उपयोग किया गया मैपिंग स्पष्ट रूप से परिभाषित है। अर्ध-अनंत समाकलन में, जहाँ निचली सीमा A है, (0, 1] के T को [A, ∞) के X पर X = A + (1 − T)/T से मैप किया जाता है। यह रूपांतरण मूल समाकलन को इकाई अंतराल पर बदल देता है, पर संशोधित इंटेग्रैंड के साथ। दो-पक्षीय अनंत सीमा के लिए भी यही विचार सममित रूप से लागू होता है: हर T के लिए फलन का मान X और −X दोनों पर लिया जाता है।

विवरण यहाँ दर्ज हैं: Piessens, Robert; de Doncker-Kapenga, Elise; Überhuber, Christoph W.; Kahaner, David (1983). QUADPACK: A subroutine package for automatic integration. Springer-Verlag. पृष्ठ 64–65.

एक-पक्षीय अनंत सीमा: [A, ∞)

प्रतिस्थापन X = A + (1 − T)/T लें, जहाँ T ∈ (0, 1]। जब T → 1, तब X → A। जब T → 0+, तब X → ∞। अवकलन करने पर dX/dT = −T^(−2) मिलता है। सीमाएँ उलटने पर ऋण चिह्न हट जाता है, और T = 0 से 1 तक का समाकलन मिलता है, जिसका इंटेग्रैंड F(A + (1 − T)/T) · T^(−2) होता है। qagie अर्ध-अनंत समाकलन को [0, 1] पर इसी मूल तरीके से मैप करती है।

दो-पक्षीय अनंत सीमाएँ: (−∞, ∞)

पूरी वास्तविक रेखा पर समाकलन के लिए, यही मैपिंग विचार सममिति के साथ लागू होता है। प्रत्येक T ∈ (0, 1] के लिए X = (1 − T)/T लिया जाता है, और इंटेग्रैंड को F(X) व F(−X) के रूप में जोड़ा जाता है: F(X) + F(−X), तथा भार T^(−2) रहता है।

व्यावहारिक समाधान: मैनुअल पुनर्पैरामीट्रीकरण

आप यही प्रतिस्थापन स्वयं भी लगा सकते हैं, [0, 1] पर समाकलन कर सकते हैं, और देख सकते हैं कि रूपांतरित इंटेग्रैंड कैसा दिखता है। यह डायग्नोस्टिक्स के लिए उपयोगी है और points जैसी उन सुविधाओं के उपयोग में भी मदद करता है जो सीधे अनंत सीमाओं के साथ संगत नहीं होतीं।

import scipy
import numpy as np
import matplotlib.pyplot as plt
def pdf_normal(z):
    return np.exp(-(z**2) / 2) / np.sqrt(2 * np.pi)
def reparam_semi_infinite(fn, a):
    """Map T in (0, 1] to X in [a, +inf) via X = a + (1 - T)/T."""
    def wrapped(t):
        return fn(a + (1 - t) / t) * (t ** -2)
    return wrapped
def reparam_full_line(fn):
    """Map T in (0, 1] to X in (-inf, +inf) using symmetry."""
    def wrapped(t):
        x = (1 - t) / t
        return (fn(x) + fn(-x)) * (t ** -2)
    return wrapped
# अर्ध-अनंत मामला: T ∈ [eps, 1] के माध्यम से 0 से +inf तक सामान्य PDF का समाकलन करें
pdf_mapped = reparam_semi_infinite(pdf_normal, a=0.0)
# रूपांतरित इंटेग्रैंड T=0 पर सिंगुलर है (यह X = +inf के अनुरूप है),
# इसलिए 0 से थोड़ा हटकर शुरू करें।
t0 = 1e-8
res, err, meta = scipy.integrate.quad(pdf_mapped, t0, 1.0, full_output=True)
print(res)
# मूल और रूपांतरित इंटेग्रैंड का दृश्यांकन
xs1 = np.linspace(0.0, 5.0, 100)
plt.title("Normal PDF in X domain")
plt.plot(xs1, pdf_normal(xs1))
plt.show()
xs2 = np.linspace(t0, 1.0, 100)
plt.title("Transformed integrand in T domain")
plt.plot(xs2, pdf_mapped(xs2))
plt.show()

T के संदर्भ में समाकलन T = 0 पर परिभाषित नहीं है, क्योंकि यह प्रतिस्थापन के तहत X = ∞ को दर्शाता है; इसलिए T की निचली सीमा शून्य से थोड़ी ऊपर लेनी पड़ती है। इस सेटअप में आप [t0, 1] पर समाकलन करेंगे और मूल अर्ध-अनंत समाकलन जैसा ही मान प्राप्त होगा।

यदि आपको इंटेग्रैंड के लिए ब्रेकपॉइंट देने हों, तो अब domain [0, 1] होने की वजह से points आर्ग्युमेंट का उपयोग किया जा सकता है। मूल X-डोमेन में परिभाषित ब्रेकपॉइंट्स को सूत्र t = 1 / (1 − a + x) के ज़रिए T में ले जाया जा सकता है, जहाँ a मूल समाकलन की निचली सीमा है।

यह क्यों मायने रखता है

यह मैपिंग तीन व्यावहारिक पहलुओं को स्पष्ट करती है। पहला, T = 0 पर जो भी singularity दिखती है, वह रूपांतरण की उपज है और मूल चर में अनंत पर व्यवहार को दर्शाती है; इससे इंटीग्रेटर की कठिनाई समझ आती है या T = 0 को बिल्कुल टालने की ज़रूरत का कारण मिलता है। दूसरा, रूपांतरित इंटेग्रैंड का निरीक्षण करने से पता चल सकता है कि endpoints के पास आकार कितनी तेज़ी से बदल रहा है और इंटीग्रेटर को कहाँ संघर्ष करना पड़ सकता है। तीसरा, जब समाकलन [0, 1] पर आ जाता है, तो points जैसी सुविधाएँ उपलब्ध हो जाती हैं; दिए गए सूत्र से ब्रेकपॉइंट्स को स्थानांतरित करने के बाद, तीखे बदलावों के आसपास विश्वसनीयता बेहतर की जा सकती है।

निष्कर्ष

quad अनंत दायरों को [0, 1] पर पुनर्पैरामीट्राइज़ करके संभालता है। [A, ∞) के लिए यह X = A + (1 − T)/T और भार T^(−2) का उपयोग करता है; (−∞, ∞) के लिए वही भार रखते हुए दोनों टेल्स को सममित रूप से F(X) + F(−X) के रूप में जोड़कर मान लेता है। आप इस रूपांतरण को मैन्युअल रूप से दोहरा कर इंटेग्रैंड को देख सकते हैं, t = 1/(1 − a + x) से ब्रेकपॉइंट्स रख सकते हैं, और समझ सकते हैं कि संख्यात्मक चुनौतियाँ कहाँ उभर सकती हैं। यदि T के सिरों के पास दिक्कत दिखे, तो याद रखें कि वह मूल चर में अनंत पर के व्यवहार का प्रतिबिंब है—तदनुसार अपनी सीमाएँ या डायग्नोस्टिक्स समायोजित करें।