संक्षिप्त सारांश: मैं जल्दी से दो सरणी के सीमित संकल्प की गणना कैसे करूं?Riemann से Artefacts scipy.signal.convolve
समस्या का वर्णन
मैं के दो कार्य f (x), जी परिमित घुमाव के प्राप्त करने के लिए कोशिश कर रहा हूँ (एक्स)
द्वारा परिभाषित किया गया इस लक्ष्य को हासिल करने के लिए, मैं ले लिया है असतत नमूने कार्यों की और उन्हें लंबाई steps की सरणियों में बदल गया:
xarray = [x * i/steps for i in range(steps)]
farray = [f(x) for x in xarray]
garray = [g(x) for x in xarray]
मैं तो रूपा की गणना करने की कोशिश की scipy.signal.convolve फ़ंक्शन का उपयोग करके olution। यह फ़ंक्शन एक ही परिणाम देता है जैसे एल्गोरिदम conv सुझाए गए here। हालांकि, परिणाम विश्लेषणात्मक समाधान से काफी अलग हैं। Trapezoidal नियम का उपयोग करने के लिए एल्गोरिदम conv संशोधित वांछित परिणाम देता है।
इसे दर्शाने के लिए, मैं
f(x) = exp(-x)
g(x) = 2 * exp(-2 * x)
जाने परिणाम हैं:
यहाँ Riemann एक सरल Riemann कुल योग होता है, trapezoidal Riemann एल्गोरिथ्म का उपयोग करने का एक संशोधित संस्करण है trapezoidal नियम, scipy.signal.convolve scipy समारोह है और analytical विश्लेषणात्मक दृढ़ संकल्प है।
अब g(x) = x^2 * exp(-x) जाने और परिणाम हो:
यहाँ 'अनुपात' scipy से विश्लेषणात्मक मूल्यों के लिए प्राप्त मूल्यों का अनुपात है। उपर्युक्त दर्शाता है कि अभिन्न अंगों को पुनर्निर्मित करके समस्या हल नहीं की जा सकती है।
सवाल
यह scipy की गति का उपयोग, लेकिन एक समलम्बाकार शासन के बेहतर परिणाम को कायम रखने में मैं एक सी विस्तार वांछित परिणाम प्राप्त करने के लिए लिखने के लिए क्या है करना संभव है?
एक उदाहरण
की प्रतिलिपि बना समस्या मैं का सामना कर रहा हूँ देखने के लिए नीचे कोड पेस्ट करें। steps परिवर्तनीय को बढ़ाकर दो परिणामों को निकट समझौते में लाया जा सकता है। मेरा मानना है कि समस्या दाएं हाथ से कलाकृतियों के कारण है, रिमेंन रकम क्योंकि अभिन्न अंग बढ़ रहा है जब यह बढ़ रहा है और विश्लेषणात्मक समाधान के रूप में फिर से घट रहा है।
संपादित: मैं अब एक तुलना जो scipy.signal.convolve समारोह के रूप में एक ही परिणाम देता है के रूप में मूल एल्गोरिथ्म 2 शामिल किया है।
import numpy as np
import scipy.signal as signal
import matplotlib.pyplot as plt
import math
def convolveoriginal(x, y):
'''
The original algorithm from http://www.physics.rutgers.edu/~masud/computing/WPark_recipes_in_python.html.
'''
P, Q, N = len(x), len(y), len(x) + len(y) - 1
z = []
for k in range(N):
t, lower, upper = 0, max(0, k - (Q - 1)), min(P - 1, k)
for i in range(lower, upper + 1):
t = t + x[i] * y[k - i]
z.append(t)
return np.array(z) #Modified to include conversion to numpy array
def convolve(y1, y2, dx = None):
'''
Compute the finite convolution of two signals of equal length.
@param y1: First signal.
@param y2: Second signal.
@param dx: [optional] Integration step width.
@note: Based on the algorithm at http://www.physics.rutgers.edu/~masud/computing/WPark_recipes_in_python.html.
'''
P = len(y1) #Determine the length of the signal
z = [] #Create a list of convolution values
for k in range(P):
t = 0
lower = max(0, k - (P - 1))
upper = min(P - 1, k)
for i in range(lower, upper):
t += (y1[i] * y2[k - i] + y1[i + 1] * y2[k - (i + 1)])/2
z.append(t)
z = np.array(z) #Convert to a numpy array
if dx != None: #Is a step width specified?
z *= dx
return z
steps = 50 #Number of integration steps
maxtime = 5 #Maximum time
dt = float(maxtime)/steps #Obtain the width of a time step
time = [dt * i for i in range (steps)] #Create an array of times
exp1 = [math.exp(-t) for t in time] #Create an array of function values
exp2 = [2 * math.exp(-2 * t) for t in time]
#Calculate the analytical expression
analytical = [2 * math.exp(-2 * t) * (-1 + math.exp(t)) for t in time]
#Calculate the trapezoidal convolution
trapezoidal = convolve(exp1, exp2, dt)
#Calculate the scipy convolution
sci = signal.convolve(exp1, exp2, mode = 'full')
#Slice the first half to obtain the causal convolution and multiply by dt
#to account for the step width
sci = sci[0:steps] * dt
#Calculate the convolution using the original Riemann sum algorithm
riemann = convolveoriginal(exp1, exp2)
riemann = riemann[0:steps] * dt
#Plot
plt.plot(time, analytical, label = 'analytical')
plt.plot(time, trapezoidal, 'o', label = 'trapezoidal')
plt.plot(time, riemann, 'o', label = 'Riemann')
plt.plot(time, sci, '.', label = 'scipy.signal.convolve')
plt.legend()
plt.show()
आपके समय के लिए धन्यवाद!
अगर आप प्रदान करते हैं यह उपयोगी हो सकता है एक [ पूर्ण न्यूनतम उदाहरण] (http: // sscce।संगठन /) जो समस्या को पुन: उत्पन्न करता है, छोटी त्रुटियों को बाहर करने के लिए जैसे कि पूर्णांक विभाजन का उपयोग किया जाता है जहां सच्चे विभाजन का उपयोग किया जाना चाहिए। – jfs
यदि आप 'विज्ञान' सरणी को दाईं ओर (एक चरण से) स्थानांतरित करते हैं और इसे सामान्यीकृत करते हैं तो समाधान [समान दिखते हैं] (http://i403.photobucket.com/albums/pp111/uber_ulrich/convolve.png): [' sci = np.r_ [0, विज्ञान [: चरण -1]] * 0.86'] (https://gist.github.com/ac45fb5d117d8ecb66a3)। ऐसा लगता है कि 'convolve() 'साधन (एफएफटी, परिपत्र) के विभिन्न परिभाषाएं हैं, [नम्पी/Scipy में कनवॉल्यूशन computations] में चर्चा देखें (http://stackoverflow.com/q/6855169/4279)। – jfs
संकल्प धागे के संदर्भों के लिए धन्यवाद। सही बदलाव समस्या को ठीक करने के एक दिलचस्प तरीके की तरह दिखता है। क्या आप मुझे बता सकते हैं कि आपने 0.86 के सामान्यीकरण कारक को कैसे प्राप्त किया? मैंने मूल रिमैन योग एल्गोरिदम को शामिल किया है ताकि यह समझाया जा सके कि मेरा मानना है कि यह संकल्प का अर्थ क्या है, इसकी एक अलग परिभाषा के बजाय यह संख्यात्मक आर्टेफैक्ट है। –