0

Given that $$\frac{d^2}{dx^2} = \lim_{\delta x\to 0} [f(x-\delta x) - 2f(x) + f(x+\delta x))]\cdot\frac1{\delta x^2}$$ find an appropriate matrix that could represent such a derivative operator, in a form analogous to the first derivative operator matrix.

So i know that first derivative has 0 on the central diagonal , and $1/2\delta x$ on the diagonals on either side of the central diagonal. and zero everywhere else.

Can I just multiply two first order derivatives together to get the second order? Because order matters generally in matrix multiplication i feel as though this is what is required. So I was thinking the answer would be

$1/4\delta x^2$ on the central diagonal. Does this make sense?

1 Answers1

1

This is not exactly physics... Still. The matrix corresponding to the differential operators is not unique. The matrix form depends among other on:

  • boundary conditions (for example Dirichlet $f(x=0)=0$ or Neumann $f'(x=0) =0$)
  • approximation order (please google for "difference scheme" in google).

The first derivative has $-1/(2\delta x)$ on first superdiagonal and $1/(2\delta x)$ on first subdiagonal, this is the second order difference scheme for first derivative: $$ f'(x) = \frac{f(x+h)-f(x-h)}{2h} + O(h^2)$$

For the second order you can indeed square the first order derivative matrix - this is a linear operator after all.

But the "standard" for numerics is approximating the second-derivative by "five point stencil" (https://en.wikipedia.org/wiki/Five-point_stencil) . This is often best balance between accuracy $O(h^4),$ speed, and numerical stability.

The most simple discretization of the $\frac{d^2}{dx^2}$ operator is then:

$$ \left[\begin{array}{ccccccccc} -2& 1 & 0& \ldots &0 &0 & 0\\ 1 &-2 & 1& \ldots &0 &0 & 0\\ 0 & 1 &-2& \ldots &0 &0 & 0\\ 0 & 0 & 1& \ldots &0 &0 & 0\\ \vdots & \vdots & \vdots& \ddots& \vdots& \vdots & \vdots \\ 0& 0 &0 & \ldots & 1 & 0& 0 \\ 0& 0 &0 & \ldots & -2& 1& 0 \\ 0& 0 &0 & \ldots & 1 &-2& 1 \\ 0& 0 &0 & \ldots & 0 & 1&-2 \end{array}\right], $$

But: - this is lowest possible order approximation. In practice it means that if you take a given interval, and for example study spectrum of $\frac{d^2}{dx^2}$ operator then the eigenvalues will be $\lambda_n = \lambda_\infty + \frac{c}{n^2}.$ One is typically interested in $\lambda_\infty$ rather than $\lambda_n$ and therefore one should either use large $n$ to get $\lambda_\infty$ (or try to study asymptotics of $\lambda_n$, try to deretmine $c$ and get $\lambda_\infty$ as $\lambda_n-\frac{c}{n^2}$) - this assumes open boundary conditions, that is that $u_0=0$ and $u_L=0.$ One can do for example Dirichlet boundary conditions, then $f'(0)=0$ and $u_0-u_{1}=0.$ In the latter case the proper matrix would be:

$$ \left[\begin{array}{ccccccccc} -1& 1 & 0& \ldots &0 &0 & 0\\ 1 &-2 & 1& \ldots &0 &0 & 0\\ 0 & 1 &-2& \ldots &0 &0 & 0\\ 0 & 0 & 1& \ldots &0 &0 & 0\\ \vdots & \vdots & \vdots& \ddots& \vdots& \vdots & \vdots \\ 0& 0 &0 & \ldots & 1 & 0& 0 \\ 0& 0 &0 & \ldots & -2& 1& 0 \\ 0& 0 &0 & \ldots & 1 &-2& 1 \\ 0& 0 &0 & \ldots & 0 & 1&-1 \end{array}\right], $$

(the above assumes Dirichlet boundary conditions at both ends).

Boundary conditions (with a fourth order difference scheme)

If you discretize an operator (for example $\frac{d^2}{dx^2}$ on a finite interval $[0,L]$, then the question is how to include a boundary condition. For example, while solving for a particle in a box the proper boundary condition is $f(0)=0,$ then if you set $u_i = f(i*h),$ then we have $u_0=0$. This means that the differential operator would be:

$$ \begin{eqnarray} \frac{1}{12h^2}(-u_5+16u_4-30u_3+16u_2-u_1) = f''(3h)\\ \frac{1}{12h^2}(-u_4+16u_3-30u_2+16u_1-u_0) = f''(2h)\\ \frac{1}{12h^2}(-u_3+16u_2-30u_1+16u_0-u_{-1}) = f''(h) \end{eqnarray} $$

And the problem is what is $u_{-1}.$ The solution is undefined for $x<0,$ so in principle there is no reason to set $u_{-1}$ as 0. For second order differential operator we have another boundary condition at $x=L,$ and fixing any value of $u_{-1}$ would be too restrictive.

The trick to be used here is: - we want f'(0) be well defined - we extend a function $f$ from $[0,L]$ to $[-L,L]$ by defining $f(-x)=-f(x)$ - this means $u_{-1}=-u{1}$

$$ \begin{eqnarray} \frac{1}{12h^2}(u_5+16u_4-30u_3+16u_2-u_1) = f''(3h)\\ \frac{1}{12h^2}(u_4+16u_3-30u_2+16u_1-u_0) = f''(2h)\\ \frac{1}{12h^2}(u_3+16u_2-30u_1+16u_0+u_{1}) = f''(h)\\ \frac{1}{12h^2}(u_2+16u_1-30u_0-16u_1+u_{2}) = f''(0)\\ \frac{1}{12h^2}(u_1+16u_0+30u_{-1}-16u_{-1}+u_{3}) = f''(-h) \end{eqnarray} $$

this gives straighforwardly the second derivative operator matrix as:

$$ \left( \begin{array}{cccccc} -30 & 17 & -1 & 0 & 0 & \ldots \\ 16 & -30 & 16 & -1 & 0 & \ldots \\ -1 & 16 & -30 & 16 & 1 & \ldots \\ 0 & -1 & 16 & -30 & 16 & \ldots \\ \vdots & \vdots & \vdots & \vdots & \vdots & \end{array}\right) $$

Very important the above matrix operates in $(u_1,u_2,\ldots,...)$. This means that if one computes eigenfunctions of the above matrix one will get $u_i \approx \sin(i \delta x \pi/L), i=1,2,\ldots,L-1,$ assuming close to $L,$ one chooses the same boundary conditions.

  • Isn't the simplest solution to put -2 on the diagonal and +1 on the first off diagonals? Where do all those 30 ad 16's come from? – DanielSank Feb 17 '16 at 15:42
  • It is simplest, but this is known as "second order difference scheme". If you consider numerically say an eigenproblem of finding a lowest eigenvalue, then you estimate of a lowest eigenvalue converges to the right value as $1/n^2$ with the scheme above it converges as $1/n^4$.

    So in short: to get "equally good" approximations it is often sufficient to use less grid points in the five-point stencil than in the (-1 2 -1) scheme.

     –  Feb 17 '16 at 16:14
  • Oh I understand now. Might be good to explain the simplest case first. Optimized solutions are often un-necessary, and even when they are, they're usually easier to understand after understanding the simplest non-optimized case. – DanielSank Feb 17 '16 at 16:48