tl;dr

Using a Taylor-matched method to find coefficients for the discretized equation $$\mathbf{A} \vec{f}'' = \mathbf{B} \vec{f}$$, a Fortran code has been implemented to find the second derivative and tested for $$f=\sin$$, with bad results on non-uniform grids. Which is the cause for that?

Methodology

I am new to non-uniform grids. Using the concepts described in the following papers:

I started implementing a Fortran module to solve differential equations on non-uniform grids. The idea behind these method is to discretize the (second-order, in my case) derivative with the following matricial equation: $$\mathbf{A} \vec{f}'' = \mathbf{B} \vec{f}$$. The matrices are defined as multi-diagonal, where the number of secondary diagonals reflects the formal error of the method. For example, in Gamet the general form $$\alpha_i f''_{i-1} + f''_i + \beta_i f''_{i+1} = A_i f_{i+1} + B_i f_{i-1} + C_i f _{i+2} + D_i f_{i-2} + E_i f_i$$ is used and the matrix $$\mathbf{A}$$ will be tridiagonal, while $$\mathbf{B}$$ will be pentadiagonal (in Lele, a similar form is defined, but with different symbols and concerning just uniform grids). $$f_i$$ and $$f''_i$$ can be expanded with a Taylor series on the grid points; as an example,

$$f''_{i+n} = f_i'' + nhf_i^{(3)} + \frac{(nh)^2}{2!}f_i^{(4)} + \mathcal{O}(h^3f^{(6)})$$

In order to determine the value of the coefficients, these expansions can be substituted in the above equation, then providing a system of equations matching terms of the same order (see below for how to solve this system).

The method is closed at the boundaries imposing the following forms:

$$i=1 \Rightarrow f''_1 + \alpha f''_2 = Af_1 + Bf_2 + Cf_3 + Df_4$$

$$i=2 \Rightarrow \alpha f_1'' + f_2'' + \beta f_3'' = A f_1 + Bf_2 + Cf_3 + Df_4$$

Strictly speaking, it follows that the matrices can have a higher number of diagonals at the boundary points.

In order to find the second derivative, the matrix $$\mathbf{A}$$ is inverted and $$\vec{f}'' = \mathbf{A}^{-1} \mathbf{B} \vec{f}$$. To do so, the Lapack library is used, as shown in http://fortranwiki.org/fortran/show/Matrix+inversion.

Grid

Different grids have been considered. The Fortran code which defines them is here reported (corresponing lines should be uncommented):

do ig=1,ng
!! uniform
!grid(ig,1) = a + (b-a)*dble(ig)/ng
!grid(ig,2) = (b-a)/ng
!! Shukla
!grid(ig,1) = (b+a)/2. + (b-a)*asin(-alpha*cos(PI_MATH*ig/ng))/(2.*asin(alpha))
!grid(ig,2) = (b-a)*( asin(-alpha*cos(PI_MATH*ig/ng))/(2.*asin(alpha)) &
!  - asin(-alpha*cos(PI_MATH*(ig-1)/ng))/(2.*asin(alpha)) )
!! random
!grid(ig,1) = a + (b-a)*(dble(ig)-1+rand())/ng
!if(ig.eq.1) then
!  grid(ig,2) = (b-a)/dble(ng)
!else
!  grid(ig,2) = grid(ig,1) - grid(ig-1,1)
!endif
end do


In particular, Shukla refers to the grid defined in the article by Shukla, Very high-order compact finite difference schemes on non-uniform grids for incompressible Navier–Stokes equations, 2007. The important thing to notice here is that this scheme has a parameter, alpha, which varies from 1 (uniform grid) to $$\alpha\rightarrow 0$$ (Chebyschev grid).

Testing

As a first test, the second derivative of a sine function has been evaluated and compared with the analytical second derivative:

$$x\in [0,1];\; f(x) = \sin(2\pi x); f''(x) = -(2\pi)^2\sin(2\pi x)$$

Using a uniform grid, the results are fine for Lele

For Gamet, using a uniform grid and a Shukla grid with alpha = 1 the results are approximately as good as with Lele (the two grids are of course the same; I compared them in order to be sure about the results):

Here is a comparison of the errors

Problems arise when the grid is imposed as non-uniform: for alpha = 0.9 and for alpha = $$10^{-10} \sim 0$$

Comparison of errors for different values of alpha:

It appears that the problem is due to the boundaries, but it is not clear to me why it happens. I double-checked the code in search of typos (but having solved the system with Python they should be ruled out) and I observe the same problems with a different method (cfr. Shukla 2007). I would like to know if this behaviour is to be expected and what I could do about that.

Solving the system for the coefficients

The cited papers provide sets of coefficients for different orders of the formal error. However, some situations could require to develop a custom method, for example setting some coefficients to zero, so here is a possible way to solve the system, using Sympy, a library for Python. The second argument of sym.solve is the list of variables for which to solve the system.

If a lower order is desired, it is sufficient to remove one of them (which can be then fixed to an arbitrary value, or in order to match a given method, e.g. for adapting Gamet's method to Lele's one) and comment a corresponding number of equations in the system (starting of course from the bottom one).

If some constraint is desired, a corresponding equation can be substituted in the system (again, starting from the bottom one).

import sympy as sym
a,b,c,d,e,alpha,beta = sym.symbols("a b c d e alpha beta")
# hmm = h_{i-2}; h0 = h_i; hp = h_{i+1}; etc.
hmm,hm,h0,hp,hpp   = sym.symbols("hmm hm h0 hp hpp")
solution = sym.solve((
a+b+c+d+e, # order -1
hp*a-h0*b+(hpp+hp)*c-(h0+hm)*d, # order 0
hp**2*a+h0**2*b+(hpp+hp)**2*c+(h0+hm)**2*d-factorial(2)/factorial(0)*(1+alpha+beta), # order 1
hp**3*a-h0**3*b+(hpp+hp)**3*c-(h0+hm)**3*d-factorial(3)/factorial(1)*(hp*beta-h0*alpha), # order 2
hp**4*a+h0**4*b+(hpp+hp)**4*c+(h0+hm)**4*d-factorial(4)/factorial(2)*(hp**2*beta+h0**2*alpha), # order 3
hp**5*a-h0**5*b+(hpp+hp)**5*c-(h0*hm)**5*d-factorial(5)/factorial(3)*(hp**3*beta-h0**3*alpha), # order 4
hp**6*a+h0**6*b+(hpp+hp)**6*c+(h0+hm)**6*d-factorial(6)/factorial(4)*(hp**4*beta+h0**4*alpha), # order 5
), (a,b,c,d,e,alpha,beta))
print("Gamet bulk:")
for sol in solution:
print(sol)
print(solution[sol])
print()