Polar Laplacian

Excuse my n00bness but I was hoping I could bet some pointers as to how to compute the Laplacian in polar coordinates. I'm using Trefethen's Chebyshev grid.

Below is how I compute the gradient. This seems to work fine

N = 31;
[D,r] = cheb(N); %function found at %https://people.maths.ox.ac.uk/trefethen/cheb.m
N2 = (N-1)/2;

M = 40;
dt = 2*pi/M; 
t = dt*(1:M)';
M2 = M/2;

[tt, rr] = meshgrid(t([M 1:M]), r(1:N2+1));

ff = rr.^3; %test function

[xx,yy] = pol2cart(tt,rr);

%first derivative d/dr+1/r*d/dt
[dudt, dudr] = gradient(ff, t([M 1:M]), r(1:N2+1));

V1 = dudr + dudt./rr;

clf, mesh(xx,yy,V1)

Now computing the Laplacian is proving more problematic. I've tried computing it from V1 using `gradient' again but I get this sort of `lip' near the boundary and it also messes up near the centre for more exotic functions.

[d2udt2, d2udtdr] = gradient((dudt./rr), t([M 1:M]), r(1:N2+1));
[d2udrdt, d2udr2] = gradient(dudr, t([M 1:M]), r(1:N2+1));

V2 = d2udr2 + dudr./rr + d2udt2./(rr.^2);

clf, mesh(xx, yy, V2)

Any help would be much appreciated. Thanks!

Answer 1

[Andrew Newell]

The Laplacian is close to correct except at the boundaries (r=0 and r=1). Try plotting, for example,

plot(rr(:,1),V2(:,1),rr(:,1),9*rr(:,1))

(the correct Laplacian being 9*r, of course). The problem is that gradient does a one-sided approximation to the derivative at the boundaries, and the errors are compounded when you take the gradient of the gradient. There doesn't seem to be a polar Laplacian in the File Exchange, so you may need to write your own on the lines of del2, which calculates centered second differences in the interior and then extrapolates to the boundaries.

EDIT: It occurred to me that the gradients might be more accurate if r were evenly spaced. It is better in the interior - but the outer lip folds over!

EDIT 2: If you happen to know that your solution depends on r only, you could use this little cheat to calculate the Laplacian:

d2udr2 = del2(ff, t([M 1:M]), r(1:N2+1))*4;
V2 = d2udr2 + dudr./rr;

EDIT 3: I only have time to sketch the more general approach you could use. Basically, del2 calculates the second derivative with respect to the first variable in the first loop, then the second derivative with respect to the second variable. It permutes the array so it can use the same equations for both. You could assign the results to d2udr2 and d2udt2 instead of adding them both to v. Cross-derivatives, unfortunately, don't fit in this framework. You'll need to add some more code that implements the finite difference formulae for the cross derivatives.

[Jarrod Rivituso]

I just want to say I found this question and answer very interesting :)

[Andrew Newell]

Thank you, Jarrod. Since you're a MATLAB insider, maybe you could push for support for higher order approximations in numerical derivatives.

[Tom Ashbee]

Thanks for the response Andrew. I figured the problem was something like that.

Are you suggesting I alter del2.m to make l~(1/y^2)*d^2/dx^2+d^2/dy^2 (and then get the (1/r)*d/dr term from the first calculation)? I've had a look at the del2 source code and can't really make any sense out of it (I'm not particularly experienced with MATLAB).

Wrt your edit, I've also tried a regular grid (I've tried many things!) and found the same thing. Also I do want to use a cheb gird in the rest of the program.

Also thanks for the comment Jarrod. I'm pretty surprised at how tricky this is proving.

[Tom Ashbee]

OK so after trying a few other things (e.g. computing the function as a sum of basis functions then obtaining the Laplacian from the eigenvalue equation and directly using Trefethen's differentiation matrix) I am back here.

Having looked at the del2 source code I am inclined to think I would be better off starting from scratch.

I'll post any progress I make here.

Thanks for the help.

Answer 2

[Andrew Newell]

If you can provide a function to calculate your potential (?) at any point, you can do much better. First, download John D'Errico's differentiation package. Then you can do the following:

N = 31;
r = cos(pi*(0:N)/N)'; % This is the relevant part of cheb(N)

N2 = (N-1)/2;
M = 40;
dt = 2*pi/M; 
t = dt*(1:M)';

[tt, rr] = meshgrid(t([M 1:M]), r(1:N2+1));
[xx,yy] = pol2cart(tt,rr);

Provide a function f(y), where y = [theta radius]:

f = @(y) y(2).^3; % This is your example function

Calculate all the derivatives:

grad = zeros(numel(tt),2);
hess = zeros(2,2,numel(tt));
for i=1:numel(tt)
  grad(i,:) = gradest(f,[tt(i) rr(i)]);
  hess(:,:,i) = hessian(f,[tt(i) rr(i)]);
end

This is the interpretation of the output:

dudt = reshape(grad(:,1),size(rr));
dudr = reshape(grad(:,2),size(rr));
d2udt2 = reshape(hess(1,1,:),size(rr));
d2udtdr = reshape(hess(1,2,:),size(rr));
d2udrdt = reshape(hess(2,1,:),size(rr));
d2udr2 = reshape(hess(2,2,:),size(rr));

Now you can plot the gradient and Laplacian:

V1 = dudr + dudt./rr;
mesh(xx,yy,V1)

V2 = d2udr2 + dudr./rr + d2udt2./(rr.^2);
figure, mesh(xx, yy, V2)

This is really accurate - often good to machine precision.

[Tom Ashbee]

Thanks very much for this Andrew.

I have one further query (that demonstrates my inexperience) about the function being differentiated. In the example I just chose f=r^3, but say I have an array of values of f already and want to compute their Laplacian. For example suppose I solved a PDE del^2u = f, then to check the solution I want compute the Laplacian of u and compare it to f.

When using del2 in a rectangular domain, one can just give the function a numerical array to differentiate, and something like that for this circular domain is what I am ultimately trying to do.

Thanks again for the help.

[Andrew Newell]

In that case you want to use the other approach I have described.