Derivation of numerical scheme for linear transport equation on a variable

This topic is a copy of that one on Mathematica.StackExchange. The question is about automatical derivation of coefficients of numerical scheme on a variable stencil. So, lets consider linear 1d transport equation

(1)\qquad u_t+u_x=0.

To discretize (1), let’s consider the following stencil:

(2)\qquad (t^{n+1},x_{j-1}),\,(t^{n},x_{j-1}),\,(t^{n},x_{j}).

Then, let’s consider the following approximation of (1) on (2):

(3)\qquad u_j^{n+1} = a_1u_{j-1}^n+a_2u_j^n,

where a_1,\,a_2 to be determined by Taylor series expansion. My first step seems to be as follows:

eqn = D[u[x, t], t] + D[u[x, t], x] == 0
approx = Series[
u[x, t + \[Tau]], {t, 0, 2}] - (a1 Series[u[x - h, t], {x, 0, 2}] + a2 Series[u[x, t], {x, 0, 2}])

The question is how to derive differential consequence from (1) to obtain system of two equation on a_i,\,i=1,2? Means, that we hold stencil (2) and approximation (3), consequently. The question is how to obtain a_1,a_2 by unknow coefficient method with the usage of differential

Here copy of the accepted answer given on the linked above MSE thread:


Here’s how I go about doing this type of algebra with Mathematica:

eqn = D[u[x, t], t] + D[u[x, t], x] == 0;

First, expand the function u[x, t + dt] to first order in dt around 0 and the function u[x -h, t] to first order in h around 0:

eqn1 = u[x, t + dt] == Series[u[x, t + dt], {dt, 0, 1}] // Normal
eqn2 = u[x - h, t] == Series[u[x - h, t], {h, 0, 1}] // Normal

enter image description here

We want to eliminate the partial derivatives, so we solve for them:

sols = First@Solve[{eqn1, eqn2}, {D[u[x, t], t], D[u[x, t], x]}]

enter image description here

Finally, we use the fact that the partial derivatives are negative each other because of the original differential equation to eliminate them from the equation:

differenceEquation = eqn /. sols

enter image description here

Finally, we solve for the next time step and collect terms on the right hand side

Collect[Equal @@ Solve[differenceEquation, u[x, dt + t]][[1, 1]], u[x, t]]

enter image description here

From there, you can read off a1 and a2.