Fast way to calculate source response using Green's function I am a newbie in the scientific computing in python. The function I am trying to calculate is the response of the wave equation for a given source term (For reference, See: Equation 11.67 in https://webhome.phy.duke.edu/~rgb/Class/phy319/phy319/node75.html). 

$$\phi(\pmb{x},t) = \int_V G_+(\pmb{x}, t; \pmb{x'}, t)\rho(\pmb{x'}, t')d^3x'dt'$$

To give context, the term phi in the equation 11.67 is the displacement and the term rho can be thought as a source of disturbance. 

Now, in my case, the problem is constructed in spatial dimension of 2 (x-y). Thus, I have to iterate the equation for grid points in x, y and t. This makes the overall calculation extremely time-consuming. Below is my lines of code:

import numpy as np

# Assigning grid points

Nx = 100
Ny = 100
Nt = 100

x = np.linspace(-Nx, Nx, 100)
y = np.linspace(-Ny, Ny, 100)
t = np.linspace(0, Nt-1, 100)

# Equation to solve
# phi(x,y,t) = ∫∫∫ G(x,y,t; x',y',t') . Source(x',y',t') dx' dy' dt'
# G(x,y,t; x',y',t') = Green's Function

# phi = displacement by the wave
phi = np.zeros((Nt,Ny,Nx))


# Define Function to realize Green's Function for Wave Equation
def gw(xx, yy, tt, mm, nn, ss):

    numer = (tt-ss) - np.sqrt((xx-mm)**2+(yy-nn)**2)
    denom = ((tt-ss)**2-(xx-mm)**2-(yy-nn)**2)

    if denom &lt; 0:
        denom = 0
    else:
        denom = np.sqrt(denom)

    kk = np.heaviside(numer,1)/(2*np.pi*denom+1)
    return (kk)

# Define Function to realize Gaussian Disturbance
def g_source(xx, yy, tt):

    spatial_sigma  = 5   # spatial width of the Gaussian 
    temporal_sigma = 30  # temporal width of the Gaussian 
    onset_time     = 20  # time when Gaussian Disturbance reaches its peak

    kk = np.exp((-(np.sqrt((xx)**2+(yy)**2)/spatial_sigma)**2)))*np.exp(-((tt-onset_time)/temporal_sigma)**2)
    return (kk)


for k in range (Nt):
    for j in range (Ny):
        for i in range (Nx):

            # this is the Green's function evaluated in each grid points                
            green = np.array([gw(x[i],y[j],t[k],i_grid,j_grid,k_grid) for k_grid in t for j_grid in y for i_grid in x])
            green = green.reshape(Nt, Ny, Nx)

            # this is the source term (rho), here it is a Gaussian Disturbance
            gauss = np.array([g_source(i_grid,j_grid,k_grid) for k_grid in t for j_grid in y for i_grid in x])
            gauss = gauss.reshape(Nt, Ny, Nx)


            # performing integration
            phi[k,j,i] = sum(sum(sum(np.multiply(green,gauss))))


Please let me know if there is a faster way to do this integration. Thanks in advance. I deeply appreciate your patience and support.
Performance wise

In general:

You use numpy, but you write it almost like in fortran. Python with numpy is good for scientific programing and computing as long as you don't do many loops


read this first A beginners guide to using Python for performance computing
if you really need to do tight loops than use
cython.
But most of the time you can avoid that by expressing the loop as some functional operation on the whole array 


most numpy functions can operator on whole array in one pass, use it.


always prefer y[:]=np.sqrt(x[:]) and avoid




          `for i in xrange(N): 
                y[i]=np.sqrt(x[i])`



operation like dot,convolve,einsum,ufunc.reduce are very potent to express most sciencetific algorithms in more abstract way
use advanced array slicing


you can express most inner if by boolen mask arrays or integer index array


e.g. f[ f&gt;15.0 ] = 15.0 clamp all elements in array f which are &gt;15.0 to 15.0; you can also store the boolean mask mask = (f&gt;15.0) and than use it like f[mask]+=g[mask]. This way you can express branching as fully functional programs/expressions with arrays.
Do not construct new np.array (np.zeros, np.ones etc.) too often (i.e. inside tight loops). Optimal performance is obtained if you prepare all arrays you need at the beggining. Avoid frequent conversion between list and array




To address your code example in particular:


Green's functions are basically convolutions. I'm pretty sure you can express it using e.g. scipy.ndimage.filters.convolve
if your convolution kernel is large (i.e. pixels interact more than with few neighbors) than it is often much faster to do it in Fourier space (convolution transforms as multiplication) and using np.fftn with O(nlog(N)) cost.
cannot def gw(xx, yy, tt, mm, nn, ss): operate on the whole array,
rather than individual numbers?
if denom &lt; 0: can be expressed using boolean mask array like mask = denom &gt;0; denom[mask] = np.sqrt(denom[maks])
Never do this if you can:


  for k in range (Nt):
     for j in range (Ny):
        for i in range (Nx):


the operations seems to be possible to rewrite that operate on the whole x,y space at once   


This is how you do convolution-like operations in numpy (from here)


           # The actual iteration
           u[1:-1, 1:-1] = ((u[0:-2, 1:-1] + u[2:, 1:-1])*dy2 +
                            (u[1:-1,0:-2] + u[1:-1, 2:])*dx2)*dnr_inv



this is so horrible:


green = np.array([gw(x[i],y[j],t[k],i_grid,j_grid,k_grid) for k_grid in t for j_grid in y for i_grid in x])



list comprehesion is relatively fast, but still much slower than numpy array operations (which are implemented in C).
do not create temporary list and convert it to temporary array, you loose lot time doing that.
why you cannot just prealocate Green[Nx,Ny,:] and  Gauss[Nx,Ny,:] and use it as a whole?