Finite Difference Matrix for 1D and 2D problem.

Central finite difference matrix for the estimation of n-th derivative of function f, i.e., d^n f/dx^n with arbitrary order of accuracy.

The code was written in 2017 when I was in Berlin. The purpose was to set up numerical equations for solving partial differential equations using finite difference. The original package includes some functionalities, such as grid transformation and matrix-free finite difference schemes, which may require customised parameter-tuning for different problems so I am not yet brave enough to publish them. The functions uploaded here only generate the "classic" finite difference matrices with (possibly) arbitrary order of accuracy. They should be enough to solve some basic boundary-value problems, such as the 1D and 2D Helmholtz equations.

The user can use
* getNonCompactFDmatrix
* getNonCompactFDmatrix2D
to generate central finite difference matrix for 1D and 2D problems, respectively. The user needs to specify
1, number of points
2, spatial step
3, order of derivative
4, the order of accuracy (an even number) of the finite difference scheme.
The stencils at the boundary are non-symmetric but have the same order of accuracy as the central finite difference.

The core function used by the above functions to generate the finite difference weightings is getNonCompactFDMWeights, which works for arbitrary stencils (symmetric and non-symmetric stencils). The mathematical basics of code is documented in the attached readme.pdf. The user may use this core function to generate single-sided finite difference matrices.

Two example files are enclosed. Here is another simpler example:

fun = @(x) sin(pi*x)+exp(2*x); % Test function
d2fundx = @(x) -pi^2*sin(pi*x)+4*exp(2*x); % 2nd derivative

% Prepare for finite difference matrix (FDM):
npx = 23; % Number of points
xVec = linspace(-1,1.5,npx).'; % x vector
dx = diff(xVec([1 2])); % Spatial step
n = 2; % derivative order
ooa = 6; % order of accuracy of the FDM
D2x = getNonCompactFDmatrix(npx, dx, n, ooa); % Here is the FDM

plot(xVec, d2fundx(xVec), xVec, D2x*fun(xVec), 'o')
xlabel('x'); ylabel(func2str(d2fundx))
legend('Exact', sprintf('Numerical. %d derivative.', n))