Thomas–Hopfield मॉडल के लिए SciPy quad से भरोसेमंद समाकलन

Thomas–Hopfield मॉडल में SciPy quad को points और epsabs/epsrel के साथ ट्यून कर 0 से ∞ तक समाकल स्थिर करें: r_max पर निर्भरता घटाएँ, t=0 पर सिद्धांत से मेल

जब किसी भौतिकी मॉडल, जिसमें दो अपूर्ण समाकल हों, को आप कोड में बदलते हैं, तो सूक्ष्म संख्यात्मक बारीकियाँ आपके ग्राफ को बना भी सकती हैं, बिगाड़ भी सकती हैं। यही बात Thomas–Hopfield मॉडल के साथ scipy.integrate.quad पर आती है: आउटपुट ऊपरी समाकलन-सीमा पर अत्यधिक निर्भर लगता है, और प्लॉटिंग के लिए वेक्टराइज़ करना उलझन बढ़ा देता है। अच्छी बात यह है कि यह व्यवहार अनुकूली क्वाड्रेचर के काम करने के तरीके से सीधे निकलता है, और इसे सिद्धांतसम्मत तरीके से ठीक किया जा सकता है।

समस्या को दोहराना

यह मॉडल r = 0 से अनंत तक दो कर्नेल्स का समाकल करता है और उनसे समय-निर्भर सिग्नल बनाता है। सीधे अनंत तक समाकल करना व्यावहारिक नहीं होता, इसलिए एक सीमित कटऑफ r_max लिया जाता है। r_max के प्रति संवेदनशीलता ही वह लक्षण है जिस पर ध्यान देना चाहिए।

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import quad
# पुनर्संयोजन दर (दूरी के साथ घातीय रूप से घटती है)
def rate_fn(r, w0, r0):
    return w0 * np.exp((-2 * r) / r0)
# मॉडल का पहला कर्नेल
def kern_a(r, w0, r0, t):
    return rate_fn(r, w0, r0) * np.exp(-rate_fn(r, w0, r0) * t) * r**2
# मॉडल का दूसरा कर्नेल (घातांक के अंदर प्रयुक्त)
def kern_b(r, w0, r0, t):
    return (np.exp(-rate_fn(r, w0, r0) * t) - 1) * r**2
# दो समाकलों से बना थॉमस–हॉपफील्ड मॉडल
def hopfield_fn(t, w0, r0, dens, amp, bias):
    r_cap = 250e-7  # अनंत के स्थान पर व्यवहारिक कटऑफ
    val1 = quad(kern_a, 0, r_cap, args=(w0, r0, t))[0]
    val2 = quad(kern_b, 0, r_cap, args=(w0, r0, t))[0]
    return amp * ((4 * np.pi * dens * val1) * np.exp(4 * np.pi * dens * val2)) + bias
# t के array को स्वीकार करने वाला वेक्टराइज़्ड रैपर
hopfield_vec = np.vectorize(hopfield_fn)
# उदाहरण स्वरूप मान
w0 = 1e6
r0 = 1.2e-9
dens = 3e-19
amp = 1
bias = 0
# साहित्यिक सूत्र से t=0 का सैद्धांतिक मान
teor_t0 = dens * np.pi * w0 * r0**3
calc_t0 = hopfield_fn(0, w0, r0, dens, amp, bias)
print('Theoretical value for t=0:', teor_t0)
print('Calculated value for t=0:', calc_t0)
# लॉग-लॉग ग्रिड पर प्लॉट
T = np.linspace(0, 100000000, 100000)
Y = hopfield_vec(T, w0, r0, dens, amp, bias)
plt.plot(T, Y, color='r', linestyle='-', linewidth=0.7)
plt.yscale('log')
plt.xscale('log')
plt.show()

r_cap बदलने से वक्र खिंच जाता है, और t=0 पर मेल कभी-कभी बड़े r_cap की तुलना में छोटे r_cap पर बेहतर दिखता है। यह यादृच्छिक नहीं है; जब अनुकूली इंटीग्रेटर को यह नहीं बताया गया हो कि इंटीग्रैंड का कठिन हिस्सा कहाँ है, तो यही अपेक्षित व्यवहार है।

वास्तव में quad क्या करता है और यह क्यों मायने रखता है

एक नए उपयोगकर्ता के तौर पर जो पहली बार scipy.integrate.quad इस्तेमाल कर रहा है, मैं खुद से पूछ रहा हूँ कि quad परिणाम कैसे निकालता है (यह किस क्रम में एरे लेता है, समाकल करता है और मान लौटाता है)।

