Instantly share code, notes, and snippets.

# tillahoffmann/conovleproblem.py Created Jan 17, 2012

What would you like to do?
Problem with numerical convolution
 import numpy as np import scipy.signal as signal import matplotlib.pyplot as plt import math import datetime from matplotlib.font_manager import FontProperties 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 = 7 #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] exp2 = [t ** 2 * math.exp(-t) for t in time] #Calculate the analytical expression #analytical = [2 * math.exp(-2 * t) * (-1 + math.exp(t)) for t in time] analytical = [1. / 3 * math.exp(-t) * t ** 3 for t in time] #Calculate the trapezoidal convolution trapezoidal = convolve(exp1, exp2, dt) #trapezoidal2 = convolve(exp1, sqexp, 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 sci = np.r_[0, sci[:steps - 1]] ratio = [a / b for a, b in zip(sci, analytical)] #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, sci, '.', label = 'signal.convolve') plt.plot(time, ratio, '.', label = 'ratio: signal.convolve/analytical') fontP = FontProperties() fontP.set_size('small') plt.legend(prop = fontP) plt.show()