1
$\begingroup$

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 :enter image description here

Matrix $b$ is $20*1$ matrix

enter image description here

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

$\endgroup$
  • $\begingroup$ That matrix is $20\times21$ it seems to me. $\endgroup$ – Kirill Dec 13 '15 at 9:32
  • $\begingroup$ @Kirill Sorry Its my mistake CD schemes so for I have seen doesn't has $0^{th}$ (mid point) so I didn't mention it. I tried to type the full matrix but latex didn't support that so I truncated it and posted as image. Please note that there is no 0 in -10 -9 -8 ... rows . $\endgroup$ – AGN Dec 13 '15 at 9:37
2
$\begingroup$

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\}. $$

$\endgroup$
  • $\begingroup$ Thanks for the answer. Sorry, it is little bit difficult for me to understand the notation used here. If you provide some external link where they solved using this style that will be helpful for me. I have some suggestions:1) Most of the CD schemes (so for I have seen) don't has mid point but this stencil has that. 2) I used $h^{-20}$ =sum(1,-20,190....184756) (Please correct me weather am I right). and cross checked the stencil by $A*x=b$ but i didn't get $b$ matrix. please check this when you are free. $\endgroup$ – AGN Dec 13 '15 at 9:18
  • $\begingroup$ @ArunGovindNeelanA Sorry, I misread the question, the coefficients should be correct now, you're asking about 1st derivative, not 20th. $\endgroup$ – Kirill Dec 13 '15 at 9:38
  • $\begingroup$ Thanks , I guess the answer is 3 or 4 decimal accuracy. Could you provide some link or refrence to calculate the co-efficients using this procedure $\endgroup$ – AGN Dec 13 '15 at 10:00
  • $\begingroup$ @ArunGovindNeelanA I'm not sure I know a reference, but I think it's standard. The coefficients are exact, what do you mean by accuracy? If you mean the residual $Ax-b$, then that can be expected as your (Vandermonde) matrix is very ill-conditioned. $\endgroup$ – Kirill Dec 13 '15 at 10:12
  • 1
    $\begingroup$ @ArunGovindNeelanA I used (x**10 * log(x)).taylor(x, 1, 20).coefficients() in sage to get the coefficients. $\endgroup$ – Kirill Dec 14 '15 at 1:29
1
$\begingroup$

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!

$\endgroup$
  • $\begingroup$ Thanks for the answer. I think step 2 need inverse calculation $C=A\I$, where I got "badly scaled" warning. Could you give me any suggestions. $\endgroup$ – AGN Dec 13 '15 at 9:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.