簡単な要約: 2 つの配列の有限畳み込みをすばやく計算するにはどうすればよいですか?
問題の説明
で定義される 2 つの関数 f(x)、g(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。この関数は、ここでconv提案されているアルゴリズムと同じ結果をもたらします。ただし、結果は分析ソリューションとはかなり異なります。台形規則を使用するようにアルゴリズムを変更すると、望ましい結果が得られます。conv
これを説明するために、
f(x) = exp(-x)
g(x) = 2 * exp(-2 * x)
結果は次のとおりです。
ここでRiemannは、単純なリーマン和を表し、trapezoidalは台形規則を使用するようにリーマン アルゴリズムを修正したバージョンであり、scipy.signal.convolveは scipy 関数でanalyticalあり、 は解析的畳み込みです。
結果は次のようg(x) = x^2 * exp(-x)になります。
ここで「比率」は、scipy から取得した値と分析値の比率です。上記は、積分を再正規化しても問題を解決できないことを示しています。
質問
scipy の速度を使用して、台形規則のより良い結果を保持することは可能ですか? または、目的の結果を得るために C 拡張機能を作成する必要がありますか?
例
以下のコードをコピーして貼り付けるだけで、発生している問題を確認できます。2 つの結果は、変数を大きくすることで、より一致させることができstepsます。積分が増加すると過大評価され、減少すると再び解析解に近づくため、問題は右手リーマン和からのアーティファクトによるものだと思います。
編集：関数と同じ結果を与える比較として、元のアルゴリズム2scipy.signal.convolveを含めました。
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()
お時間をいただきありがとうございます！
