Forgive my ignorance, but the function H is actually an infinite sum (including -n terms). I'm hoping to truncate at larger values on the order of the 100s ideally. My code seems to work but I am not sure if it is actually summing over the values of n described.
Code
import numpy as np
from scipy.integrate import trapz
tvals = [1, 2, 3, 4, 5] # fixed values of t
xvals = [0, 0.2, 0.4, 0.6, 0.8, 1.0] # fixed values of x
xi = np.linspace(0, 1, 100)
def H(xi, x, t):
for n in range(-2, 2):
return 0.5 * np.sin(np.pi * xi) * np.exp(-(x - 2 * n - xi)**2 / 4 * t) / np.sqrt(np.pi * t)
# Everything to the right of the sin term is part of a sum
# xi is the integral variable!
TrapzH = trapz(H(xi, xvals[1], tvals[0]), x=None, dx=0.1, axis=-1)
print(TrapzH)
I've tested this against a version where n = 1 and this has range (0, 1) and they seem to have different values but I am still unsure.
Your function H does not iterate over the specified n range since it exits the function on the first return encountered (at n=-2). Maybe you are looking for something like
def H(xi, x, t):
sum = np.zeros(xi.shape)
for n in range(-2, 2):
sum += np.exp(-(x - 2 * n - xi)**2 / 4 * t) / np.sqrt(np.pi * t)
return 0.5 * np.sin(np.pi * xi) * sum