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?
This is not exactly physics... Still. The matrix corresponding to the differential operators is not unique. The matrix form depends among other on:
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.