4
$\begingroup$

The question is about automatica derivation of coefficients of numerical scheme on a variable stencil.

So, lets consider 1d transport equation \begin{equation} (1)\qquad u_t+u_x=0. \end{equation} To discretize (1), let's consider the following stencil: \begin{equation} (2)\qquad (t^{n+1},x_{j-1}),\,(t^{n},x_{j-1}),\,(t^{n},x_{j}). \end{equation} 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 + τ], {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$?

$\endgroup$
  • $\begingroup$ So what is the question? And do you want to do this with the computing software Mathematica? If so, please explain what specifically you are having difficulties with in implementing your algorithm in Mathematica. If not, your question belongs elsewhere, perhaps math.se. $\endgroup$ – march Nov 3 '15 at 19:29
  • $\begingroup$ Sure, edited just now. $\endgroup$ – Oleg Kravchenko Nov 3 '15 at 19:37
  • $\begingroup$ I'm still a little confused. Do you want to use Mathematica to show that the discretization of Eq. 1 results in Eq. 3, and as a consequence derive what $a_1$ and $a_2$ should be? $\endgroup$ – march Nov 3 '15 at 19:40
  • $\begingroup$ I guess, my question is the second one. 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 consequences of (1). $\endgroup$ – Oleg Kravchenko Nov 3 '15 at 19:44
6
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Ok, but how such derivation could be adopted on the following scheme, for example, $$ u_j^{n+1} = a_1u_{j-1}^n + a_2 u_{j}^n + a_3 u_{j+1}^n, \quad \text{(Lax-Wendroff scheme)} $$ or $$ u_j^{n+1} = a_1u_{j-1}^n + a_2 u_{j}^n + a_3 u_{j+1}^n + a_4 u_{j+2}^n? $$ $\endgroup$ – Oleg Kravchenko Nov 23 '15 at 18:26
  • 1
    $\begingroup$ @olekravchenko. I can't think of a way of automating derivations of the different schemes, because they use different choices for the finite-difference approximations of the derivatives. I would just generalize what I did above, but in general my method requires knowledge of how the derivation is done in the first place. Nonetheless, in the meantime, I'll think on it. You might consider asking a new question about how to automate the derivation of general schemes on general stencils. I just answered this question originally assuming that all you wanted was the particular one shown. $\endgroup$ – march Nov 23 '15 at 19:00
  • $\begingroup$ In fact I wrote on a variable stencil to ask question in more wide sence. Anyway, let me follow your advice and ask new question about automatical derivation on quite general stencil. $\endgroup$ – Oleg Kravchenko Nov 23 '15 at 19:39

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.