CODEWITHSUNDEEP
pythonnumpy
ari
ari

Reputation: 4513

Take laplacian of unevenly sampled data in python

I have a bunch of unevenly sampled (1-D) data in the form of x and y values that I would like to take the laplacian of. Is there an easy way to do this in numpy/scipy?

I was thinking of taking the gradient of the gradient, but wouldn't that introduce artifacts, even in evenly sampled data? eg, given the data 0,0,0,4,8,8,8 and 1d laplace operator 1,-2,1, the laplacian would be 0,0,4,0,-4,0,0, but using the gradient of the gradient would yield 0,1,2,0,-2,-1,0

Upvotes: 1

Views: 790

Answers (1)

Jaime
Jaime

Reputation: 67427

One way of arriving at the [1, -2, 1] operator for evenly spaced grids is to either compute the first derivatives with a forward difference scheme, and the second with a backwards difference scheme, or viceversa, since both yield the same result. So a good option for your uneven grid is to do it both ways and then average the results.

Lets say we are looking at a grid point which is hb ahead of the previous grid point, and hf behind the next one:

--+----+----+---
 x-hb  x   x+hb

If I didn't mess my algebra, your Laplacian at x computed as the average described above would be:

(f(x+hf)*(1+hb/hf) + f(x)(2+hb/hf+hf/hb) + f(x-hb)(1+hf/hb)) / (2*hf*hb)

You can compute this on a 1D array doing:

def laplacian(x, y):
    x = np.asarray(x)
    y = np.asarray(y)
    hb = x[1:-1] - x[:-2]
    hf = x[2:] - x[1:-1]
    y_hb = y[:-2]
    y_hf = y[2:]
    hb_hf = hb / hf
    hf_hb = hf / hb
    return (y_hf*(1+hb_hf) - y[1:-1]*(2+hb_hf+hf_hb) +
            y_hb*(1+hf_hb)) / 2 / hb / hf

>>> laplacian(range(7), [0,0,0,4,8,8,8])
array([ 0.,  4.,  0., -4.,  0.])

Upvotes: 1

Related Questions