# How to derive finite-difference scheme automatically on a quite general stencil

Info

This question is a generalization of the following one Derivation of numerical scheme for linear transport equation on a variable stencil.

Statement of a problem

Linear scalar hyperbolic (transport/advection) equation is considered $$u_t + u_x = 0, \tag{0a}$$ here $u = u(x,t)$. If it is necessary, one could also add into a consideration the initial condotion as well $$\left.u(x,t)\right|_{t=0}=u_0(x). \tag{0b}$$

So, let's consider quite general stencil, such as: $$(t^{n+1}, x_j),\,(t^{n}, x_{j-1}),\,(t^{n}, x_{j}),\,(t^{n}, x_{j+1}),\,(t^{n}, x_{j+2}), \tag{1}$$ so, 4 points on $n$-th time step and 1 point on $n+1$ time step. The explicit approximation of $u_j^{n+1}$ to be as follows: $$u_j^{n+1} = a_1u_{j-1}^n + a_2u_{j}^n + a_3u_{j+1}^n + a_4u_{j+2}^n. \tag{2}$$

Notice

Stencil (1) could be different, for example, one more point $(t^n,\,x_{j+3})$ could be added.

Question

The question is how to obtain undetermined coefficients $a_i,~i=\overline{1,4}$ automatically wether is it possible or not?

Possible ways of solution ( ... to be updated)

In general notation Eq. (2) to be read as $$u_j^{n+1} = \sum_{i=-k_L}^{k_L}a_i u_{j+i}^n,\quad i=-k_L,-k_L+1,\ldots,0,\ldots,k_R-1,k_R. \tag{3}$$ For example, if $k_L=1,k_R=2$, then obtain particular case $$u_j^{n+1} = \sum_{i=-1}^2a_iu_{j+1}^n=a_{-1}u_{j-1}^n+a_0u_j^n+a_1u_{j+1}^n+a_2u_{j+1}^n. \tag{4}$$

Obviously, that there are several question to be solved:

1. How to set stencil (1) in WM.

2. How to set the note and how to expand all terms in right hand side of (3) around the choosen node.

3. How to use a differential consequences (DC) of (0a), i.e. \begin{align} u_t &= -u_x, \\ u_{tx} &= -u_{xx}, \\ u_{tt} &= -u_{xt} = u_{xx}, \\ u_{txx} &= - u_{xxx}, \\ u_{ttt} &= -u_{xtt} = u_{xxt} = -u_{xxx}, \\ &\ldots \tag{5} \end{align}

DC (5) are applied in expansion of left hand side of (3) around $(x_j,t^n)$ node as follows: \begin{align} & u_j^{n+1} = u_j^n + \tau\bigl(u_t\bigr)_j^n + \frac{\tau^2}{2!}\bigl(u_{tt}\bigr)_j^n + \frac{\tau^3}{3!}\bigl(u_{ttt}\bigr)_j^n + O(\tau^4) = \\ =\,& u_j^n - \tau\bigl(u_x\bigr)_j^n + \frac{\tau^2}{2!}\bigl(u_{xx}\bigr)_j^n - \frac{\tau^3}{3!}\bigl(u_{xxx}\bigr)_j^n + O(h^4). \tag{6} \end{align}

Semi-automatic solution

WM code:

    {
eq = D[u[x, t], t] + D[u[x, t], x] == 0,
eqdc01 = D[eq, t],
eqdc02 = D[eq, t, t],
eqdc03 = D[eq, x],
sol00 = Solve[eq, D[u[x, t], t]] // First,
sol01 = Solve[eqdc01, D[u[x, t], {x, 0}, {t, 2}]] // First,
sol02 = Solve[eqdc03, D[u[x, t], {x, 1}, {t, 1}]] // First,
sol03 = Solve[eqdc02, D[u[x, t], {x, 1}, {t, 2}]] // First,
sol04 = Solve[eqdc02, D[u[x, t], {x, 0}, {t, 3}]] // First,
sol05 = D[sol02, {t, 1}],
sol06 = D[sol00, {x, 2}],
sol07 = sol05 /. sol06,
sol08 = sol04 /. sol07
} // Column


    {
se01 = Series[u[x, t + τ], {τ, 0, 3}] // Normal,
lhs = se01 /. sol00 /. sol01 /. sol02 /. sol08 // Normal,
rhs = Sum[
Subscript[a,
i] (Series[u[x + i h, t], {h, 0, 3}] // Normal), {i, -1, 2, 1}]
} // Column


   {
eqnSys = lhs - rhs;
eqnSysGathered =
Collect[eqnSys, Table[D[u[x, t], {x, k}], {k, 0, 3}]]
} // Column
eqnSys2 =
Table[Coefficient[eqnSysGathered /. τ -> σ h,
D[u[x, t], {x, k}]] == 0, {k, 0, 3}];
eqnSys2Sol =
Solve[eqnSys2, Table[Subscript[a, i], {i, -1, 2, 1}]] // First;
% // MatrixForm


   Subscript[u, j]^(n + 1) ==
Sum[Subscript[a, i] Subscript[u, j + i]^n, {i, -1, 2, 1}] /. eqnSys2Sol


Tracking

UPD1. 28/11/15: Section "Possible ways of solution" was added.

UPD2. 28/11/15: Section "Semi-automatic solution" was added.

• The earlier comment by @marsh is correct. Basically, any set of as that sum to zero is a valid four-point finite difference approximation to D[u[x, t], x], although some are better than others. To derive the coefficients automatically requires three additional constraints. – bbgodfrey Nov 23 '15 at 20:30
• From memory these are derived (somewhere) in the extensive NDSolve documentation. – Mike Honeychurch Nov 23 '15 at 22:20
• @MikeHoneychurch You probably are referring to Spatial Derivative Approximations. – bbgodfrey Nov 23 '15 at 23:02
• @bbgodfrey yes I was. thanks. I thought they were derived actually but in any case the tables are useful and I don't think the derivations are difficult in any case (used to do this stuff but those brain cells have long since been archived). – Mike Honeychurch Nov 23 '15 at 23:36
• @MikeHoneychurch The derivations of a few of the many discretizations do appear just above the tables. Fornberg has an extensive article on this topic. – bbgodfrey Nov 23 '15 at 23:55