Hello,
I am working on a summer project at my university which involves simulating a pair of reaction-diffusion equations. I have chosen to approximate the differentials with the finite difference method. I have successfully applied the method to the laplacian operator (shown below), however, I appear to be having issue with the bi-harmonic operator (shown below).
Laplacian_2D_1stOrder
function out = Laplacian2D(in, res)
h = res;
k = res;
f_xx = (circshift(in,[ 0, 1]) - 2 .* circshift(in,[ 0, 0]) + circshift(in,[ 0, -1])).*h;
f_yy = (circshift(in,[ 1, 0]) - 2 .* circshift(in,[ 0, 0]) + circshift(in,[ -1, 0])).*k;
out = f_xx + f_yy;
end
BiHarmonic_2D_1stOrder
function out = Biharmonic2D(in, res)
h = res;
k = res;
f_xxxx = (circshift(in,[ 0, -2]) - 4.*circshift(in,[ 0, -1]) + 6.*circshift(in,[ 0, 0]) - 4.*circshift(in,[ 0, 1]) + circshift(in,[ 0, 2])).*res;
f_yyyy = (circshift(in,[ -2, 0]) - 4.*circshift(in,[ -1, 0]) + 6.*circshift(in,[ 0, 0]) - 4.*circshift(in,[ 1, 0]) + circshift(in,[ 2, 0])).*res;
f_xxyy = (circshift(in,[ -1, -1]) - 2.*circshift(in,[ 0, -1]) + circshift(in,[ 1, -1]) - 2.*circshift(in,[ -1, 0]) + 4.*circshift(in,[ 0, 0]) - 2.*circshift(in,[ 1, 0]) + circshift(in,[ -1, 1]) - 2.*circshift(in,[ 0, 1]) + circshift(in,[ 1, 1])).*res;
out = f_xxxx + f_xxyy + f_xxyy + f_yyyy;
end
Quoting from Wikipedia:
In mathematics, finite-difference methods (FDM) are numerical methods for solving differential equations by approximating them with difference equations, in which finite differences approximate the derivatives. FDMs are thus discretization methods.
The x and y resolutions are the same. I found that different resolutions led to spatial stretching. (I think this is obvious but would be happy to learn if there is a reason why or a situation where one would have different x and y spatial resolutions?)
When I call my functions, they appear to work, but the Laplacian appears far better behaved than the bi-harmonic operator. I tested both on the MATLAB Peaks function and compared them to MATHEMATICA's built in laplacian and hiharmonic operator functions and they returned the same results (roughly, I assume the difference is between my approximation and MATHEMATICA's more accurate differentiation).
My questions to the community are:
- Have I implemented these correctly, specifically the resolutions, but also the addition at the end of the functions to return the overall approximation?
- Is there a more efficient way to implement them within MATLAB? I know that MATLAB has a laplacian operator but I was unable to find a Bi-Harmonic function, hence deriving it myself.
- When applying these to PDEs, is the error order acceptable? The PDEs will be applied to a pair of Cahn-Hilliard equations:
Should they therefore be of a higher order for a stability or a convergence reason?
Thank you for your assisstance. Regards, Matthew
