Numerical approximation of the 2nd derivative across a diffuse interface

[Dilon]

Imagine you create a diffuse interface in space and determine which side of the interface you are on by a local scalar value that can be between 0 and 1. We could create a circle, centered in a rectangular ynum-by-xnum grid, with such a diffuse interface with the following MATLAB code:

xnum = 100; % grid nodes in x direction
ynum = 100; % grid nodes in y direction

xgrid       =  (-(xnum/2-1/2):1:(xnum/2-1/2)) ;  x = repmat(xgrid,[ynum 1]); % grid distance in x direction
ygrid       =  (-(ynum/2-1/2):1:(ynum/2-1/2))';  y = repmat(ygrid,[1 xnum]); % grid distance in y direction
r =  sqrt(x.^2+y.^2); % distance from central node

L = 10; % thickness of diffuse interface
radii = 40; % radius of circle

F = 1-(1+sin((pi/L)*(r-radii)))/2; % create diffuse interface
F(r<(radii-L/2))       =  1; % create diffuse interface
F(r>(radii+L/2))       =  0; % create diffuse interface

What I want to do is, use numerical methods on the regular grid to get the best approximation of the second derivative of F in the direction normal to the interface, for every node within the interface. To start, I get surface normal components as

dfdy       =  abs(circshift(F,[1 0])-circshift(F,[-1 0]))/2; % first derivative of f in y direction
dfdx       =  abs(circshift(F,[0 1])-circshift(F,[0 -1]))/2; % first derivative of f in x direction
normy = dfdy./sqrt(dfdy.^2+dfdx.^2); % surface normal, y component
normx = dfdx./sqrt(dfdy.^2+dfdx.^2); % surface normal, x component

However, even with this information, I do not see an easy way to compute the derivatives across the interface in the normal direction using the gridded data. My best effort so far has been to project a line from each node, in the direction of the surface normal, until it intersects the edge of the local cell (see below). At this position the value of F is estimated by linear interpolation between nodes a and b. The same is done in the backward direction and using the known distance, the 2nd derivative at node i is found as

dfdn2 = (F0+F1-F*2)/d^2

This part of the code is a big loop so I just explain it in this illustration:

I've compared my numerical prediction to the true value for the second derivative across the interface, which is:

dfd2actual     = (1-2.*F).*((pi/L)^2/2)

My approach works okay, but it is not very precise, with errors that can be quite large (>5%).

Is there perhaps a better way to do this?? I planned to try cubic interpolation instead of linear interpolation to improve things, but maybe a different approach is better.

 

[Twigg]

Ok, I'm a little lost in the specifics, but I have a couple of diagnostic questions.

Why not do this with a polar coordinate grid? In this case, it would make the normal derivatives a heck of a lot simpler. Is this a test case for a bigger project that requires cartesian coordinates?

How does the error you reported (the ~5% one) scale with the mesh size? Vary the mesh size and see how your error scales. If this is a simple truncation error, it should scale down with a finer mesh, and you can brute force your way through with a higher mesh count at the cost of more RAM usage. (Note: by mesh size I'm referring to what you call xnum and ynum.)

 

[Stephen Tashi]

What I want to do is, use numerical methods on the regular grid to get the best approximation of the second derivative of F in the direction normal to the interface, for every node within the interface.

Your goal isn't clear. You illustrate taking a known function on the grid and approximating a set of directional second derivatives of that function. In that case "best approximation" is clearly defined because we can compare the approximation method to the actual directional second derivatives of the function.

However, it's possible that your example is intended as a test of approximation methods that will be used when you don't have a known function and only have data on grid points - or have a known function that is very time consuming to compute and intend only to compute it at grid points.

If your are developing an approximation for an abitrary unknown function, it isn't clear what "best approximation" means. There may not be a best approximation method or even a good approximation method that works for any arbitrary function. I think answering your question depends on narrowing down what kind of functions you want to approximate.