Inverse of "diagonally not dominant matrix"

Question

I want to frame a higher order Central difference scheme of about $20^{th}$ order for first derivative. I'm using $20^{th}$ order because I need one scheme with good modified wave number. To find the co-efficient matrix (eg. for CD-2 it is [-1/2,1/2]), I have to solve $Ax=b$ equation. $x=A^{-1}*b$. I tried to solve that equation using

Build-in inverse (inv(A)) command of matlab,

"solve" command in R,

Gauss Seidel algorithm,

Gauss Jacobi

Conjugate gradient method

Gauss Seidel and Jacobi may fail because of the matrix is not diagonally dominant. Matlab gave answer with a warning and that answer is wrong. R- gave an error message because of poor conditional number.

Matrix $A$ is $20*20$ matrix :

Matrix $b$ is $20*1$ matrix

Is there any sophisticated algorithm or build-in command of any programming language to solve this? Please help me out.

Answer 1

One easy way to derive finite-difference approximations to derivatives is as follows. To find the coefficients $c_k$ corresponding to order-$s$ derivative on points $uh,(u+1)h,\ldots,vh$ ($u=-10,v=10,s=1$ in your case), write the condition defining, using Taylor theorem, as $$ \sum_{j=u}^{v}c_k e^{j h \partial} = \partial^s, $$ which should hold asymptotically as $h\to0$.

Introducing a change of variable $z=e^{h\partial}$, this is equivalent to $$z^u p(z) = h^{-s}(\log z)^s, $$ where $p$ is the polynomial $p(z) = c_u + c_{u+1}z + \cdots + c_v z^{v-u}$. Thus $$ p(z) = h^{-s}\frac{(\log z)^s}{z^u}. $$

So just compute the approximating Taylor polynomial around $z=1$ for $z^{-u}(\log z)^s$ to order $v-u$, this will be the polynomial with the finite-difference coefficients.

With $u=-10$, $v=10$, $s=1$, this polynomial is easily computed. The coefficients are: $$ h^{-1}\left\{\frac{1}{1847560},-\frac{5}{415701},\frac{5}{38896},-\frac{15}{17017},\frac{5}{1144},-\frac{12}{715},\frac{15}{286},-\frac{20}{143},\frac{15}{44},-\frac{10}{11},0,\frac{10}{11},-\frac{15}{44},\frac{20}{143},-\frac{15}{286},\frac{12}{715},-\frac{5}{1144},\frac{15}{17017},-\frac{5}{38896},\frac{5}{415701},-\frac{1}{1847560}\right\}. $$

Answer 2

I would not generally expect a "20th order" derivative estimate to typically be very stable/reliable/useful (e.g. due to well known artifacts of high-order polynomial interpolation).

That said, a general procedure for deriving finite-difference stencils is to solve an appropriate polynomial interpolation problem. I will use Matlab-style notation, as you mention that program. In 1D, if you have a "stencil-grid" $x=-m:m$, then you can fit an interpolating polynomial $p(x)$ of degree $2m$ to any set of function values $f(x)$. The derivatives $\partial_x^kp$ evaluated at $x=0$ are then central-difference estimates of the corresponding derivatives of $f(x)$.

For a particular set of function values this would be done via $polyfit$ and $polyder$ in Matlab. To solve for the stencil coefficients themselves, this is indeed one of the few times I have found a need to actually compute a matrix inverse. The procedure is:

1) assemble the Vandermonde matrix $A$ for the stencil $x$.

2) Invert this matrix via $C=A\backslash I$, where $I$ is the identity matrix (i.e. solve a set of $2m+1$ interpolation problems, where the columns of $I$ are your "$b$" vectors).

3) The rows of $C$ now give "filters" for the $2m+1$ coefficients of interpolating the polynomial (i.e. the coefficient vector would be $c=Cf$ for given data $f$). So, the second row of $C$ gives the coefficients of your 1st derivative stencil (divide this by $\Delta x$ for a grid with non-unit spacing).

Note that to get the stencils for the first $2m$ derivatives, you would just compute $C(k+1,:)/(k!\Delta x^{k-1})$ for the $k$th derivative. Similarly, $x=-m:m$ assumes central differences on a uniform grid, but an arbitrary $x$ can be used, e.g. to compute off-center derivatives on a non-uniform grid. And you can get moving-least-squares stencils by having a rectangular $A$ matrix (note that the Matlab "$\backslash$" command automatically solves the normal equations in this case).

Hope this helps!