quad, Gauss–Kronrod आधारित अनुकूली योजना का उपयोग करता है। यह नोड्स के एक सेट पर इंटीग्रैंड का मान लेता है, मौजूदा उप-अंतराल पर निम्न-क्रम और उच्च-क्रम दोनों आकलन बनाता है, उनके अंतर से स्थानीय त्रुटि का अंदाजा लगाता है, और केवल वहीं परिशोधन करता है जहाँ यह अंतर बताता है कि कुछ रोचक हो रहा है। यदि नमूने लिए गए नोड्स शिखर या तीव्र परिवर्तन वाले क्षेत्र के पास नहीं पड़ते, तो एल्गोरिद्म गलत निष्कर्ष निकाल सकता है कि वहाँ फ़ंक्शन सपाट है और आगे बढ़ जाता है।

एक छोटा उदाहरण इसे ठोस बनाता है। 0 केंद्र वाली मानक नार्मल pdf को चौड़े सममित अंतराल पर समाकल करने पर परिणाम 1 ही मिलता है, जैसा अपेक्षित है।

import scipy
import numpy as np
pdf0 = scipy.stats.norm(loc=0, scale=1)
val, err = scipy.integrate.quad(pdf0.pdf, -1000, 1000)
print('result', val, 'error', err)

अब उसी वितरण को 100 केंद्र पर खिसका दें और वही सीमाएँ रखें। सैद्धांतिक रूप से कुछ नहीं बदलना चाहिए। व्यवहार में, यदि अनुकूली सैम्पलर x=100 के पास झाँक ही नहीं पाया, तो वह हर जगह फ़ंक्शन को नगण्य मानकर शून्य लौटा सकता है।

pdf100 = scipy.stats.norm(loc=100, scale=1)
val, err = scipy.integrate.quad(pdf100.pdf, -1000, 1000)
print('result', val, 'error', err)

इसे ठीक करने के लिए आपको quad को बताना पड़ता है कि फ़ंक्शन की संरचना कहाँ है। points आर्ग्यूमेंट निर्दिष्ट स्थानों के आसपास विभाजन को बाध्य करके ठीक यही काम करता है।

val, err = scipy.integrate.quad(pdf100.pdf, -1000, 1000, points=[100])
print('result', val, 'error', err)

अब Thomas–Hopfield मॉडल पर लौटें। इन इंटीग्रैंड्स में छोटे r पर, मोटे तौर पर 1e-9 से 1e-7 के आसपास, पैमाने-संवेदनशील संरचना होती है। यदि इंटीग्रेटर को वहाँ देखने के लिए उकसाया न जाए, तो परिणाम मनमाने सीमाबदल पर बहुत निर्भर हो जाते हैं। यहाँ एक दूसरी बात भी महत्वपूर्ण है: quad का डिफ़ॉल्ट निरपेक्ष सहनशीलता (epsabs) मान क्रम-एक परिमाणों के लिए बनाया गया है। अत्यंत छोटे मानों के साथ यह epsabs बहुत बड़ा पड़ जाता है। epsabs को शून्य कर देना और एक सापेक्ष सहनशीलता देना इंटीग्रेटर को सही लक्ष्य-सटीकता देता है।

सीमित सीमा वाले मामले के लिए भरोसेमंद समाधान

यह संस्करण quad को साफ-साफ बताता है कि “एक्शन” कहाँ है, और सहनशीलताओं को कड़ा करता है। व्यवहार में एकमात्र बदलाव बेहतर सटीकता और स्थिरता है।

def hopfield_fn(t, w0, r0, dens, amp, bias):
    r_cap = 1.0
    anchors = np.array([1e-9, 1e-8, 1e-7])
    val1 = quad(kern_a, 0, r_cap, args=(w0, r0, t), points=anchors, epsabs=0, epsrel=1.49e-8)[0]
    val2 = quad(kern_b, 0, r_cap, args=(w0, r0, t), points=anchors, epsabs=0, epsrel=1.49e-8)[0]
    return amp * ((4 * np.pi * dens * val1) * np.exp(4 * np.pi * dens * val2)) + bias

इससे तुरंत ही सीमित अंतराल के भीतर r_cap पर निर्भरता गायब हो जाती है, और जहाँ डोमेन पर समाकल पर्याप्त सटीक है, t=0 का मान सैद्धांतिक अभिव्यक्ति से मेल बैठता है।

[0, ∞) पर समाकलों तक विस्तार

आधारभूत लाइब्रेरी में अनंत सीमाओं के साथ points आर्ग्यूमेंट काम नहीं करता। इसका आसान उपाय यह है कि समाकल को दो हिस्सों में बाँट दें: एक सीमित खंड, जहाँ points दिए जाएँ, और एक टेल, जो points के बिना अनंत तक जाए। व्यवहार में कटऑफ वहीं चुनें जहाँ टेल नगण्य हो जाए या कम-से-कम चिकनी हो।

