Strange scaling issue with derivative implementation

Hello,

I am trying to implement a second derivative (3d Laplacian) in MATLAB to solve differential equations. I am getting the following strange scaling issue.

As a test, I am running my code on the function psi=x.^2+y.^2+z.^2, which should return 6 identically on non-extreme cells, i.e., inside the 3d matrix. If I run it on the unit cube with a spatial resolution of .1, it works fine. If I run it on the unit cube with spatial resolution .01, it return .0006, as if it doesn't get multiplied by the right scaling factor. Here's the function:

function psi2=secondDerivative(psi,dr)
% compute second derivative of 3d matrix

% get size
s=size(psi);
scale=1./(dr.^2);

%compute each side's differentiation matrix
DiffMatx=diffMat(s(1));
DiffMaty=diffMat(s(2));
DiffMatz=diffMat(s(3));

% compute dimensional permutations
psiY=permute(psi,[2 1 3]);
psiZ=permute(psi,[3 1 2]);

psi2=scale(1).*(ThreeDMultiply(DiffMatx,psi))...
    +scale(2).*permute(ThreeDMultiply(DiffMaty,psiY),[2 1 3])...
    +scale(3).*permute(ThreeDMultiply(DiffMatz,psiZ),[2 3 1]);
end

You can see the scaling factor near the top. DiffMat and ThreeDMultiply compute a derivative matrix (3 point symmetric) and extend matrix multiplication to a rank2 times rank3, respectively. I attached them for your convenience, but I tested them extensively.

function D=diffMat(K)

% generate differentiation matrix of kxk elemenets
A=sparse(1:K,1:K,-2*ones(1,K),K,K);
B=sparse(1:K-1,2:K,ones(1,K-1),K,K);
C=sparse(2:K,1:K-1,ones(1,K-1),K,K);
D=A+B+C;
end

function B=ThreeDMultiply(M,A)
%%return B=MA, where
% M is a rank 2 tensor
% A is a rank 3 tensor
m=size(M);
a=size(A);

A2=reshape(A,[a(1),a(2)*a(3)]);
B=reshape(M*A2,a);

Here's what I do:

[x y z]=meshgrid(0:.1:1,0:.1:1,0:.1:1);
psi=x.^2+y.^2+z.^2;
dr=.1*ones(3,1);
out=secondDerivative(psi,dr);

That returns an 11x11x11 matrix such that all entries out_{i,j,k}=6 if none of i,j,k are 1, as expected. The following code, however

[x y z]=meshgrid(0:.01:1,0:.01:1,0:.01:1);
psi=x.^2+y.^2+z.^2;
dr=.01*ones(3,1);
out=secondDerivative(psi,dr);

returns a 101x101x101 matrix such that all entries out_{i,j,k}=0.0006, unless one of i,j,k are 1. This is off by a factor of 10000, i.e., it seems to ignore the multiplication by a scaling factor. I did make double-triple sure that I am using the correct dr for the correct meshgrid.

Thanks for your help!

Answer 1

[David Goodmanson]

Hi Nimrod,

probably a bit late for an answer, especially one that doesn't really change anything. I believe that your code is working correctly, and outputs 6 where you expect it to in both cases. Is it possible that in the second case you were misled by the forest of .0006's that you get when displaying e.g. out(23,:,:) and did not notice the factor of 1.0e+04 in front?