回答此问题可获得 20 贡献值,回答如果被采纳可获得 50 分。
<p><strong>简短摘要</strong>:如何快速计算两个数组的有限卷积?在</p>
<h2>问题描述</h2>
<p>我试图得到由定义的两个函数f(x),g(x)的有限卷积</p>
<p><img src="https://i.stack.imgur.com/wqGyY.png" alt="finite convolution"/></p>
<p>为了实现这一点,我对函数进行了离散采样,并将它们转换为长度为<code>steps</code>的数组:</p>
<pre><code>xarray = [x * i / steps for i in range(steps)]
farray = [f(x) for x in xarray]
garray = [g(x) for x in xarray]
</code></pre>
<p>然后我尝试使用<code>scipy.signal.convolve</code>函数计算卷积。此函数给出的结果与<code>conv</code>建议的算法<a href="http://www.physics.rutgers.edu/~masud/computing/WPark_recipes_in_python.html" rel="nofollow noreferrer">here</a>相同。然而,结果与解析解有很大不同。修改算法<code>conv</code>以使用梯形规则,可以得到预期的结果。在</p>
<p>为了说明这一点,我让</p>
^{pr2}$
<p>结果是:</p>
<p><img src="https://i.imgur.com/u5Kle.png" alt="enter image description here"/></p>
<p>这里<code>Riemann</code>表示一个简单的Riemann和,<code>trapezoidal</code>是使用梯形规则的Riemann算法的改进版本,<code>scipy.signal.convolve</code>是scipy函数,<code>analytical</code>是解析卷积。在</p>
<p>现在让<code>g(x) = x^2 * exp(-x)</code>得到的结果是:</p>
<p><img src="https://i.stack.imgur.com/4Xntp.png" alt="enter image description here"/></p>
<p>这里的“比率”是从scipy获得的值与分析值的比率。上述结果表明,该问题不能通过重新规范化积分来解决。在</p>
<h2>问题</h2>
<p>是否可以使用scipy的速度,但保留梯形规则的更好结果,还是必须编写一个C扩展来获得所需的结果?在</p>
<h2>一个例子</h2>
<p>只需复制并粘贴下面的代码,以查看我遇到的问题。增加<code>steps</code>变量可以使这两个结果更接近一致。我认为这个问题是由于右Riemann和的伪影造成的,因为积分在增大时被高估了,当积分减小时又接近解析解。在</p>
<p><strong>编辑</strong>:我现在将原始算法<a href="http://www.physics.rutgers.edu/~masud/computing/WPark_recipes_in_python.html" rel="nofollow noreferrer">2</a>作为比较,其结果与<code>scipy.signal.convolve</code>函数相同。在</p>
<pre><code>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.<a href="https://www.cnpython.com/list/append" class="inner-link">append</a>(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()
</code></pre>
<p>感谢您抽出时间!在</p>