def hopfield_fn(t, w0, r0, dens, amp, bias):
    r_cap = 1.0
    anchors = np.array([1e-9, 1e-8, 1e-7])
    a_finite = quad(kern_a, 0, r_cap, args=(w0, r0, t), points=anchors, epsabs=0, epsrel=1.49e-8)[0]
    a_tail   = quad(kern_a, r_cap, np.inf, args=(w0, r0, t), epsabs=0, epsrel=1.49e-8)[0]
    b_finite = quad(kern_b, 0, r_cap, args=(w0, r0, t), points=anchors, epsabs=0, epsrel=1.49e-8)[0]
    b_tail   = quad(kern_b, r_cap, np.inf, args=(w0, r0, t), epsabs=0, epsrel=1.49e-8)[0]
    a_val = a_finite + a_tail
    b_val = b_finite + b_tail
    return amp * ((4 * np.pi * dens * a_val) * np.exp(4 * np.pi * dens * b_val)) + bias

इस बदलाव के साथ, [0, ∞) पर t=0 की सैद्धांतिक जाँच इस तरह मेल खाती है।

Theoretical value for t=0: 1.628601631620949e-39
Calculated value for t=0: 1.628601631620949e-39

np.vectorize और उसके विकल्पों के बारे में

समाकलन की समस्या की जड़ np.vectorize नहीं है। यह बस एक सुविधाजनक रैपर है जो t पर लूप लगाकर परिणाम इकट्ठा करता है। एक साधारण Python लूप समकक्ष है और अक्सर डिबग करना आसान रहता है।

grid_t = np.linspace(0, 1e8, 1000)
vals = np.array([hopfield_fn(tt, w0, r0, dens, amp, bias) for tt in grid_t])

quad स्वयं एक स्केलर इंटीग्रेटर है: इसका इंटीग्रैंड स्केलर r पर ही आंका जाता है, इसलिए t को एक-एक करके पास करना ही सही तरीका है।

वैकल्पिक: “दिलचस्प” r स्वतः ढूँढना

जहाँ उपयुक्त एंकर-पॉइंट्स पहले से अंदाज़ा लगाना आसान नहीं, वहाँ एक प्रोग्रामेटिक तरीका यह है कि इंटीग्रैंड्स के गैर-तुच्छ स्थिरांक बिंदु खोजे जाएँ और उन्हीं से डोमेन विभाजित किया जाए। नीचे का प्रोग्राम प्रत्येक कर्नेल का प्रतीकात्मक अवकलज निकालता है, r=0 से दूर प्रथम अवकलज का शून्य ढूँढता है, और उसी बिंदु के आसपास एंकर रखता है। फिर इन्हीं एंकरों के साथ कड़ी सहनशीलताओं में quad चलाया जाता है। यह तरीका कुछ अधिक जटिल है और यहाँ एक स्वचालित विकल्प के रूप में दिखाया गया है।

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import quad
import sympy
import mpmath
import scipy
import math
from collections import defaultdict
mpmath.mp.dps = 50
def rate_fn(r, w0, r0):
    return w0 * np.exp((-2 * r) / r0)
def kern_a(r, w0, r0, t):
    return rate_fn(r, w0, r0) * np.exp(-rate_fn(r, w0, r0) * t) * r**2
def kern_b(r, w0, r0, t):
    return (np.exp(-rate_fn(r, w0, r0) * t) - 1) * r**2
w0_symb = sympy.Symbol('w0')
r0_symb = sympy.Symbol('r0')
t_symb = sympy.Symbol('t')
r_symb = sympy.Symbol('r')
w0c = 1e6
r0c = 1.2e-9
def rate_expr():
    return w0_symb * sympy.exp((-2 * r_symb) / r0_symb)
def kern_a_expr():
    return rate_expr() * sympy.exp(-rate_expr() * t_symb) * r_symb**2
def kern_b_expr():
    return (sympy.exp(-rate_expr() * t_symb) - 1) * r_symb**2
def expr_derivative(expr, w0, r0, n=1):
    expr = expr.subs(w0_symb, w0)
    expr = expr.subs(r0_symb, r0)
    deriv = sympy.diff(expr, r_symb, n)
    deriv = sympy.simplify(deriv)
    return sympy.lambdify((r_symb, t_symb), deriv, modules=['mpmath'])
def expand_bracket(func, cond, args, start, mult, maxiter):
    x = start
    if cond(func(x, *args)):
        raise Exception('invalid start point')
    for _ in range(maxiter):
        last_x = x
        x *= mult
        if cond(func(x, *args)):
            return last_x, x
    raise Exception(f'Unable to find starting area for {func} with {args}')
