# 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()