def safe_underflow_wrapper(func):
    def inner(x, *args):
        y_mpf = func(x, *args)
        y_f = float(y_mpf)
        if y_f == 0 and y_mpf != 0:
            y_f = float(mpmath.sign(y_mpf)) * math.ulp(0)
        return y_f
    return inner
def root_after(function, derivative, t, lower_bound):
    t = np.clip(t, 1e-30, np.inf)
    if derivative(lower_bound, t) > 0:
        cond = lambda x: x <= 0
    else:
        cond = lambda x: x > 0
    lo, hi = expand_bracket(derivative, cond, (t,), lower_bound, 5, 100)
    rng = hi - lo
    xtol = rng * 1e-3
    res = scipy.optimize.bisect(
        safe_underflow_wrapper(derivative), lo, hi, args=(t,), xtol=xtol
    )
    return res
ka_d0 = expr_derivative(kern_a_expr(), w0c, r0c, n=0)
ka_d1 = expr_derivative(kern_a_expr(), w0c, r0c, n=1)
kb_d0 = expr_derivative(kern_b_expr(), w0c, r0c, n=0)
kb_d1 = expr_derivative(kern_b_expr(), w0c, r0c, n=1)
last_root_a = None
last_root_b = None
def pick_points_a(t):
    r_star = root_after(ka_d0, ka_d1, t, lower_bound=1e-15)
    return [r_star * 0.01, r_star * 0.9, r_star, r_star * 1.3, r_star * 20]
def pick_points_b(t):
    r_star = root_after(kb_d0, kb_d1, t, lower_bound=1e-15)
    return [r_star * 0.01, r_star * 0.03, r_star, r_star * 1.4, r_star * 20]
def hopfield_fn_auto(t, w0, r0, dens, amp, bias):
    r_cap = 1.0
    pts_a = pick_points_a(t)
    pts_b = pick_points_b(t)
    a_val, _ = quad(kern_a, 0, r_cap, args=(w0, r0, t), points=pts_a, epsabs=0, epsrel=1.49e-8)
    b_val, _ = quad(kern_b, 0, r_cap, args=(w0, r0, t), points=pts_b, epsabs=0, epsrel=1.49e-8)
    return amp * ((4 * np.pi * dens * a_val) * np.exp(4 * np.pi * dens * b_val)) + bias
dens = 3e-19
amp = 1
bias = 0
teor_t0 = dens * np.pi * w0c * r0c**3
calc_t0 = hopfield_fn_auto(0, w0c, r0c, dens, amp, bias)
print('Theoretical value for t=0:', teor_t0)
print('Calculated value for t=0:', calc_t0)
T = np.linspace(0, 100000000, 1000)
Y = np.vectorize(hopfield_fn_auto)(T, w0c, r0c, dens, amp, bias)
plt.plot(T, Y, color='r', linestyle='-', linewidth=0.7)
plt.yscale('log')
plt.xscale('log')
plt.show()

यह जानना क्यों उपयोगी है

अनुकूली क्वाड्रेचर शक्तिशाली है, सर्वज्ञ नहीं। यदि आप सैम्पलर को पैमाने-संवेदनशील क्षेत्रों की ओर नहीं ले जाते या सहनशीलताएँ गणना के परिमाणों से मेल नहीं खातीं, तो कटऑफ पैरामीटर पर दिखने वाली मनमानी निर्भरताएँ उभर आती हैं। points और tolerances का सही उपयोग आपको भौतिकी को बरकरार रखते हुए कृत्रिम कलाकृतियों को हटाने देता है।

व्यावहारिक सार

Thomas–Hopfield जैसे मॉडलों के कर्नेल्स का समाकल करते समय, points देकर इंटीग्रेटर को वहीं देखने को कहें जहाँ संरचना है। बहुत छोटे परिमाणों के लिए epsabs को शून्य करके कड़ा करें और सापेक्ष सहनशीलता पर भरोसा रखें। यदि 0 से अनंत तक का सही अपूर्ण समाकल चाहिए, तो उसे दो हिस्सों में बाँटें—एक सीमित भाग points के साथ, और एक अनंत टेल points के बिना। np.vectorize का उपयोग ठीक है; एक साधारण Python लूप भी उतना ही मान्य है और कई बार अधिक पारदर्शी। अंत में, इसी तरह के प्रसंगों में छोटे x के लिए exp(x) - 1 जैसी अभिव्यक्तियाँ निकालते समय catastrophic cancellation से बचने के लिए np.expm1(x) पर विचार करें।

इन समायोजनों के साथ आपके प्लॉट स्थिर हो जाते हैं, t=0 की जाँच सैद्धांतिक मान से मेल खाती है, और आप इंटीग्रेटर से जूझने के बजाय निश्चिंत होकर मॉडल फिट करने की ओर बढ़ सकते हैं।