Moving Node Method for Differential Equations

Abstract

The chapter contains information about new approaches to solving boundary value problems for differential equations. It introduces a new method of moving nodes. Based on the approximation of differential equations (by the finite difference method or the control volume method), introducing the concept of a moving node, approximately analytical solutions are obtained. To increase the accuracy of the obtained analytical solutions, multipoint moving nodes are used. The moving node method is used to construct compact circuits. The moving node method allows you to investigate the diskette equation for monotonicity, as well as the approximation error of the differential equation. Various test problems are considered.

1. Introduction

Methods for solving problems of mathematical physics can be divided into the following four classes [1, 2, 3, 4, 5, 6, 7].

Analytical methods (the method of separation of variables, the method of characteristics, the method of Green’s functions [8], etc.) have a relatively low degree of universality, i.e. focused on solving rather narrow classes of problems. As a result of applying these methods, a solution is obtained in the form of analytical formulas. The use of these formulas for the implementation of the calculation may require the solution of auxiliary computational problems (solution of nonlinear equations, calculation of special functions, numerical integration, summation of an infinite series). Nevertheless, in a number of cases, the application of these methods makes it possible to quickly and with high accuracy calculate the desired solution.

Approximate analytical methods (projection, variational methods, small parameter methods, operational methods, and various iterative methods [4, 9]) are more universal than analytical ones. The use of such methods involves modifying the original problem or changing the problem statement in such a way that the new problem can be solved by the analytical method, and its solution itself differs little enough from the solution of the original problem.

Numerical methods (finite difference method, method of lines, control volume method, finite element method, etc. [1, 2, 5, 6, 7, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34]) are very universal methods. Often used to solve nonlinear problems of mathematical physics, as well as linear problems with variable operator coefficients.

Probabilistic methods (Monte Carlo methods) are highly versatile. It can be used to calculate discontinuous solutions. However, they require large amounts of calculations and, as a rule, lose with the computational complexity of the above methods when solving such problems to which these methods are applicable.

Comparing methods for solving problems of mathematical physics, it is impossible to give unconditional superiority to any of them. Any of them may be the best for solving problems of a certain class. At the same time, when characterizing a specific method, it is advisable to highlight those features that often determine its advantages or disadvantages in practical application compared to an alternative method.

The advantages of the finite difference method include its high universality, for example, much higher than that of analytical methods. The application of this method is often characterized by the relative simplicity of constructing a decision algorithm and its software implementation. Often it is possible to parallelize the decision algorithm.

The shortcomings of the method include: the problematic nature of its use on irregular grids; a very rapid increase in computational complexity with an increase in the dimension of the problem (an increase in the number of unknown variables); the complexity of the analytical study of the properties of the difference scheme.

The proposed method of moving nodes combines numerical and analytical methods [7, 8, 13, 35, 36, 37, 38]. In this case, we can obtain, on the one hand, an approximate analytical solution to the problem, which is not related to the methods listed above. On the other hand, this method allows one to obtain compact discrete approximations of the original problem. Note that obtaining an approximate analytical solution to differential equations is based on numerical methods. The nature of numerical methods also makes it possible to obtain an approximate analytical expression for solving differential equations. For this, a so-called “movable node” is introduced.

The aim of the study is to develop a computing technology based on the proposed method of moving nodes, develop a two-point convective-diffusion problem an analytical method generated by numerical methods based on the method of moving nodes, and give test examples.

2. Derivation of approximate analytical solutions of differential equations by the moving nodes method

This chapter introduces the concept of a roaming node and provides approximate solutions to simple problems using a moving node. We also studied the derivation algorithm for nonstationary and two-dimensional problems.

Note that the concept of a movable node in this context is considered for the first time.

2.1 The concept of a moving node

The solution of differential equations (DE) (ordinary or partial derivatives) by the method of finite differences is based on a finite-difference approximation of derivatives. When applying the finite difference method to the solution of DE, there is a transition from a continuous region to a finite difference one. A grid of “nodal points” is introduced into the solution area. Representing the derivatives in a finite difference form, they bring it to the form of a difference equation. The difference equation is written for all grid nodes and results in a system of algebraic equations [4, 36].

Most of the DEs found in the equations of mathematical physics contain only partial derivatives of the first and second orders, while for the approximation of the derivatives they try to use no more than three nodes of the difference grid (in the case of ordinary DEs) (Figure 1). Let the node W and E be considered fixed, and the node x changes into segments (W, E). Then the approximation of the derivatives (first or second order) also changes based on the location of the node. The node x is said to be movable.

Figure 1.

One moving node.

You can increase the number of moved nodes. Let us select additional moving nodes as follows: x₁ = (W + x)/2, x₂ = (W + x)/2. When node x changes its position, x₁ and x₂ automatically change their positions (Figure 2). In this way, you can increase the number of moved nodes. The increase in the number of moved nodes is related to the accuracy of the difference equations.

Figure 2.

Three moving node.

The displacement of nodal points is not only related to finite-difference equations, this approach can be successfully applied when discretizing differential equations using the control volume method.

2.2 Obtaining an approximate analytical solution with one moving node

Let, it is necessary to find $Ф\left(х\right)$ a solution to the DE in the region $W\le x\le E$ with the corresponding boundary conditions. Let us take an arbitrary point $x\in \left(W,E\right).$ We have three nodes: $W,E$ boundary nodes and an internal node $x$. The position of a point inside the region is determined by the node being moved $x$. The difference equation is usually written for an arbitrary node, $x$. When approximating differential operators, the first derivatives on the moving node are approximated by different relations:

$$\frac{\mathit{dФ}\left(х\right)}{\mathit{dx}}\approx \frac{U\left(x\right)-U\left(W\right)}{x-W},$$E1

$$\frac{\mathit{dФ}\left(х\right)}{\mathit{dx}}\approx \frac{U\left(E\right)-U\left(x\right)}{E-x},$$E2

$$\frac{\mathit{dФ}\left(х\right)}{\mathit{dx}}\approx \frac{U\left(E\right)-U\left(W\right)}{E-W}.$$E3

The approximation of the derivative by (1) and (2) is called the approximation of this derivative using a one-sided difference, and (3) is the approximation using the central difference.

The second derivative on the moving node is approximated as follows [4] (similarly to the approximation of the second derivative in a non-uniform grid):

$$\frac{d^2Ф\left(х\right)}{{\mathit{dx}}^2}\approx \frac{2}{E-W}\left(\frac{U\left(E\right)-U\left(x\right)}{E-x}-\frac{U\left(x\right)-U\left(W\right)}{x-W}\right)$$E4

Let us consider some model problems of applying the moving nodes method (MNM) to obtain an analytical solution.

2.2.1 Flow in a flat pipe

The flow of a viscous fluid in a flat pipe in a one-dimensional formulation is described by the equation

$$\frac{d^{2}U}{{\mathit{dy}}^2}=-\frac{\mathrm{\Delta }p}{\mu l}$$E5

where $U$ is the fluid velocity, is the vertical coordinate perpendicular to the flow, $\mathrm{\Delta }p/l$ is the pressure drop (const), $\mu$ is the viscosity. Let $y=0$ and $y=h$ motionless walls.

We average (5) over the liquid volume: $\left[y/2,\left(h-y\right)/2\right]$, here “y” is a moving node (Figure 3). Then we have

$$\underset{y/2}{\overset{\left(h+y\right)/2}{\int }}\frac{d^{2}U}{{\mathit{dy}}^2}\mathit{dy}=\underset{y/2}{\overset{\left(h+y\right)/2}{\int }}\left(-\frac{\mathrm{\Delta }p}{\mu l}\right)\mathit{dy}$$

Figure 3

Control volume.

From here

$${\left.\frac{\mathit{dU}}{\mathit{dy}}\right|}_{\left(h+y\right)/2}-{\left.\frac{\mathit{dU}}{\mathit{dy}}\right|}_{y/2}=\left(-\frac{\mathrm{\Delta }p}{\mu l}\right)\frac{h}{2}$$E6

We replace the derivatives in (6) with the difference relation:

$${\left.\frac{\mathit{dU}}{\mathit{dy}}\right|}_{\left(h+y\right)/2}\approx \frac{u\left(h\right)-u\left(y\right)}{h-y},{\left.\frac{\mathit{dU}}{\mathit{dy}}\right|}_{y/2}\approx \frac{u\left(y\right)-u\left(0\right)}{y}.$$E7

Here $u\left(y\right)$ is an approximate value $U\left(y\right)$. Thus, approximation (5) with respect to the moving node has the form:

$$\frac{u\left(h\right)-u\left(y\right)}{h-y}-\frac{u\left(y\right)-u\left(0\right)}{y}=\left(-\frac{\mathrm{\Delta }p}{\mu l}\right)\frac{h}{2}.$$E8

Hence, taking into account the no-slip condition (u(h) = u(0) = 0)

$$u\left(y\right)=-\frac{\mathrm{\Delta }p}{2\mu l}y\left(h-y\right).$$

Here $u\left(y\right)$ is the average solution. For this problem, the averaged solution coincides with the exact solution.

This means that the approximation (7) for this problem is exact. The reason for the coincidence of the solution obtained with the help of the MNM with one node and the exact solution is explained by the following fact.

Lagrange’s mean value theorem states that if a function $f\left(x\right)$ is continuous on an interval $\left[a,b\right]$ and differentiable on an interval $\left(a,b\right)$, then in this interval there is at least one $x=\xi$ point such that

$$\frac{f\left(b\right)-f\left(a\right)}{b-a}=f^{\prime }\left(\xi \right).$$E9

It is easy to check that if $f\left(x\right)$ represents a parabola, then in (9) $\zeta =\left(a+b\right)/2.$ The exact solution (5) is a parabola. Integrating (5) over the control volume $\left[x/2,\left(h+x\right)/2\right]$, we obtain

$$\underset{y/2}{\overset{\left(h+y\right)/2}{\int }}\frac{d^{2}u}{{\mathit{dy}}^2}\mathit{dy}={\left.\frac{\mathit{du}}{\mathit{dy}}\right|}_{\left(h+y\right)/2}-{\left.\frac{\mathit{du}}{\mathit{dy}}\right|}_{y/2}=\underset{y/2}{\overset{\left(h+y\right)/2}{\int }}\frac{1}{\mu }\frac{\mathrm{\Delta }p}{l}\mathit{dy}.$$

Since $u\left(y\right)$ there is a parabola, therefore

$${\left.\frac{\mathit{du}}{\mathit{dy}}\right|}_{\left(h+y\right)/2}=\frac{u\left(h\right)-u\left(y\right)}{h-y},{\left.\frac{\mathit{du}}{\mathit{dy}}\right|}_{y/2}=\frac{u\left(y\right)-u\left(0\right)}{y-0},$$

and (8) is the exact difference analog of (5).

2.2.2 Heat distribution in the plate

Heat propagation in the plate is described by the equation

$$\frac{d^{2}T}{{\mathit{dx}}^2}+\frac{q}{k}=0,\frac{\mathit{dT}\left(0\right)}{\mathit{dx}}=0,T\left(1\right)=1$$E10

where k is the thermal conductivity and q is the heat release per unit volume (k and q = const). It is assumed that the source does not depend on temperature. Replacing (10) with a difference equation with a moving node, we have

$$\frac{2}{1-0}\left[\frac{T\left(1\right)-T\left(x\right)}{1-x}-\frac{T\left(x\right)-T\left(0\right)}{x-0}\right]+\frac{q}{k}=0$$E11

Solving Eq. (11), we obtain

$$T\left(x\right)=1+\frac{q}{2k}\left(1-x^2\right)$$E12

Solution (12) coincides with the exact solution. Note that the exact solution is obtained not only for the Dirichlet problem but as for the problem of flow in a flat pipe. Here the boundary conditions are of mixed type.

2.2.3 Magnetohydrodynamic Couette flow

Consider the Couette flow, when a conducting fluid flows in a uniform magnetic field between two plates, one of which is stationary, and the other moves in its own plane at a constant speed. Based on the Navier-Stokes equation, taking into account the magnetic field and taking into account the one-dimensionality of the flow, it can be written in a dimensionless form as follows:

$$\frac{d^{2}u}{{\mathit{dy}}^2}-M^{2}u=P$$E13

Boundary conditions

$$u\left(0\right)=0,u\left(1\right)=1$$E14

Here, $u$ is the dimensionless flow velocity and $y$ is the dimensionless coordinate. Dimensionless quantities $M$ – Hartmann number, $P$ – pressure coefficient (М and Р = const).

Replacing the second-order derivative in (13) with a difference relation similar to (7), and considering the boundary condition (14), we can obtain an approximate solution

$$u_1\left(y\right)=\frac{2y-\mathit{Py}\left(1-y\right)}{2+M^{2}y\left(1-y\right)}$$E15

This solution comes close to the exact solution (Figure 4).

Figure 4.

Comparison of exact and approximate solutions (M = 2, P = 4). The solid line is the exact solution, the dotted line is according to (5).

2.2.4 The method of moving nodes for the convection-diffusion equation

Consider the transport equation

$$\frac{\mathit{dФ}}{\mathit{dx}}=\frac{1}{\mathit{Pe}}\frac{d^2Ф}{{\mathit{dx}}^2}+S\left(x\right),$$E16

Here, $Ф$ the unknown function, $S\left(x\right)$ the source, $\mathit{Pe}$ is the Peclet number. The equation is considered under appropriate boundary conditions.

The convective term of Eq. (16) is approximated by (1), and the diffusion term by (4). Consider (16) into segments with boundary conditions $Ф\left(0\right)=0,Ф\left(1\right)=1$ and $S\left(x\right)=0$. Then, using the upwind scheme, we replace Eq. (16) with a difference equation that looks like this:

$$\frac{U\left(x\right)}{x}=\frac{2}{\mathit{Pe}}\left(\frac{1-U\left(x\right)}{1-x}-\frac{U\left(x\right)}{x}\right)$$E17

From here, we can easily determine $U\left(x\right)$:

$$U\left(x\right)=\frac{2x}{2+\mathit{Pe}\left(1-x\right)}$$E18

Figure 5 shows a comparison of the exact and approximate solutions. The solid line corresponds to the exact solution, and the dotted line corresponds to the solution (8). It can be seen from the graph that numerical diffusion takes place.

Figure 5.

Comparison of exact and approximate solutions. Pe = 10.

For $Ф\left(0\right)=0,Ф\left(1\right)=1$ and $S\left(x\right)=5cos4x$ $\mathit{Pe}=5$, the results of the exact and approximate solutions are shown in Figure 6. It can be seen from the graph that there are large errors. Here the Peclet number plays an important role. Indeed, for $Ф\left(0\right)=0,Ф\left(1\right)=1$ and $S\left(x\right)=5cos4\mathit{x Pe}=0,1$, we obtain solutions shown in Figure 7, which shows that the approximate and exact solutions are close.

Figure 6.

Comparison of exact and approximate solutions. Pe = 5.

Figure 7.

Comparison of exact and approximate solutions. Pe = 0.1.

2.2.5 Equation with variable coefficient

Consider the equation

$$\epsilon u^{″}\left(x\right)+2x{u}^{\prime }\left(x\right)=0,$$E19

into segments (−1,1) with boundary $u\left(-1\right)=-1,u\left(1\right)=2$ conditions $u\left(-1\right)=-1,u\left(1\right)=2.$ The exact solution is determined through the error functions:

$$u\left(x\right)=\frac{\mathit{erf}\left(1/\sqrt{\epsilon }\right)+3\mathit{erf}\left(x/\sqrt{\epsilon }\right)}{2\mathit{erf}\left(1/\sqrt{\epsilon }\right)}.$$

The difference scheme with a moving node for (19) has the form (upwind scheme):

$$\epsilon \left[\frac{2-U\left(x\right)}{1-x}-\frac{U\left(x\right)+1}{x+1}\right]+\frac{1}{2}\left(2x-\left|2x\right|\right)\frac{U\left(x\right)+1}{1+x}+\frac{1}{2}\left(2x+\left|2x\right|\right)\frac{2-U\left(x\right)}{1+x}=0.$$

Solving this equation with respect to $U\left(x\right)$, we obtain an approximate analytical solution. Figure 8 compares the solutions of the exact and approximate analytical solution; the solid line corresponds to the exact solution, and the dotted line corresponds to the approximate one.

Figure 8.

Solution comparisons.

Remark 1. In the given examples, the convective term is approximated by the upwind scheme. Other approximations can be used to improve.

Remark 2. In the above examples, the approximation of the term with the source is carried out constant in the considered moving segment. For improvement, other approximations can be used to obtain an improved solution.

2.3 Obtaining an analytical solution with several moving nodes

2.3.1 Moving nodes method for a one-dimensional convective-diffusion problem

Due to the importance of convective-diffusion problems, we will apply multipoint MNM to such problems [14]:

$$\frac{\mathit{dФ}}{\mathit{dx}}=\frac{1}{\mathit{Pe}}\frac{d^2Ф}{{\mathit{dx}}^2}+S\left(x\right).$$E20

Let us take an arbitrary one node inside the segment $x\in \left(W,E\right)$.

Let us consider a difference analog of Eq. (20), in which the convective term is approximated by a one-sided difference relation.

Then the upwind scheme has the form:

$$\mathit{Pe}\frac{U^1-U_W^1}{x-W}=\frac{2}{\left(E-W\right)}\left(\frac{U_E^1-U^1}{E-x}-\frac{U^1-U_W^1}{x-W}\right)+\mathit{Pe}\cdot S\left(x\right).$$E21

This schema can be rewritten like this:

$$a_P^1{U}^1=a_E^1{U}_E^1+a_W^1{U}_W^1+F^1\left(x\right),$$E22

Here

$$\begin{array}{l}{a}_E^1=\frac{2}{\left(E-W\right)\left(E-x\right)},a_W^1=\frac{\mathit{Pe}}{\left(x-W\right)}+\frac{2}{\left(E-W\right)\left(x-W\right)},a_P^1=a_E^1+a_W^1,\\ F^1\left(x\right)=\mathit{Pe}\cdot S\left(x\right)\end{array}$$

Hence, we have

$$U^1=\frac{2\left(x-W\right)U_E^1+\left(E-x\right)\left(2+\mathit{Pe}\left(E-W\right)\right)U_W^1}{\left(E-W\right)\left(2+\mathit{Pe}\left(E-x\right)\right)}+\frac{\left(x-W\right)\left(E-x\right)}{2+\mathit{Pe}\left(E-x\right)}\mathit{Pe}\cdot S\left(x\right)$$E23

When $x\in \left(W,E\right)$ changes its position (let us make it moveable within the interval $\left(W,E\right)$), based on (23) we obtain the values of the unknown function in each position. In other words, $U^1$ obtained with the help of (23), will give us an approximate solution to the problem. Note that in this case, $U_W^1=Ф\left(W\right),U_E^1=Ф\left(E\right).$ The superscript corresponds to the number of nodes being moved.

Adding additional moving nodes $x_1=\frac{x+W}{2},x_2=\frac{x+E}{2}.$.

Now we have three moving nodes $x,x_1,x_2.$ Note that if x changes its position, then $x_1$ and $x_2$ also changes its position.

A scheme of type (21) for a segment $\left[W,x\right]$ has the form:

$$\mathit{Pe}\frac{U_1^3-U_W^3}{\left(x-W\right)/2}=\frac{2}{\left(x-W\right)}\left(\frac{U^3-U_1^3}{x-x_1}-\frac{U_1^3-U_W^3}{x_1-W}\right)+\mathit{Pe}\cdot S\left(x_1\right).$$E24

Here $U_1^3=U^3\left(x_1\right)$.

A scheme of type (21) for a segment $\left[x,E\right]$ has the form

$$\mathit{Pe}\frac{U_2^3-U^3}{\left(E-x\right)/2}=\frac{2}{\left(E-x\right)}\left(\frac{U_E^3-U_2^3}{E-x_2}-\frac{U_2^3-U^3}{x_2-x}\right)+\mathit{Pe}\cdot S\left(x_2\right).$$E25

Scheme upstream for a segment $\left[x_1,x_2\right]:$

$$\mathit{Pe}\frac{U^3-U_1^3}{x-x_1}=\frac{2}{\left(x_2-x_1\right)}\left(\frac{U_2^3-U^3}{x_2-x}-\frac{U^3-U_1^3}{x-x_1}\right)+\mathit{Pe}\cdot S\left(x\right).$$E26

Here $U_2^3=U^3\left(x_2\right)$.

In (26) we exclude $U_1^3,U_2^3$ using (24) and (25). Then we get the following diagram:

$$\mathit{Pe}\frac{U^3-U_W^3}{\frac{\left(x-W\right)}{2}\cdot \left(1+{\tau }_1\right)}=\frac{4}{\left(E-W\right)}\left(\frac{U_E^3-U^3}{\frac{E-x}{2}\cdot \left(1+{\gamma }_1\right)}-\frac{U^3-U_W^3}{\frac{x-W}{2}\cdot \left(1+{\tau }_1\right)}\right)+F^3\left(x\right)$$E27

Here we have introduced the notation.

$${\tau }_1=2/\left(2+\sigma \right),{\gamma }_1=\left(2+\theta \right)/2,\sigma =\mathit{Pe}\left(x-W\right),\theta =\mathit{Pe}\left(E-x\right),$$

$$F^3\left(x\right)=\mathit{Pe}\cdot S\left(x\right)+\frac{4+\mathit{Pe}\cdot \left(E-W\right)}{E-W}\cdot \frac{1-{\tau }_1}{1+{\tau }_1}\cdot S\left(x_1\right)+\frac{4}{E-W}\cdot \frac{{\gamma }_1-1}{{\gamma }_1+1}\cdot S\left(x_2\right).$$

And $U_W^3=Ф\left(W\right),U_E^3=Ф\left(E\right).$

(25) can be rewritten as follows:

$$a_P^3{U}^3=a_E^3{U}_E^3+a_W^3{U}_W^3+F^3\left(x\right),$$E28

where

$a_E^3=\frac{8}{\left(E-W\right)\left(E-x\right)\left(1+{\gamma }_1\right)},a_W^3=\frac{2\mathit{Pe}}{\left(x-W\right)\left(1+{\tau }_1\right)}+\frac{8}{\left(E-W\right)\left(x-W\right)\left(1+{\tau }_1\right)},$$a_P^3=a_W^3+a_E^3.$

Increase the number of moved nodes:

$x_1^{-}=\frac{x_1+W}{2}=\frac{x+3W}{4},$$x_1^{+}=\frac{x_1+x}{2}=\frac{3x+W}{4},$

$x_2^{-}=\frac{x_2+x}{2}=\frac{3x+E}{4},$$x_2^{+}=\frac{x_2+E}{2}=\frac{x+3E}{4}.$

In the difference scheme (28), the unknown function appears at three nodes: W, x, E. The function S is calculated at points $x_1,x,x_2.$ Let us write a scheme of type (28) for each of the segments $\left[W,x\right]$ and $\left[x_1,x_2\right].$

The scheme of type (28) for a segment has the form:

$$a_{x_1}^3{U}_{x_1}^3=a_x^3{U}_x^3+a_{W^{-}}^3{U}_W^3+F_{-}^3\left(x_1\right),$$E29

where

$a_x^3=\frac{8}{\left(x-W\right)\left(x-x_1\right)\left(1+{\gamma }_1^{-}\right)},a_{W-}^3=\frac{2\mathit{Pe}}{\left(x_1-W\right)\left(1+{\tau }_1^{-}\right)}+\frac{8}{\left(x-W\right)\left(x_1-W\right)\left(1+{\tau }_1^{-}\right)},$

$a_{x_1}^3=a_x^3+a_{W^{-}}^3$, $F_{-}^3\left(x_1\right)=\mathit{Pe}\cdot S\left(x_1\right)+\frac{4+\mathit{Pe}\cdot \left(x-W\right)}{x-W}\cdot \frac{1-{\tau }_1^{-}}{1+{\tau }_1^{-}}\cdot S\left(x_1^{-}\right)+\frac{4}{x-W}\cdot \frac{{\gamma }_1^{-}-1}{{\gamma }_1^{-}+1}\cdot S\left(x_1^{+}\right)$,

${\tau }_1^{-}=2/\left(2+{\sigma }^{-}\right),{\gamma }_1^{-}=\left(2+{\theta }^{-}\right)/2$, ${\sigma }^{-}=\mathit{Pe}\left(x_1-W\right),{\theta }^{-}=\mathit{Pe}\left(x-x_1\right).$

Similarly, we write a scheme of type (29) for the segments $\left[x,W\right]$ and $\left[x_1,x_2\right].$ Excluding the obtained three systems of equations $U_{x_1}^3$ and $U_{x_2}^3$ obtain a scheme with seven movable nodes:

where

$\begin{array}{l}{a}_E^7=\frac{2^5\left(1-{\gamma }_2\right)}{\left(E-W\right)\left(E-x\right)\left(1-{\gamma }_2^4\right)},a_W^7=\frac{4\mathit{Pe}\left(1-{\tau }_2\right)}{\left(x-W\right)\left(1-{\tau }_2^4\right)}+\frac{2^5\left(1-{\tau }_2\right)}{\left(E-W\right)\left(x-W\right)\left(1-{\tau }_2^4\right)},a_P^7=a_W^7+a_E^7.\\ {\tau }_2=4/\left(4+\sigma \right),{\gamma }_2=\left(4+\theta \right)/4\\ F^7\left(x\right)=\mathit{Pe}\cdot S\left(x\right)+\frac{8+\mathit{Pe}\cdot \left(x-W\right)}{x-W}\cdot \frac{{\left(1-{\tau }_2\right)}^2}{1-{\tau }_2^4}\cdot \sum _{j=1}^3\sum _{i=1}^j{\tau }_2^{i-1}S\left(W+j\frac{x-W}{4}\right)-\\ \frac{8}{E-W}\cdot \frac{{\left(1-{\gamma }_2\right)}^2}{1-{\gamma }_2^4}\cdot \sum _{j=1}^3\sum _{i=1}^j{\gamma }_2^{i-1}S\left(x+\left(4-j\right)\frac{E-x}{4}\right).\end{array}$

Continuing in this way, we can get a scheme with $2^k-1$ moving nodes

$$a_P^{\left(2^k-1\right)}{U}^{\left(2^k-1\right)}=a_E^{\left(2^k-1\right)}{U}_E^{\left(2^k-1\right)}+a_W^{\left(2^k-1\right)}{U}_W^{\left(2^k-1\right)}+F^{\left(2^k-1\right)}\left(x\right),$$E30

where $a_E^{\left(2^k-1\right)}=\frac{2^{2k+1}\left(1-{\gamma }_k\right)}{\left(E-W\right)\left(E-x\right)\left(1-{\gamma }_k^{2^k}\right)},a_W^{\left(2^k-1\right)}=\frac{2^{2k+1}\mathit{Pe}\left(1-{\tau }_k\right)}{\left(x-W\right)\left(1-{\tau }_k^{2^k}\right)}+\frac{2^{2k+1}\left(1-{\tau }_k\right)}{\left(E-W\right)\left(x-W\right)\left(1-{\tau }_k^{2^k}\right)},$ $a_P^{\left(2^{л}-1\right)}=a_W^{\left(2^{л}-1\right)}+a_E^{\left(2^{л}-1\right)}.{\tau }_k=2^k/\left(2^k+\sigma \right),{\gamma }_k=\left(2^k+\theta \right)/2^k$,

$$F^{\left(2^k-1\right)}\left(x\right)=\mathit{Pe}\cdot S\left(x\right)+\frac{2^{k+1}+\mathit{Pe}\cdot \left(E-W\right)}{E-W}\frac{{\left(1-{\tau }_k\right)}^2}{1-{\tau }_k^{2^k}}\sum _{j=1}^{2^k-1}\sum _{i=1}^j{\tau }_k^{i-1}\cdot S\left(x+j\frac{x-W}{2^k}\right)-$$

$$\frac{2^{k+1}}{E-W}\frac{{\left(1-{\gamma }_k\right)}^2}{1-{\gamma }_k^{2^k}}\sum _{j=1}^{2^k-1}\sum _{i=1}^j{\gamma }_k^{i-1}\cdot S\left(x+\left(2^k-j\right)\frac{E-x}{2^k}\right).$$

Figures 9 and 10 show graphs of approximate solutions to the problem (18), obtained by (30) for $W=0,E=1,$with different moving nodes.

Figure 9.

Pe = 20, ${Ф}_W=0,{Ф}_E=1,S\left(x\right)=0$. Approximate solutions of the problem. Dotted—at k = 1, dotted—k = 2, dotted-dotted—k = 3, long dotted—k = 4, rarely dotted—k = 5. The solid line is the exact solution.

Figure 10.

Pe = 20, ${Ф}_W=0,{Ф}_E=0,S\left(x\right)=x,$. Approximate solutions of the problem. Dotted—at k = 1, dotted—k = 2, dotted-dotted—k = 3, long dotted—k = 4, rarely dotted—k = 5. The solid line is the exact solution.

It can be seen from the graphs that the approximate solutions give good results.

Remark. When obtaining many point-moving nodes, we proceeded from the upwind scheme. It was possible to proceed from the other three-point schemes.

2.3.2 Analytical control volume method for a one-dimensional convective-diffusion problem

It is known that differential equations are obtained on the basis of the integral conservation law. Therefore, discretization of the equations can be carried out using the approximation of integral conservation laws. This method is called the Finite Volume Method. Another name for the method is integro-interpolation.

Consider a one-dimensional convective-diffusion equation on a finite interval with boundary conditions in the form:

$$\frac{d}{\mathit{dx}}\left(\mathit{\rho uФ}\right)=\frac{d}{\mathit{dx}}\left(Г\frac{\mathit{dФ}}{\mathit{dx}}\right)+S\left(x\right)$$E31

$$Ф\left(W\right)={Ф}_W,Ф\left(E\right)={Ф}_E$$E32

where u is the flow velocity in the x direction, $\rho$ is the flow density, $Г$ is the diffusion coefficient, $S\left(x\right)$ is a given function (source), $Ф$ an unknown function. It follows from the continuity equation that $F=\mathit{\rho u}=\mathit{\text{const}}$.

Consider Eq. (31) into segments $\left[W,E\right]$. To obtain an approximate analytical solution to the problem using the control volume method, we take an arbitrary point $x\in \left[W,E\right]$ and control volume [w, e] (Figure 11). Let us assume that the face w is located in the middle between the points W and x, and the face e is in the middle between the points x and E. Integrating Eq. (31) over the control volume and replacing the derivatives with the upwind scheme, we obtain the zero approximation.

Figure 11.

Control volume [w, e].

$$\left({а}_{Е}+{а}_W\right){Ф}^0={а}_{Е}{Ф}_E^0+{а}_W{Ф}_W^0+\frac{E-W}{2}\cdot S\left(x\right)$$E33

Here $a_E=\frac{{Г}_{е}}{Е-х}+max\left(-F_e,0\right);a_W=\frac{{Г}_w}{x-W}+max\left(F_w,0\right).$ Since $x\in \left[W,E\right]$ an arbitrary point from (33), we can determine ${Ф}^0$ and obtain an approximate analytical solution to problem (31).

Note that, from (33) it follows that, in the absence of a source ( $S\left(x\right)\equiv 0$) on the segments $\left[W,E\right]$, the function is monotone.

To improve the approximate solution, we take additional nodes: $x_1=\frac{x+W}{2},x_2=\frac{x+E}{2}$. Let us write an upwind scheme of type (33) for the segment $\left[W,x\right],\left[x_1,x_2\right]$, and $\left[x,E\right]$. We get a system of three equations. We exclude the resulting system ${Ф}^1\left(x_1\right)$, ${Ф}^1\left(x_2\right)$ and as a result, we get an improved scheme:

$$\begin{array}{l}\left[\frac{{\beta }_1^{+}}{1+{\tau }_1}+\frac{{\alpha }_1^{-}}{1+{\gamma }_1}\right]{Ф}^1=\frac{{\beta }_1^{+}}{1+{\tau }_1}{Ф}_W^1+\frac{{\alpha }_1^{-}}{1+{\gamma }_1}{Ф}_E^1+\frac{E-W}{4}\cdot S\left(x\right)+\\ \frac{1}{1+{\tau }_1}\cdot \frac{x-W}{2}\cdot S\left(W+\frac{x-W}{2}\right)+\frac{1}{1+{\gamma }_1}\cdot \frac{E-x}{2}\cdot S\left(x+\frac{E-x}{2}\right).\end{array}$$E34

where ${\tau }_1=\frac{{\beta }_1^{-}}{{\beta }_1^{+}},{\gamma }_1=\frac{{\alpha }_1^{+}}{{\alpha }_1^{-}},{\beta }_1^{-}=2D_W+F^{-},{\beta }_1^{+}=2D_W+F^{+},$${\alpha }_1^{-}=2D_E+F^{-},$ ${\alpha }_1^{+}=2D_E+F^{+},$ $D_E=Г/\left(E-x\right),D_W=Г/\left(x-W\right),$$F^{-}=max\left(-F,0\right),$$F^{+}=max\left(F,0\right)$.

In (34), ${Ф}^1$ is the improved value of the unknown function at the nodal point x$\left({Ф}_W^1\equiv {Ф}_W,{Ф}_E^1\equiv {Ф}_E\right)$.

where in (34), the improved value of the unknown function at the nodal point is x.

Solving (34) with respect to, we obtain an improved analytical solution. Again, to improve the solution, we proceed in a similar way: we write the scheme (34) for the segment $\left[W,x\right],\left[x_1,x_2\right]$ and $\left[x,E\right]$, and eliminate the unknowns at the points $x_1$ and $x_2$, and so on. Continuing this process, we get.

$$\begin{array}{l}\left[\frac{\left(1-{\tau }_k\right){\beta }_k^{+}}{1-{\tau }_k^{2^k}}+\frac{\left(1-{\gamma }_k\right){\alpha }_k^{-}}{1-{\gamma }_k^{2^k}}\right]{Ф}^k=\frac{\left(1-{\tau }_k\right){\beta }_k^{+}}{1-{\tau }_k^{2^k}}{Ф}_W^k+\frac{\left(1-{\gamma }_k\right){\alpha }_k^{-}}{1-{\gamma }_k^{2^k}}{Ф}_E^k+\frac{E-W}{2^{k+1}}\cdot S\left(x\right)\\ +\frac{1-{\tau }_k}{1-{\tau }_k^{2^k}}\cdot \frac{x-W}{2^k}\cdot \sum _{j=1}^{2^k-1}\sum _{i=1}^j{\tau }_k^{i-1}S\left(W+j\frac{x-W}{2^k}\right)\\ +\frac{1-{\gamma }_k}{1-{\gamma }_k^{2^k}}\cdot \frac{E-x}{2^k}\cdot \sum _{j=1}^{2^k-1}\sum _{i=1}^j{\gamma }_k^{i-1}S\left(x+\left(2^k-j\right)\frac{E-x}{2}\right).\end{array}$$E35

Here ${\tau }_k=\frac{{\beta }_k^{-}}{{\beta }_k^{+}},{\gamma }_k=\frac{{\alpha }_k^{+}}{{\alpha }_k^{-}},{\beta }_k^{-}=2^k{D}_W+F^{-},{\beta }_k^{+}=2^k{D}_W+F^{+},$${\alpha }_k^{-}=2^k{D}_E+F^{-},$ ${\alpha }_k^{+}=2^k{D}_E+F^{+}$.

In (35), ${Ф}^k$ is the improved value of the unknown function at the nodal point x$\left({Ф}_W^k\equiv {Ф}_W,{Ф}_E^k\equiv {Ф}_E\right)$. Solving Eq. (35) with respect to ${Ф}^k$, we obtain an approximate analytical solution of the original problem.

Examples.

Figure 12 shows solutions to the problem (31) $Г=\mathit{\text{const}},R=\mathit{\rho u}/Г=20,S\left(x\right)=0$ for segments $\left[0;1\right]$ with boundary conditions ${Ф}_W=0,{Ф}_E=1.$ Figure 13 shows solutions to the problem (31) $R=50,S\left(x\right)=5cos4x$ for segments $\left[0;1\right]$ with boundary conditions ${Ф}_W=0,{Ф}_E=1.$ The graph shows that, as k increases, the approximate solutions approach the exact one.

Figure 12.

Comparison of the approximate solutions for S(x) = 0. Continuous line is exact, point—k = 0, dotted line—k = 1, dot-dotted line—k = 2, long dotted line—k = 4, rare dotted line—k = 6.

Figure 13.

Comparison of the approximate solutions for S(x) = 5cos4x. Continuous line is exact, point—k = 0, dotted line—k = 1, dot-dotted line—k = 2, long dotted line—k = 4, rare dotted line—k = 6.

It can be seen from the graphs that, starting from k = 6, the exact and approximate solutions visually coincide.

It is interesting to compare the analytical solution obtained by the finitely different method (30) and the control volume method (R = 10, S(x) = 0).

From Figures 14 and 15, it can be seen that the solution obtained by the control volume method is preferable.

Figure 14.

Approximate solution for k = 0. Solid curves are exact solutions, point curves are the finite difference method; dotted lines—control volume method.

Figure 15.

Approximate solution for k = 1. Solid curves are exact solutions, point curves are the finite difference method; dotted lines—control volume method.

2.3.3 Improving accuracy with Richardson extrapolation

Using the method described, we can improve the accuracy of approximate solutions to the problem [39]. Linear combination $Q^3\left(x\right)=-\frac{1}{3}{U}^1\left(x\right)+\frac{4}{3}{U}^3\left(x\right)$ is more accurately approximates the solution. With a linear combination of $U^1\left(x\right),U^3\left(x\right)$ and $U^7\left(x\right)$ in the form $Q^7\left(x\right)=\frac{1}{45}{U}^1\left(x\right)-\frac{4}{9}{U}^3\left(x\right)+\frac{64}{45}{U}^7\left(x\right)$, we obtain a more refined solution to the problem [39].

Figure 16 shows graphs of approximate solutions to the problem (31) obtained by Richardson’s extrapolation for $W=0,E=1.$ The solid line in Figure 16–19 is the exact solution.

Figure 16.

${Ф}_W=0,$ ${Ф}_E=1,S\left(x\right)=0,$ $\mathit{Pe}=20.$ Comparisons of solutions. Dotted line is $U^3\left(x\right)$, point line—$Q^3\left(x\right)$, dot-dotted line—$U^7\left(x\right)$, long dotted line—$Q^7\left(x\right)$.

Figure 17.

${Ф}_W=0,$ ${Ф}_E=0,S\left(x\right)=x,$ $\mathit{Pe}=20.$ Comparisons of solutions. Dotted line is $U^3\left(x\right)$, point line—$Q^3\left(x\right)$, dot-dotted line—$U^7\left(x\right)$, long dotted line—$Q^7\left(x\right)$.

Figure 18.

${Ф}_W=0,$ ${Ф}_E=1,S\left(x\right)=0,$ $\mathit{Pe}=20.$ Comparisons of solutions. Dotted line—$U^{15}\left(x\right)$, dotted line—$Q^{15}\left(x\right)$.

Figure 19.

${Ф}_W=0,$ ${Ф}_E=0,S\left(x\right)=x,$ $\mathit{Pe}=20.$ Comparisons of solutions. Dotted line—$U^{15}\left(x\right)$, dotted line—$Q^{15}\left(x\right)$.

Figures 16–19 allow us to state that Richardson’s extrapolation makes it possible to obtain a more refined solution to the problem.

2.4 Moved node method for non-stationary problems

In the previous paragraphs, the application of the MNM for ordinary differential equations has been considered. Here we consider the application of the MNM for parabolic equations.

An example of a problem that leads to a parabolic partial differential equation is the problem of heat transfer along a long rod, described by the heat transfer (or diffusion) equation.

The problem is to find a function U(x,t) in the region Ω = {(x,t) | W ≤ x ≤ E, 0 ≤ t ≤ T} satisfying the equation.

$$\frac{\partial U}{\partial t}=A\frac{{\partial }^{2}U}{\partial x^2}+f\left(x,t\right),A>0$$E36

initial condition

$$U\left(x,0\right)=U^0\left(x\right)$$

and boundary conditions of the first kind

$$U\left(W,t\right)=U_W\left(t\right);U\left(E,t\right)=U_E\left(t\right).$$

Let us take an arbitrary point Ω in the area $\left(x,t\right)\in \mathrm{\Omega }$ (Figure 20). We will accept this point as moving. We approximate (36) by the implicit scheme

Figure 20.

The region of solution.

$$\frac{Y\left(x,t\right)-U^0\left(t\right)}{t}=A\frac{2}{E-W}\left(\frac{U_E\left(t\right)-Y\left(x,t\right)}{E-x}-\frac{Y\left(x,t\right)-U_W\left(t\right)}{x-W}\right)+f\left(x,t\right),$$E37

In (37), $Y\left(x,t\right)$ is an approximate analytical solution. When the point runs through $\mathrm{\Omega }$, we get a solution in the area under consideration. From (37), we get

$$\begin{array}{l}Y\left(x,t\right)=\frac{\left(E-x\right)\left(x-W\right)}{2\mathit{At}+\left(E-x\right)\left(x-W\right)}{U}^0\left(t\right)+\frac{2\mathit{At}\left[U_E\left(t\right)\left(x-W\right)+U_W\left(t\right)\left(E-x\right)\right]}{2\mathit{At}+\left(E-x\right)\left(x-W\right)}+\\ \frac{\left(E-x\right)\left(x-W\right)t}{2\mathit{At}+\left(E-x\right)\left(x-W\right)}f\left(x,t\right).\end{array}$$E38

Consider examples.

2.4.1 Test problems

Let us consider Eq. (36) $0<x<1$ with conditions $U^0\left(x\right)=x,$ $U_W\left(t\right)=0,U_E\left(t\right)=e^{-t},$ $f\left(x,t\right)=-x{e}^{-t}.$ Exact solution of problem is $U\left(x,t\right)=x{e}^{-t}.$ Figure 21 presents a comparison of the exact and approximate solutions for the cross-section $x=0,5$ and $x=0,2$. The solid lines are the exact solution. Figure 21 shows the closeness of the exact and approximate solutions.

Figure 21.

Solution comparison of the exact and approximate solutions for the sections x = 0, 5 and x = 0.2.

Let us consider Eq. (36)$0<x<1$ with conditions $U^0\left(x\right)=sin\mathit{\pi x}+x^2,$ $U_W\left(t\right)=0,U_E\left(t\right)=1,$ $f\left(x,t\right)=-sin\mathit{\pi x}{e}^{-t}+{\pi }^{2}sin\mathit{\pi x}{e}^{-t}-2.$ Exact solution of problem presents a comparison of the exact and approximate solutions for the sections $x=0,1$, $x=0,5$, and $x=0,8$. Figure 22 shows the closeness of the exact and approximate solutions.

Figure 22.

Solution comparison of the exact and approximate solutions for the sections x = 0, 1, x = 0.5 and x = 0.8.

2.4.2 Unsteady flow of a viscous fluid between parallel walls

As a practical example, consider an unsteady flow of a viscous fluid between parallel walls. Let a viscous fluid fill the entire space between horizontal planes located at a certain distance from each other. Let the lower plane be stationary all the time, and the upper one starts to move to the right at a constant speed. We neglect the action of gravity and assume that the pressure is constant everywhere. The flow is assumed to be directed parallel to the x-axis. Then the equation of motion of a viscous fluid in dimensionless variables has the form.

$$\frac{\mathit{\partial u}}{\mathit{\partial t}}=\frac{{\mathit{\partial }}^{2}u}{\mathit{\partial }{y}^2}$$E39

The exact solution of the equation under the conditions:

$$u\left(0,y\right)=0,u\left(t,0\right)=0,u\left(t,1\right)=1$$

looks like:

$$u=y+\frac{2}{\pi }\sum _{k=1}^{\infty }\frac{{\left(-1\right)}^k}{k}sin\left(\mathit{k\pi y}\right)exp\left(-k^2{\pi }^{2}t\right)$$E40

Let us replace (39) with the difference relation:

$$\frac{u_1\left(t,y\right)-u_1\left(0,y\right)}{t-0}=\frac{2}{1-0}\left[\frac{u_1\left(t,1\right)-u_1\left(t,y\right)}{1-y}-\frac{u_1\left(t,y\right)-u_1\left(t,0\right)}{y-0}\right]$$

From here, considering the boundary conditions, we obtain an approximate solution in the form [15]:

$$u_1\left(t,y\right)=\frac{2\mathit{ty}}{y\left(1-y\right)+2t}$$E41

Note that, $\underset{t\to \infty }{lim}u\left(t,y\right)=\underset{t\to \infty }{lim}{u}_1\left(t,y\right)=y.$

There is another approach to obtaining an approximate solution to Eq. (39). We replace (39) with the following equation:

$$\frac{{\mathit{du}}_2}{\mathit{dt}}=\frac{2}{1-0}\left[\frac{u_2\left(t,1\right)-u_2\left(t,y\right)}{1-y}-\frac{u_2\left(t,y\right)-u_2\left(t,0\right)}{y-0}\right]$$E42

Considering Eq. (42) y as a parameter, and solving it, we get

$$u_2\left(t,y\right)=y\left(1-exp\left(\frac{2t}{y\left(1-y\right)}\right)\right)$$E43

This shows that partial approximation gives the best result (Figures 23 and 24).

Figure 23.

Comparison of (39) and (41) (red light approx. solution).

Figure 24.

Comparison of (39), (41), and (43) on the section y = 0.8. Blue line on (41) black on (43), red fine.

2.4.3 Non-stationary convection-diffusion differential equation

Consider the equation

$$\frac{\partial Ф}{\partial t}+\frac{\partial Ф}{\partial x}=\frac{1}{\mathit{Pe}}\frac{{\partial }^2Ф}{\partial x^2}+f\left(x,t\right),$$E44

Under appropriate boundary and initial conditions. We approximate Eq. (44) as follows

$$\begin{array}{l}\frac{U\left(x,t\right)-U\left(x,0\right)}{t}+\frac{U\left(x,t\right)-U\left(W,t\right)}{x-W}=\\ \frac{1}{\mathit{Pe}}\frac{2}{\left(E-W\right)}\left(\frac{U\left(E,t\right)-U\left(x,t\right)}{E-x}-\frac{U\left(x,t\right)-U\left(W,t\right)}{x-W}\right)+f\left(x,t\right),\end{array}$$E45

(45) is an implicit difference scheme with a moving node. In this case, the convective term was approximated by the scheme against the flow, and the diffusion term, as usual, with the second order of accuracy.

The comparison obtained with the help of (45) of the approximate solution with the exact solution (44) under the conditions $U^0\left(x\right)=x^2+x,$ $U_W\left(t\right)=0,$ $U_E\left(t\right)=1+e^{-\mathit{Pe}t},$ $f\left(x,t\right)=\left(1+\mathit{Pe}x\right)e^{-\mathit{Pe}t}+2x-2/\mathit{Pe}$ is shown in Figure 25. Exact solution $(W=0,$$E=1)Ф\left(x,t\right)=x^2+e^{-\mathit{Pe}t}x$. In Figure 25, the solid curves are the exact solution, the dotted curves are the approximate solution, and the graphs correspond to the sections $x=0,1;x=0,5;x=0,8$. Figure 26 same results corresponding to $t=1;t=5;t=10.$

Figures 25 and 26 show the acceptability of the approximate solution for the MNM.

Figure 25.

Comparison solution $\mathit{Pe}=0,5$.

Figure 26.

Comparison solution $\mathit{Pe}=0,5$.

It should be noted that with increasing $\mathit{Pe}$ the discrepancy between the exact and approximate solutions increases. On Figures 27 and 28 compare the same problem with $\mathit{Pe}=2$.

Figure 27.

Comparison solution. $\mathit{Pe}=2$.

Figure 28.

Comparison solution. $\mathit{Pe}=2$.

Thus, the MMN makes it possible to obtain an approximate analytical solution.

2.5 MNM for two-dimensional boundary value problems

Now let us consider the application of MMN to two-dimensional boundary value problems to obtain rough approximate solutions of DE.

Consider a convex closed two-dimensional region (Figure 29). P point inside area. If P changes position inside the region, the boundary points $E,N,W,S$ change their positions while being on the border of the region.

Figure 29.

The convex closed two-dimensional region.

When studying stationary processes of various physical natures (oscillations, heat conduction, diffusion, hydrodynamics, etc.), one usually leads to equations of the elliptic type. The most common equation of this type is the Poisson equation.

There are various approximate-analytical and numerical methods for the equation of mathematical physics.

Consider the two-dimensional Poisson equation in a rectangle $\left(x,y\right)\in \left[W,E\right]\times \left[S,N\right]$

$$\mathrm{\Delta }U\left(x,y\right)=f\left(x,y\right),$$E46

with boundary conditions.

$$U\left(W,y\right)=U_W\left(y\right),U\left(E,y\right)=U_E\left(y\right),U\left(x,S\right)=U_S\left(x\right),U\left(x,N\right)=U_N\left(x\right).$$E47

Take an arbitrary point in the rectangle approximation of second-order partial derivatives:

$$\frac{{\partial }^{2}U}{\partial x^2}\approx \frac{2}{E-W}\left[\frac{U_E\left(y\right)-u\left(x,y\right)}{E-x}+\frac{u\left(x,y\right)-U_W\left(y\right)}{x-W}\right],$$E48

$$\frac{{\partial }^{2}U}{\partial y^2}\approx \frac{2}{N-S}\left[\frac{U_N\left(x\right)-u\left(x,y\right)}{N-y}+\frac{u\left(x,y\right)-U_S\left(x\right)}{y-S}\right].$$E49

Substituting (46) and (49) into (46), and solving, the resulting equation with respect to $u\left(x,y\right),$ we have

$$\begin{array}{l}u\left(x,y\right)=\frac{1}{\left(E-x\right)\left(x-W\right)+\left(N-y\right)\left(y-S\right)}\cdot \left[\frac{\left(N-y\right)\left(y-S\right)}{E-W}\left(\left(x-W\right)U_E+\left(E-x\right)U_W\right)\right.+\\ \left.+\frac{\left(E-x\right)\left(x-W\right)}{N-S}\left(\left(y-S\right)U_N+\left(N-y\right)U_S\right)\right]+\frac{\left(E-x\right)\left(x-W\right)\left(N-y\right)\left(y-S\right)}{2\left(\left(E-x\right)\left(x-W\right)+\left(N-y\right)\left(y-S\right)\right)}f\left(x,y\right)\end{array}$$E50

This is the approximate analytical solution of the Poisson equation in a rectangle. (49) satisfies the boundary conditions. Due to the fact that (48) and (49) is an approximate relation for the approximation of the second derivatives (50) is an approximate solution. Nevertheless, (50) gives an acceptable solution to many practical problems.

Consider examples.

2.5.1 Test problems

1. Consider the Laplace equation in a rectangle $\left[0,1\right]\times \left[0,1\right]$ with boundary conditions $U\left(0,y\right)=0,U\left(1,y\right)=y,U\left(x,0\right)=0,U\left(x,1\right)=x$. The exact solution to this problem is $U\left(x,y\right)=\mathit{xy}$. If we use the approximate solution (50), we obtain. In this case, the approximate solution coincides with the exact solution.

2. The function $U\left(x,y\right)=x^2+y^2$, with boundary conditions $U\left(0,y\right)=-y^2,$ $U\left(1,y\right)=1-y^2,U\left(x,0\right)=x^2,U\left(x,1\right)=x^2-1$ satisfies the Laplace equation. Relation (50) gives us an identical result.

3. The function $U\left(x,y\right)=ln\left(x^2+y^2\right)$, with boundary conditions $U\left(1,y\right)=ln\left(1-y^2\right),$$U\left(2,y\right)=ln\left(4+y^2\right),U\left(x,0\right)=ln\left(x^2\right),U\left(x,1\right)=ln\left(x^2+1\right)$ in the region $\left[1,2\right]\times \left[0,1\right]$, satisfies the Laplace equation. An approximate solution based on (50) gives $\begin{array}{l}u\left(x,y\right)=\frac{1}{\left(2-x\right)\left(x-1\right)+y\left(1-y\right)}\left[y\left(1-y\right)\right.\left(\left(x-1\right)ln\left(4+y^2\right)+\left(2-x\right)ln\left(1+y^2\right)\right)+\left.\left(2-x\right)\left(x-1\right)\left(yln\left(1+x^2\right)++\left(1-y\right)ln\left(x^2\right)\right)\right]\end{array}$If we compare the exact and approximate solutions in the area under consideration at points $x_i=1+\mathit{ih},y_j=\mathit{jh},i,j=1,2,\dots n$ with a step $h=0,1$ for the maximum difference, we obtain $0,0011$.

4. The function $U\left(x,y\right)=x^3-y^3$, with boundary conditions, $U\left(0,y\right)=-y^3,U\left(1,y\right)=1-y^3,$$U\left(x,0\right)=x^3,U\left(x,1\right)=x^3-1$ satisfies Eq. (46) for $f\left(x,y\right)=6x-6y$. Based on (50), the approximate solution has the form:

$$u\left(x,y\right)=\frac{\mathit{xy}\left(y-x\right)+y^4\left(1-y\right)+x^4\left(x-1\right)}{y\left(y-1\right)+x\left(x-1\right)}$$

The maximum absolute difference between the exact and approximate solutions calculated by points $x_i=1+\mathit{ih},y_j=\mathit{jh},i,j=1,2,\dots n,h=0,1$ is 0.048.

If we approximate the right side based on the control volume [35], the approximate solution has the form:

$$u\left(x,y\right)=\frac{\mathit{xy}\left(x-y\right)\left(1-3\left(x+y\right)+3\mathit{xy}\right)+2y^4\left(1-y\right)+2x^4\left(x-1\right)}{2\left[y\left(y-1\right)+x\left(x-1\right)\right]}$$

and the maximum absolute difference between the exact and approximate solutions calculated by points $x_i=1+\mathit{ih},y_j=\mathit{jh},i,j=1,2,\dots n,h=0,1$ is 0.024.

2.5.2 Flow in an ellipsoidal pipe

The equation describing the one-dimensional flow in an ellipsoidal tube of a viscous fluid has the form:

$$\frac{{\partial }^{2}U}{\partial y^2}+\frac{{\partial }^{2}U}{\partial z^2}=-\frac{\mathrm{\Delta }p}{\mu l}$$E51

Here u is the flow rate, μ is the flow viscosity, Δp/l (Δp/l = const) is the pressure drop. Eq. (51) is considered in the area $\frac{y^2}{a^2}+\frac{z^2}{b^2}\le 1$ (section of an ellipsoidal pipe, Figure 30), and the boundary condition is the no-slip condition (U = 0).

Figure 30.

Ellipsoidal pipe section.

Eq. (51) is replaced by the difference

$$\frac{2}{y_E-y_W}\left(\frac{U_E-u}{y_E-y}-\frac{u-U_W}{y-y_W}\right)+\frac{2}{z_N-z_S}\left(\frac{U_N-u}{z_N-z}-\frac{u-U_S}{z-z_S}\right)=-\frac{\mathrm{\Delta }p}{\mathit{\mu l}}.$$

Hence, given that

$$z_N=b\sqrt{1-y^2/a^2},z_S=-b\sqrt{1-y^2/a^2},y_E=a\sqrt{1-z^2/b^2},y_W=-a\sqrt{1-z^2/b^2}.$$

we get

$$u=\frac{a^2{b}^2}{2\left(a^2+b^2\right)}\left(1-\frac{y^2}{a^2}-\frac{z^2}{b^2}\right)\frac{\mathrm{\Delta }p}{\mathit{\mu l}}$$

coinciding with the exact solution.

2.5.3 Two-dimensional temperature field in a solid

This problem is reduced to solving an equation $\mathrm{\Delta }T=0,$ with boundary conditions

$$T\left(0,y\right)=0,T\left(1,y\right)=0,T\left(x,0\right)=T_0,T\left(x,1\right)=0.$$

Exact solution to the problem

$$T\left(x,y\right)=\sum _{n=1}^{\infty }{A}_{n}sin\left(\mathit{n\pi x}\right)\mathit{sh}\left[\mathit{n\pi }\left(y-1\right)\right]$$

where ${А}_n=\frac{2T_0}{\mathit{n\pi }}\frac{\left[{\left(-1\right)}^n-1\right]}{\mathit{sh}\left(\mathit{n\pi }\right)}$.

Approximate solution

$$T\left(x,y\right)=\frac{\left(1-x\right)x\left(1-y\right)T_0}{x\left(1-x\right)+y\left(1-y\right)}.$$

The maximum absolute difference between the exact and approximate solutions calculated by points $x_i=\mathit{ih},y_j=\mathit{jh},i,j=1,2,\dots n,h=0,01$ is 0.015.

Thus, the method presented here allows for obtaining solutions to Dirichlet problems. To improve the solution, the mesh refinement technique can be used.

To increase the accuracy, increase the number of moved nodes. When the number of nodes to be moved is four, we get.

$$u_4\left(x,y\right)={\left\{1-\frac{4}{A+B}\left[\frac{B^2}{8}\left(\frac{x}{{\left(\frac{1-x}{2}\right)}^2+B}+\frac{1-x}{{\left(\frac{x}{2}\right)}^2+B}\right)+\frac{A^2}{8}\left(\frac{y}{{\left(\frac{1-y}{2}\right)}^2+A}+\frac{1-y}{{\left(\frac{y}{2}\right)}^2+A}\right)\right]\right\}}^{-1}$$

$$\begin{array}{l}\times \frac{4}{A+B}\left\{\frac{B}{2}\left[\frac{1}{2}\frac{\frac{x}{2}\left({\mathit{Bu}}_b\left(y\right)+{\left(\frac{1-x}{2}\right)}^2\left({\mathit{yu}}_d\left(\frac{1+x}{2}\right)+\left(1-y\right)u_c\left(\frac{1+x}{2}\right)\right)\right)}{{\left(\frac{1-x}{2}\right)}^2+B}+\right.\right.\\ +\left.\frac{1}{4}\frac{\left(1-x\right)\left({\mathit{Bu}}_a\left(y\right)+{\left(\frac{x}{2}\right)}^2\left({\mathit{yu}}_d\left(\frac{x}{2}\right)+\left(1-y\right)u_c\left(\frac{x}{2}\right)\right)\right)}{{\left(\frac{x}{2}\right)}^2+B}\right]+\end{array}$$

$$\begin{array}{l}+\frac{A}{2}\left[\frac{1}{2}\frac{y\left(\frac{A}{2}{u}_d\left(x\right)+{\left(\frac{1-y}{2}\right)}^2\left({\mathit{xu}}_b\left(\frac{1+y}{2}\right)+\left(1-x\right)u_a\left(\frac{1+y}{2}\right)\right)\right)}{{\left(\frac{1-y}{2}\right)}^2+A}+\right.\\ +\left.\left.\frac{1}{2}\frac{\left(1-y\right)\left({\mathit{Au}}_c\left(x\right)+{\left(\frac{y}{2}\right)}^2\left({\mathit{xu}}_b\left(\frac{y}{2}\right)+\left(1-x\right)u_a\left(\frac{y}{2}\right)\right)\right)}{{\left(\frac{y}{2}\right)}^2+A}\right]\right\}\end{array}$$

The maximum absolute difference between the exact and approximate solutions, calculated by points $x_i=\mathit{ih},y_j=\mathit{jh},i,j=1,2,\dots n,h=0,1$, is 0.14 according to the formula with one moving node, and when calculating with five moving nodes, it is 0.07.

2.5.4 Flow in a rectangular pipe

Eq. (51) also describes the flow of an incompressible viscous fluid in a rectangular pipe. Let us denote the height of the rectangle parallel to the axis $\mathit{Oz}$ as $2h$, and the base parallel to the axis $\mathit{Oy}$ as – $2\mathit{\sigma h}$, where $\sigma$ is any positive constant. We draw the axis through the center of the rectangle and direct it downstream.

Let us transform Eq. (51) into a dimensionless form. For the scale of lengths, we take the height, $h$, and for the scale of speeds—the value $h^2/\mu \cdot \mathrm{\Delta }p/l$. We introduce the following dimensionless quantities:

$$Y=y/h,Z=z/h,V=\mathit{U\mu }l/\left(h^2\mathrm{\Delta }p\right)$$

Substituting into (51), we obtain

$$\frac{{\partial }^{2}V}{\partial Y^2}+\frac{{\partial }^{2}V}{\partial Z^2}=-1$$E52

Boundary conditions for (52)

$$V\left(Y,-1\right)=0,V\left(Y,1\right)=0,V\left(-\sigma ,Z\right)=0,V\left(\sigma ,Z\right)=0$$E53

Eq. (52) is replaced by a difference equation and taking into account the boundary condition (53) we have

$$\frac{2}{2\sigma }\left(\frac{-V}{\sigma -Y}-\frac{V}{Y+\sigma }\right)+\frac{2}{1+1}\left(\frac{-V}{1-Z}-\frac{V}{Z+1}\right)=-1.$$

From here we determine the approximate analytical solution:

$$V=\frac{1}{2}\frac{\left({\sigma }^2-Y^2\right)\left(1-Z^2\right)}{1-Z^2+{\sigma }^2-y^2}$$E54

The exact solution of the problem has the form:

$$u=\frac{16{\sigma }^2}{{\pi }^3}\sum _{n=0}^{\infty }\frac{{\left(-1\right)}^n}{{\left(2n+1\right)}^3}\left[1-\frac{\mathit{ch}\left(\frac{2n+1}{2}\frac{\pi }{\sigma }Y\right)}{\mathit{ch}\left(\frac{2n+1}{2}\frac{\pi }{\sigma }\right)}\right]cos\left(\frac{2n+1}{2}\frac{\pi }{\sigma }Z\right)$$

Figure 31 shows a comparison of the exact and approximate solutions on the cross-section $x=0$ for $\sigma =1$. The maximum absolute difference between the exact and approximate solutions is 0.045.

Figure 31.

Comparison of the solution on the section according to (54).

To increase the accuracy of the approximate solution in Eq. (52), we approximate only one of the terms. For example, we approximate Eq. (52) as follows:

$$\frac{2}{2\sigma }\left(\frac{-V}{\sigma -Y}-\frac{V}{Y+\sigma }\right)+\frac{{\partial }^{2}V}{\partial Z^2}=-1$$E55

We got an ordinary differential equation, we consider the variable $Y$ in Eq. (55) as a parameter. We solve Eq. (55) with constant coefficients, considering the boundary conditions, we find an approximate solution

$$V=C_{1}exp\left(\sqrt{k}Z\right)+C_{2}exp\left(-\sqrt{k}Z\right)+\frac{1}{k}.$$E56

Here $k=2/\left(\left(\sigma -Y\right)\left(Y+\sigma \right)\right),$${С}_2=-\frac{1}{k}\frac{exp\left(\sqrt{k}\right)-exp\left(-\sqrt{k}\right)}{exp\left(2\sqrt{k}\right)-exp\left(-2\sqrt{k}\right)},$ ${С}_1=-{С}_2-exp\left(2\sqrt{k}\right)-\frac{1}{k}exp\left(-\sqrt{k}\right).$

Figure 31 shows a comparison of the exact approximate solution obtained based on (56) on the cross-section $x=0$ at $\sigma =1$. A comparison of Figures 31 and 32 shows that the calculation by formula (56) gives a more accurate result. The maximum absolute difference between the exact and approximate solutions is equal to that obtained by (56) and equals 0.024. In Figures 31 and 32 solid curves are the exact solution.

Figure 32.

Comparison of the solution on the section according to (56).

2.5.5 Flow at the inlet section of the pipe

With appropriate simplifications, the flow of a viscous incompressible fluid in a dimensionless form is described by the following differential equation:

$$u\frac{\partial u}{\partial x}=-\frac{1}{Re}N+\frac{1}{Re}\left(\frac{{\partial }^{2}u}{\partial y^2}+\frac{{\partial }^{2}u}{\partial x^2}\right),$$E57

Here, N = −12 is the pressure drop, Re is the Reynolds number. The equations are considered in the area $D:-0,5<x<0.5,0<y<L.$ (Figure 33). Boundary conditions for (57):

Figure 33.

Coordinate systems and the region of solution.

$$\begin{array}{l}u\left(0,y\right)=1;u\left(L,y\right)=1,5\left(1-4y^2\right);\\ u\left(x,-0,5\right)=0;u\left(x,0.5\right)=0.\end{array}$$

The convective term is linearizable

$$u\frac{\partial u}{\partial x}\approx \frac{\partial u}{\partial x}.$$

Approximating in (57) by the liquid volume $\left(y+0,5\right)/2<y<\left(y-0,5\right)/2,$ we obtain an ordinary equation, solving which we obtain an approximate solution:

$$u=C_{1}exp\left(k_{1}x\right)+C_{2}exp\left(k_{2}x\right)-\frac{1-4y^2}{8}N,$$E58

where $k_{1,2}=\frac{Re}{2}\left(1\pm \sqrt{1+\frac{32}{{Re}^2\left(1-4y^2\right)}}\right).$

For comparison, solutions (57) were also made with the numerical method.

Figure 34 shows the velocity profiles obtained on the basis of an approximate solution. The solid curve to the section $x=0,1,$ and the pointed curve to $x=0,5,$ the dotted one corresponds to the section $x=3.$ Figure 35 shows a comparison of the approximate and numerical solution of Eq. (57). The solid lines correspond to the solution (58), and the dotted lines correspond to the numerical solution (velocity profiles are given for the cross-section $x=0,1$ and $x=0,5$).

Figure 34.

Approximate solution based on (58). Velocity profiles corresponding to sections $x=0,1;0,5;3.$ Re = 1, L = 5.

Figure 35.

Comparison of approximate and numerical solution. Re = 1, L = 5.

2.6 Solution of the flow problem in the combined region

Exact solution. Let a liquid flow in a flat pipe partially filled with a porous medium. The lower part of the horizontal pipe is filled with a porous medium of height h (pipe height H). Considering the flow to be one-dimensional and stationary, we obtain from the Rakhmatulin equation [16, 17], we obtain

$$\mu \frac{\mathit{du}}{\mathit{dy}}\left(f\frac{\mathit{du}}{\mathit{dy}}\right)-\mathit{Ku}=f\frac{\mathit{dp}}{\mathit{dx}}.$$E59

In (59) for the parameter K, we use the Kozeny-Karman relation as adopted in porous media:

$$K=\frac{\mu \cdot f^2}{k}.$$E60

where the $k=\frac{d^2{f}^3}{150{\left(1-f\right)}^2},$ permeability, d is the characteristic size of the porous medium.

Let us pass to dimensionless variables assuming $u=\overline{u}U,y=\overline{y}H,x=\overline{x}H,$ $p=\frac{\rho U^2}{Re}\overline{p}.$ Then Eq. (59) in dimensionless form for $f=\mathit{\text{const}},$ has the form:

$$\frac{d^2\overline{u}}{d{\overline{y}}^2}-A\overline{u}=\frac{d\overline{p}}{d\overline{x}}.$$E61

Here $A=180{\left(H/d\right)}^2{\left(1-f\right)}^2/f^2$.

In the free zone, the one-dimensional flow satisfies the equation

$$\frac{d^2\overline{u}}{d{\overline{y}}^2}=\frac{d\overline{p}}{d\overline{x}}.$$E62

In the future, in Eqs. (61) and (62), we release the dash above the variables.

Eq. (61) is considered when $0<y<h_0,$ and Eq. (62) $h_0<y<1.$ Equations are solved under the following boundary conditions.

No-slip conditions for Eq. (61) to the lower walls, and for Eq. (62) to the upper walls:

$$u\left(0\right)=0,u\left(1\right)=0.$$E63

In the inner boundary region, we set the conditions for the continuity of the flow and the equality of the shear stress:

$$u\left(h_0-0\right)=u\left(h_0+0\right),\frac{\mathit{du}\left(h_0-0\right)}{\mathit{dy}}=\frac{\mathit{du}\left(h_0+0\right)}{\mathit{dy}}.$$E64

It is easy to obtain an analytical solution of (61) and (62) under the given boundary conditions. Figure 36 shows an analytical solution. The dimensionless pressure difference is adopted $\frac{d\overline{p}}{d\overline{x}}=-12,$ so that it corresponds to the flow without a porous layer. The dotted line corresponds to the solution obtained with a porosity of 0.3, and the dotted-dotted line is 0.5.

Figure 36.

Exact solution: velocity distributions for different porosity values.

Numerical solution. Consider eq. (61) for the entire region and set

$$f=\left\{\begin{array}{c}\epsilon \mathit{при}0<y<h_0\\ 1\mathit{при}{h}_0\le y<1\end{array}\right..$$E65

In this case, Eq. (61) in the pure region takes the form (62). Thus, Eq. (61) can be used in the entire area, with porosity (65), while the interboundary conditions are satisfied automatically (in the porous layer, the porosity is taken equal to ε). For this purpose, a finite-difference approximation of Eq. (61) was compiled and calculated using the sweep method in the combined region. Figure 37 presents the results of numerical calculations (solid curves are the analytical solution, and point data are the numerical results). This shows that it is possible to perform a thorough calculation without highlighting the interboundary condition.

Figure 37.

Comparison of exact and numerical results for different porosity values.

Approximate analytical solution using a moving node. Using the moving node method, one can find an approximate analytical solution to the problem.

Eqs. (61) and (62) is approximated by difference relations:

$$\frac{2}{h_0}\left(\frac{u_G-u}{h_0-y}-\frac{u}{y}\right)-\mathit{Au}=\frac{\mathit{dp}}{\mathit{dx}}.$$E66

$$\frac{2}{1-h_0}\left(\frac{-u}{1-y}-\frac{u-u_G}{y-h_0}\right)=\frac{\mathit{dp}}{\mathit{dx}}.$$E67

In the difference Eqs. (66) and (67) no-slip boundary conditions are used. In Eqs. (66) and (67) $u_G$— the value of the unknown function on the inner boundary. To find $u_G$, we use the second interboundary condition (64). We put in (66) $y\to h_0-0,$ then we have

$$\frac{2}{h_0}\left({\left.\frac{\mathit{du}}{\mathit{dy}}\right|}_{h_0-0}-\frac{u}{h_0}\right)-{\mathit{Au}}_G=\frac{\mathit{dp}}{\mathit{dx}}.$$E68

If in (67) $y\to h_0+0$ then we have:

$$\frac{2}{1-h_0}\left(\frac{-u_G}{1-h_0}-{\left.\frac{\mathit{du}}{\mathit{dy}}\right|}_{h_0+0}\right)=\frac{\mathit{dp}}{\mathit{dx}}.$$E69

Using (64), we obtain

$$u_G=-\frac{h_0\left(1-h_0\right)}{2+h_0^2\left(1-h_0\right)}\frac{\mathit{dp}}{\mathit{dx}}.$$E70

Using (70) from (66) and (67) we determine the distribution of velocities in the porous

$$u=-\frac{y}{2+\mathit{Ay}\left(h_0-y\right)}\left(h_0-y+\frac{2\left(1-h_0\right)}{2+h_0^2\left(1-h_0\right)}\right)\frac{\mathit{dp}}{\mathit{dx}}.$$E71

and free zone

$$u=-\frac{1-y}{1-h_0}\left(\frac{\left(1-h_0\right)\left(y-h_0\right)}{2}+\frac{h_0\left(1-h_0\right)}{2+{\mathit{Ah}}_0^2\left(1-h_0\right)}\right)\frac{\mathit{dp}}{\mathit{dx}}.$$E72

Figure 38 compares the exact and approximate solutions (the solid line is the exact solution, and the dotted line is the approximate one obtained using (71) and (72) at $A=40000,h_0=0,2$).

Figure 38.

Solution comparison.

3. Application of the moving node method

In the first chapter, we considered MNM for some boundary value problems in order to obtain an approximate analytical solution. This chapter focuses on some uses of moved nodes. With the help of multipoint moving nodes, improved schemes are built for the convective-diffusion problem. It is proposed to improve the accuracy of schemes using the Richardson extrapolation method. Some properties of schemes are also presented for research with the help of MNM

3.1 Obtaining discrete compact schemes for the convective-diffusion problem of MNN

Of great interest is the construction and analysis of the discretization of a singularly perturbed ordinary differential equation of the second order. The equation of convection-diffusion is basic in modeling fluid flow at high Reynolds numbers and in convective mass exchange at high Peclet numbers. Many works are devoted to this subject [25, 40, 41, 42, 43, 44, 45].

Numerical solutions of the convection-diffusion equation often show numerical fluctuations. In practical calculations, many authors have observed parasitic oscillations at high Peclet numbers when the central approximation for the convective term is used. On the other hand, the upwind scheme usually leads to unpleasant artificial numerical diffusion.

This dilemma is central to the numerical solutions of convection-diffusion problems.

In order to compute approximate solutions to a partial differential equation, some form of local approximation must be used. This means that the decision values at each node are used to generate an approximate decision value. With finite differences, one usually tries to make the local area as compact as possible, for example, using only neighboring nodes when updating on a node.

If we consider the approximation of the convection-diffusion problem on a uniform grid, we can observe that most of the literature deals with the choice between schemes in a three-point pattern: (W,P,E). To obtain an approximation of a high order of accuracy, it is necessary to increase the number of points of the computational pattern.

Here we use the structure described in the first chapter to derive a new finite difference scheme [14]. Although such a procedure cannot be easily generalized to partial differential equations with variable coefficients.

Consider the DE of convection-diffusion

$$\frac{\mathit{dФ}}{\mathit{dx}}=\frac{1}{\mathit{Pe}}\frac{d^2Ф}{{\mathit{dx}}^2}+S\left(x\right),$$E73

with boundary conditions.

$$Ф\left(0\right)={Ф}_0,Ф\left(1\right)={Ф}_1$$E74

where Pe is the Peclet number $\left(\mathit{Pe}=\mathit{\rho vL}/Г\right)$, (v is the velocity, $\rho -$ density, L is the length scale, Г is the diffusion coefficient, x is the dimensionless coordinate, S(х) is the source.

On [0,1] we introduce a non-uniform grid

$$\mathrm{\Omega }=\left\{x_i,i=0,1,2,\dots ,N,0=x_0<{х}_1<\dots <x_{i-1}<x_i<x_{i+1}<\dots <x_N=1\right\}.$$

In the first chapter, with the help of moving nodes, an analytical solution to problem (73), (74) was constructed. When constructing compact circuits, we rely on a circuit against the flow, which is monotonic for any Peclet numbers.

Let us rewrite the scheme against the flow (21) for the segment (W,E)

$$\mathit{Pe}\frac{U^1-U_W^1}{x-W}=\frac{2}{\left(E-W\right)}\left(\frac{U_E^1-U^1}{E-x}-\frac{U^1-U_W^1}{x-W}\right)+\mathit{Pe}\cdot S\left(x\right).$$E75

In (75), the equation relates the unknown function at three points: $W,x,E,$ i.e. Eq. (75) is written in a three-point pattern. Now let us write Eq. (75) for an arbitrary internal node $x_i$, which is connected with neighboring nodes $x_{i-1},x_{i+1}$. Then

$$\mathit{Pe}\frac{U_i^1-U_{i-1}^1}{x_i-x_{i-1}}=\frac{2}{\left(x_{i+1}-x_{i-1}\right)}\left(\frac{U_{i+1}^1-U_i^1}{x_{i+1}-x_i}-\frac{U_i^1-U_{i-1}^1}{x_i-x_{i-1}}\right)+\mathit{Pe}\cdot S\left(x_i\right).$$E76

Here, $i=1,2,\dots ,i,\dots ,N-1$ and $U_i^1$ means the approximate value of the unknown function at the node $x_i$.

This schema can be rewritten like this:

$$a_P^1{U}_P^1=a_E^1{U}_E^1+a_W^1{U}_W^1+F_i^1,$$E77

But now

$$\begin{array}{l}{a}_E^1=\frac{2}{\left(x_{i+1}-x_{i-1}\right)\left(x_{i+1}-x_i\right)},a_W^1=\frac{\mathit{Pe}}{\left(x_i-x_{i-1}\right)}+\frac{2}{\left(x_{i+1}-x_{i-1}\right)\left(x_{i+1}-x_i\right)},\\ a_P^1=a_E^1+a_W^1,F_i^1=\mathit{Pe}\cdot S\left(x_i\right)\end{array}$$E78

To increase accuracy, based on three moving nodes (28),

$$a_P^3{U}_P^3=a_E^3{U}_E^3+a_W^3{U}_W^3+F_i^3,$$E79

here $a_E^3=\frac{8}{\left(x_{i+1}-x_{i-1}\right)\left(x_{i+1}-x_i\right)\left(1+{\gamma }_1\right)},a_W^3=\frac{2\mathit{Pe}}{\left(x_i-x_{i-1}\right)\left(1+{\tau }_1\right)}+\frac{8}{\left(x_{i+1}-x_{i-1}\right)\left(x_i-x_{i-1}\right)\left(1+{\tau }_1\right)},$ $a_P^3=a_W^3+a_E^3,$ $\theta =\mathit{Pe}\left(x_{i+1}-x_i\right),\sigma =\mathit{Pe}\left(x_i-x_{i-1}\right),$${\tau }_1=2/\left(2+\sigma \right),{\gamma }_1=\left(2+\theta \right)/2$ $F_i^3=\mathit{Pe}\cdot S\left(x_i\right)+\frac{4+\mathit{Pe}\cdot \left(x_{i+1}-x_{i-1}\right)}{x_{i+1}-x_{i-1}}\cdot \frac{1-{\tau }_1}{1+{\tau }_1}\cdot S\left(x_{i-1/2}\right)+\frac{4}{x_{i+1}-x_{i-1}}\cdot \frac{{\gamma }_1-1}{{\gamma }_1+1}\cdot S\left(x_{i+1/2}\right)$

$$x_{i-1/2}=0,5\left(x_{i-1}+x_i\right),x_{i+1/2}=0,5\left(x_i+x_{i+1}\right).$$

Based on with $2^k-1$ moving nodes (31), we have

$$a_P^{\left(2^k-1\right)}{U}_P^{\left(2^k-1\right)}=a_E^{\left(2^k-1\right)}{U}_E^{\left(2^k-1\right)}+a_W^{\left(2^k-1\right)}{U}_W^{\left(2^k-1\right)}+{F_i}^{\left(2^k-1\right)},$$E80

where

$$\begin{array}{l}{a}_E^{\left(2^k-1\right)}=\frac{2^{2k+1}\left(1-{\gamma }_k\right)}{\left(x_{i+1}-x_{i-1}\right)\left(x_{i+1}-x_i\right)\left(1-{\gamma }_k^{2^k}\right)},a_W^{\left(2^k-1\right)}=\frac{2^{2k+1}\mathit{Pe}\left(1-{\tau }_k\right)}{\left(x_i-x_{i-1}\right)\left(1-{\tau }_k^{2^k}\right)}+\\ \frac{2^{2k+1}\left(1-{\tau }_k\right)}{\left(x_{i+1}-x_{i-1}\right)\left(x_i-x_{i-1}\right)\left(1-{\tau }_k^{2^k}\right)},\end{array}$$

$$a_P^{\left(2^k-1\right)}=a_W^{\left(2^k-1\right)}+a_E^{\left(2^k-1\right)}.{\tau }_k=2^k/\left(2^k+\sigma \right),{\gamma }_k=\left(2^k+\theta \right)/2^k,$$

$${F_i}^{\left(2^k-1\right)}=\mathit{Pe}\cdot S\left(x_i\right)+\frac{2^{k+1}+\mathit{Pe}\cdot \left(x_{i+1}-x_{i-1}\right)}{x_{i+1}-x_{i-1}}\frac{{\left(1-{\tau }_k\right)}^2}{1-{\tau }_k^{2^k}}\sum _{j=1}^{2^k-1}\sum _{i=1}^j{\tau }_k^{i-1}\cdot S\left(x_i+j\frac{x_i-x_{i-1}}{2^k}\right)-$$

$$\frac{2^{k+1}}{x_{i+1}-x_{i-1}}\frac{{\left(1-{\gamma }_k\right)}^2}{1-{\gamma }_k^{2^k}}\sum _{j=1}^{2^k-1}\sum _{i=1}^j{\gamma }_k^{i-1}\cdot S\left(x_i+\left(2^k-j\right)\frac{x_{i+1}-x_i}{2^k}\right).$$

Let us consider numerical experiments.

Figures 39 and 40 show graphs for solving problem (73), (74) for $\mathit{Pe}=50$on segments $\left[0;1\right]$ with boundary conditions ${Ф}_0=0,{Ф}_1=1.$ Figure 39 corresponds to $S\left(x\right)=5cos4x,$ and the graphs in Figure 40 are obtained with

Figure 39.

Comparison of various schemes with source term $S\left(x\right)=5cos4x.$

Figure 40.

Comparison of various schemes with source term (81).

$$S\left(x\right)=\left\{\begin{array}{c}10-50x\mathit{\text{if}}х\le 0,3\\ 50x-20\mathit{\text{if}}0,3<x<0,4\\ 0\mathit{\text{if}}0,4<x\end{array}\right.$$E81

The numerical results are obtained for $h=0,1$ and the grid Peclet number is equal to 5. The solid lines are plots of the exact solutions of the problem. Circle symbols are obtained for the upwind scheme, rectangles according to the Patankar scheme, asterisks according to (79), diamonds according to (80) at $k=2$, and circles according to (80) at $k=7$.

Table 1 shows the root-mean-square errors $\sigma =\sqrt{{\sum }_1^N{\left(Ф\left(x_i\right)-U_i\right)}^2/N}$ for the considered schemes. Фi is exact solution at nodal points, Ui is numerical solution obtained by the considered schemes, N number of nodes.

Scheme	Upwind	Patankar	(79)	(80), k = 2	(80), k = 7
$S\left(x\right)=5cos4x$	0.361	0.300	0.192	0.096	0.004
$S\left(x\right)$ with (81)	0.282	0.199	0.154	0.077	0.003

Table 1.

The root-mean-square errors.

From Figures 39 and 40, and from the Table 1, it is clear that the proposed schemes give good results.

3.2 Construction of compact schemes of the convective-diffusion problem based on the finite volume method

The finite volume method is one of the methods that can give a good approximate solution to the problem. Here we explore the application of the finite volume method to solve the convection-diffusion equation for constructing compact schemes.

The basic strategy of all finite volume methods is to write the differential equation in a conservative form, integrate it over small domains (called “cells” or “finite volumes”), and transform each such integral over the cell boundary.

Our goal is to construct a qualitative scheme for the problem (82)

$$\frac{d}{\mathit{dx}}\left(\mathit{\rho uФ}\right)=\frac{d}{\mathit{dx}}\left(Г\frac{\mathit{dФ}}{\mathit{dx}}\right)+S\left(x\right)$$E82

$$Ф\left(0\right)={Ф}_0,Ф\left(1\right)={Ф}_1$$E83

based on the control volume method. The procedure for obtaining a scheme is similar to that described in paragraph 3.1.

On [0,1] we introduce a non-uniform grid

$$\mathrm{\Omega }=\left\{x_i,i=0,1,2,\dots ,N,0=x_0<{х}_1<\dots <x_{i-1}<x_i<x_{i+1}<\dots <x_N=1\right\}.$$

In the first chapter, with the help of moving nodes, an analytical solution to problem (82), (83) was constructed using the control volume method in the form.

$$\begin{array}{l}\left[\frac{\left(1-{\tau }_k\right){\beta }_k^{+}}{1-{\tau }_k^{2^k}}+\frac{\left(1-{\gamma }_k\right){\alpha }_k^{-}}{1-{\gamma }_k^{2^k}}\right]U^k=\frac{\left(1-{\tau }_k\right){\beta }_k^{+}}{1-{\tau }_k^{2^k}}{U}_W^k+\frac{\left(1-{\gamma }_k\right){\alpha }_k^{-}}{1-{\gamma }_k^{2^k}}{U}_E^k+\frac{E-W}{2^{k+1}}\cdot S\left(x\right)\end{array}$$

$$\begin{array}{l}+\frac{1-{\tau }_k}{1-{\tau }_k^{2^k}}\cdot \frac{x-W}{2^k}\cdot \sum _{j=1}^{2^k-1}\sum _{i=1}^j{\tau }_k^{i-1}S\left(W+j\frac{x-W}{2^k}\right)\\ +\frac{1-{\gamma }_k}{1-{\gamma }_k^{2^k}}\cdot \frac{E-x}{2^k}\cdot \sum _{j=1}^{2^k-1}\sum _{i=1}^j{\gamma }_k^{i-1}S\left(x+\left(2^k-j\right)\frac{E-x}{2}\right).\end{array}$$E84

Here ${\tau }_k=\frac{{\beta }_k^{-}}{{\beta }_k^{+}},{\gamma }_k=\frac{{\alpha }_k^{+}}{{\alpha }_k^{-}},{\beta }_k^{-}=2^k{D}_W+F^{-},{\beta }_k^{+}=2^k{D}_W+F^{+},$ ${\alpha }_k^{-}=2^k{D}_E+F^{-},$ ${\alpha }_k^{+}=2^k{D}_E+F^{+},$ $F^{-}=max\left(-F,0\right),$ $F^{+}=max\left(F,0\right),$ $D_E=Г/\left(E-x\right),$ $D_W=Г/\left(x-W\right)$.

Now let us write Eq. (84) for an arbitrary internal node $x_i$, which is connected with neighboring nodes $x_{i-1},x_{i+1}$.

Then

$$\begin{array}{l}\left[\frac{\left(1-{\tau }_k\right){\beta }_k^{+}}{1-{\tau }_k^{2^k}}+\frac{\left(1-{\gamma }_k\right){\alpha }_k^{-}}{1-{\gamma }_k^{2^k}}\right]{U_P}^k=\frac{\left(1-{\tau }_k\right){\beta }_k^{+}}{1-{\tau }_k^{2^k}}{U}_W^k+\frac{\left(1-{\gamma }_k\right){\alpha }_k^{-}}{1-{\gamma }_k^{2^k}}{U}_E^k+\frac{x_{i+1}-x_{i-1}}{2^{k+1}}\cdot S\left(x_i\right)+\\ \end{array}$$

$$\begin{array}{l}\frac{1-{\tau }_k}{1-{\tau }_k^{2^k}}\cdot \frac{x_i-x_{i-1}}{2^k}\cdot \sum _{j=1}^{2^k-1}\sum _{m=1}^j{\tau }_k^{m-1}S\left(x_{i-1}+j\frac{x_i-x_{i-1}}{2^k}\right)+\\ \frac{1-{\gamma }_k}{1-{\gamma }_k^{2^k}}\cdot \frac{x_{i+1}-x_i}{2^k}\cdot \sum _{j=1}^{2^k-1}\sum _{m=1}^j{\gamma }_k^{m-1}S\left(x_i+\left(2^k-j\right)\frac{x_{i+1}-x_i}{2}\right).\end{array}$$E85

What does it have to do with $D_E=Г/\left(x_{i+1}-x_i\right),D_W=Г/\left(x_i-x_{i-1}\right).$

3.3 Improving the accuracy of circuits using the Richardson extrapolation method

The Richardson extrapolation method is used to solve grid problems on a sequence of grids. The method consists in carrying out calculations for the same circuit, with different steps. Then we have several grid solutions. On the basis of the grid solutions, some linear combination is compiled. The resulting linear combination has a higher order of accuracy.

Creation of new schemes using Richardson extrapolation based on the schemes given in paragraph 3.1.

The accuracy of scheme (76), with a uniform arrangement of grid nodes, is $O\left(h\right)$. Scheme (79) has order $O\left(h/2\right)$. For a linear combination $Q^3\left(x_i\right)=-\frac{1}{3}{U}^1\left(x_i\right)+\frac{4}{3}{U}^3\left(x_i\right)$, we get an approximation error for a uniform grid $O\left(h^2\right)$. A linear combination of $U^1\left(x_i\right),U^3\left(x_i\right)$ and $U^7\left(x_i\right)$ in the form $Q^7\left(x_i\right)=\frac{1}{45}{U}^1\left(x_i\right)-\frac{4}{9}{U}^3\left(x_i\right)+\frac{64}{45}{U}^7\left(x_i\right)$ has an approximation order of $O\left(h^4\right)$. Consider $N=10$, $S\left(x\right)=x^2$, $\mathit{Pe}=30$. Table 2 shows the absolute difference between the exact and approximate solutions according to the schemes.

$x$	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9
$U^1\left(x_i\right)$	0.001	0.004	0.007	0.011	0.017	0.025	0.039	0.073	0.160
$U^3\left(x_i\right)$	0.001	0.002	0.004	0.006	0.008	0.012	0.018	0.034	0.089
$U^7\left(x_i\right)$	0.000	0.001	0.002	0.003	0.005	0.007	0.010	0.019	0.046
$Q^3\left(x_i\right)$	0.000	0.001	0.002	0.004	0.006	0.007	0.010	0.021	0.065
$Q^7\left(x_i\right)$	0.000	0.001	0.001	0.002	0.003	0.005	0.007	0.014	0.030

Table 2.

The absolute difference between the exact and approximate solutions.

Table 3 shows the root-mean-square error $\sigma =\sqrt{{\sum }_1^N{\left(Ф\left(x_i\right)-U_i\right)}^2/N}$ for the considered schemes. $Ф\left(x_i\right)$ the exact solution at the nodal points, $U_i$ is the numerical solution obtained by the considered schemes.

Schemes	$U^1\left(x\right)$	$U^3\left(x\right)$	$U^7\left(x\right)$	$Q^3\left(x\right)$	$Q^7\left(x\right)$
S = x^2, Pe = 50, ФW = 0, ФE = 1	0.047	0.023	0.011	0.015	0.006
S = 10, Pe = 50, ФW = 0, ФE = 1	0.033	0.017	0.008	0.011	0.005
S = x^2, Pe = 100, ФW = 0, ФE = 1	0.034	0.014	0.006	0.008	0.003
S = 5cos(4πx), Pe = 50, ФW = 0, ФE = 1	0.213	0.120	0.061	0.090	0.038

Table 3.

Comparison by the root-mean-square errors.

Figures 41 and 42 show numerical solutions for Ф_W = 0, Ф_E = 0.

Figure 41.

Pe = 100, S = 5cos4πx Solid curve exact solution, circle obtained by scheme U¹, circle by U³, solid rectangle by U⁷, diamond by Q³, star by Q⁷.

Figure 42.

Pe = 100, S = exp.(−4x). Solid curve exact solution, circle obtained by scheme U¹, circle by U³, solid rectangle by U⁷, diamond by Q³, star by Q⁷.

From the graphs in Figures 41 and 42, and from Tables 2 and 3, it is clear that the Richardson linear combination allows you to get a more improved circuit.

3.4 Influence of the choice of profile on the face of the control volume on the quality of difference schemes

When obtaining discrete analogs of the convective-diffusion problems given above, on the basis of multipoint PUs, it was possible to construct better compact circuits in a three-point template. However, there is another approach to improve the quality of the scheme based on the choice of the decision profile.

Since the work of Leonard, in order to improve the results of the numerical solution, attempts have been made to improve the algorithm, which is built in a five-point pattern.

In all the above schemes (except for the scheme against the flow), the conditions of boundedness and non-negativity of the coefficients are violated.

Here it is proposed to improve the scheme based on the choice of the solution profile on the edge of the control volume in the three-point template of the convective-diffusion problem. The upwind scheme with one-sided differences is taken as the initial scheme. The QUICK scheme uses quadratic upwind interpolation to determine the convective flow. Here we use the solution obtained by the upwind scheme based on the method of moving nodes.

MNM for simple cases allows one to obtain an analytical representation of the solution between the nodal points of the boundary value problem. Based on this representation, it is possible to construct a better discrete scheme.

We integrate (73) over the control volume $\left[w,e\right]$

$${Ф}_e-{Ф}_w=\frac{1}{\mathit{Pe}}{\left(\frac{\mathit{dФ}}{\mathit{dx}}\right)}_e-\frac{1}{\mathit{Pe}}{\left(\frac{\mathit{dФ}}{\mathit{dx}}\right)}_w+\underset{w}{\overset{e}{\int }}S\left(x\right)\mathit{dx}.$$

Replacing the derivatives with difference relations, we have

$${Ф}_e-{Ф}_w=\frac{1}{\mathit{Pe}}\frac{{Ф}_E-{Ф}_P}{x_E-x_P}-\frac{1}{\mathit{Pe}}\frac{{Ф}_P-{Ф}_W}{x_P-x_W}+\left(x_e-x_w\right)f_P.$$E86

Here $f_P=\frac{1}{x_e-x_w}\underset{w}{\overset{e}{\int }}S\left(x\right)\mathit{dx}$. Depending on the type of function profile $Ф$ on the control volume, different schemes are obtained.

Let the profile $Ф$ be piecewise constant in each control volume. Then, assuming ${Ф}_e={Ф}_P,{Ф}_w={Ф}_W,$ we have an upwind scheme:

$${Ф}_P-{Ф}_W=\frac{1}{\mathit{Pe}}\frac{{Ф}_E-{Ф}_P}{x_E-x_P}-\frac{1}{\mathit{Pe}}\frac{{Ф}_P-{Ф}_W}{x_P-x_W}+\left(x_e-x_w\right)f_P.$$E87

If the profile $Ф$ is linear between the nodes and the edges of the control volume are located in the middle between the node points, we have a scheme with central differences:

$$\frac{{Ф}_E+{Ф}_P}{2}-\frac{{Ф}_P+{Ф}_W}{2}=\frac{1}{\mathit{Pe}}\frac{{Ф}_E-{Ф}_P}{x_E-x_P}-\frac{1}{\mathit{Pe}}\frac{{Ф}_P-{Ф}_W}{x_P-x_W}+\left(x_e-x_w\right)f_P.$$E88

To improve the accuracy of circuits, many authors recommended various circuits. All these schemes are multipoint (more than three). Here is a way to improve three-point circuits.

From (87) we get

$$\begin{array}{l}{Ф}_P=\frac{x_P-x_W}{P\mathrm{e}\left(x_E-x_P\right)\left(x_P-x_W\right)+x_E-x_W}{Ф}_E+\frac{\left(x_E-x_P\right)\left(1+P\mathrm{e}\left(x_P-x_W\right)\right)}{P\mathrm{e}\left(x_E-x_P\right)\left(x_P-x_W\right)+x_E-x_W}{Ф}_W+\\ =\frac{\left(x_P-x_W\right){Ф}_E+\left(x_E-x_P\right)\left(1+\mathit{Pe}\left(x_P-x_W\right)\right){Ф}_W}{P\mathrm{e}\left(x_E-x_P\right)\left(x_P-x_W\right)+x_E-x_W}\end{array}$$E89

If the nodes $x_E$ and $x_W$ are fixed, and the node $x_P$ is movable, we get a profile ${Ф}_P$ between the nodes $x_E$ and $x_W$. This profile is used in (86) to determine ${Ф}_e$ and ${Ф}_w$.

To improve scheme (87), we proceed as follows. Eq. (89) connects at three nodes ($x_W,x_P$, $x_E$), if we apply Eq. (89) for nodes ( $x_W,x_w$, $x_P$), we have

$${Ф}_w=\frac{2+R_h}{4+R_h}{Ф}_P+\frac{2}{4+R_h}{Ф}_E+\frac{R_h}{2\left(R_h+4\right)}\frac{h}{4}{f}_w.$$E90

Similarly, for nodes ( $x_P,x_e$, $x_E$), we have

$${Ф}_e=\frac{2+R_h}{4+R_h}{Ф}_P+\frac{2}{4+R_h}{Ф}_E+\frac{R_{h}h}{2\left(R_h+4\right)}{f}_e.$$E91

Substituting (90) and (91) into (86) we have

$$\left[\frac{R_h^2}{4+R_h}+2\right]{Ф}_P=\left[1+\frac{2+R_h}{4+R_h}\right]{Ф}_W+\left[1-\frac{2R_h}{4+R_h}\right]{Ф}_E+h{R}_h{f}_P-\frac{h\cdot R_h^2}{2\left(4+R_h\right)}\left(f_e-f_w\right).$$E92

The condition $R_h<4$ is ensured by the positivity of the coefficients and the stability of the scheme (92).

Proceeding similarly as in the derivation of (92), but using (92) for the profile, for a uniform step we obtain

$$\begin{array}{l}\left[\frac{{\left(4+R_h\right)}^2-16}{{\left(4+R_h\right)}^2+16}\right]{Ф}_P=\left[\frac{1}{R_h}-\frac{16}{{\left(4+R_h\right)}^2+16}\right]{Ф}_E+\left[\frac{{\left(4+R_h\right)}^2}{{\left(4+R_h\right)}^2+16}+\frac{1}{R_h}\right]{Ф}_W+\\ \mathit{hS}\left(P\right)+\frac{h\left(8R_h+8R_h^2\right)}{2\left[{\left(4+R_h\right)}^2+16\right]}\cdot \left(S\left(x_w\right)-S\left(x_e\right)\right),\end{array}$$E93

Test problems

1. Consider the equation

$$\frac{\mathit{du}}{\mathit{dx}}=\frac{1}{P\mathrm{e}}\frac{d^{2}u}{{\mathit{dx}}^2}+sin\mathit{\pi x}.$$

with boundary conditions $u\left(0\right)=u\left(1\right)=0$. Table 4 shows the maximum absolute differences of the schemes calculated at the nodal points ( $u$ is the exact solution of the problem, $u_1$ is the solution obtained according to the upwind scheme, $u_2$ is according to the power law, $u_3$ according to the Leonard scheme, $u_4$ according to (92) and $u_5$- according to the scheme (93).

$P\mathrm{e}$	$R_h$	$max\left|u-u_1\right|$	$max\left|u-u_2\right|$	$max\left|u-u_3\right|$	$max\left|u-u_4\right|$	$max\left|u-u_5\right|$
100	10	0.0526	0.03770	0.1801	0.03701	0.00077
1000	100	0.0470	0.0464	0.2732	0.01607	0.00927

Table 4.

The maximum absolute differences.

2. Consider the equation

$$\frac{\mathit{du}}{\mathit{dx}}=\frac{1}{P\mathrm{e}}\frac{d^{2}u}{{\mathit{dx}}^2}+\mathrm{s}\left(x\right),$$

with boundary conditions $u\left(0\right)=0,u\left(1\right)=1$, with source

$$s\left(x\right)=\left\{\begin{array}{c}10-50x,0\le x\le 0.3,\\ 50x-20,0.3<x\le 0.4\\ 0,0.4<x\le 1\end{array}\right.,$$

Figure 43 shows that scheme (93) gives the best results. Leonard’s scheme gives an incorrect solution near the right boundary. Scheme (92) also exhibits a slight non-monotonicity. This is due to the fact that scheme (92) is stable for R_k<4.

Figure 43.

Comparison of various schemes. $P\mathrm{e}=100,R_h=5.$ The solid line is the exact solution, the circle is the upwind scheme, the circle is the Patankar scheme, the asterisk is the Leonard scheme, + is the scheme (92), the diamond is according to (93).

Figure 44 Shows that for large grid Peclet numbers, the upstream and Patankar schemes give close results. Scheme (93) gives the best results. This can also be seen in Table 5, which compares the considered schemes (SDS—central scheme).

Figure 44.

Comparison of various schemes. $Re=500,R_h=25.$ The solid line is the exact solution, the circle is the upwind scheme, the circle is the Patankar scheme, the asterisk is the Leonard scheme, + is the scheme (92), the diamond is according to (93).

Scheme	$Re$	$h$	$R_h$	$max\left|u-u_p\right|$	$\frac{\sum \left|u_i-{\left(u_p\right)}_i\right|}{\sum u_i}$
Upwind	100	1/40	2.5	0.1627	0.2116
	100	1/20	5	0.3258	0.4224
	500	1/20	25	0.3650	0.4101
Power	100	1/40	2.5	0.0833	0.1057
	100	1/20	5	0.2385	0.3025
	500	1/20	25	0.3454	0.3868
(8)	100	1/40	2.5	0.0164	0.0169
	100	1/20	5	0.0460	0.0401
	500	1/20	25	0.0531	0.0398
(12)	100	1/40	2.5	0.0129	0.00840
	100	1/20	5	0.0452	0.0358
	500	1/20	25	0.0571	0.0404
QUICK	100	1/40	2.5	0.0700	0.0701
	100	1/20	5	0.2231	0.1931
	500	1/20	25	0.3653	0.2055
CDS	100	1/40	2.5	0.1237	0.0062
	100	1/20	5	0.3033	0.0467
	500	1/20	25	0.5136	0.1355

Table 5.

Comparison of circuits with respect to grid Peclet number.

3. Two-dimensional case. Consider the equation

$$\frac{\partial g}{\partial x}=\frac{1}{Re}\left(\frac{{\partial }^{2}g}{\partial x^2}+\frac{{\partial }^{2}g}{\partial y^2}\right)+s\left(x,y\right).$$

Exact solution $g=6y^{10}\left(1-y^{10}\right)\left(1-x^3\right)+6x^{3}y\left(1-y\right)$. The equations are solved in the area $\left[0,1\right]\times \left[0,1\right]$. The source term is defined so that the given function is a solution to the equation. The boundary conditions were determined based on the exact solution. Table 6 shows the results of calculations according to the schemes.

Scheme	$Re=100,n=5,h=0,1$		$Re=500,n=5,h=0,1$		$Re=1000,n=10,h=0,1$
	$max\left|g-g_p\right|$	$\frac{\sum \left|g-g_p\right|}{\sum \left|g\right|}$	$max\left|g-g_p\right|$	$\frac{\sum \left|g-g_p\right|}{\sum \left|g\right|}$	$max\left|g-g_p\right|$	$\frac{\sum \left|g-g_p\right|}{\sum \left|g\right|}$
Upwind	0.150	0.074	0.169	0.074	0.186	0.129
CDS	0.074	0.023	0.035	0.018	0.470	0.382
Power	0.130	0.061	0.165	0.071	0.184	0.127
QUICK	0.063	0.017	0.020	0.0051	0.097	0.016
(8)	0.035	0.019	0.013	0.008	0.057	0.023
(12)	0.033	0.016	0.008	0.005	0.060	0.015
VONOS	0.055	0.016	0.019	0.005	0.073	0.015

Table 6.

Results of calculations of errors according to the schemes.

From Table 6, it is clear that the proposed schemes show the best results.

3.5 Schema improvement with flow equality

MNM can improve the quality of the scheme. We demonstrate this method based on the upwind scheme (87) written in the form:

$$\frac{{Ф}_P-{Ф}_W}{x_P-x_W}=\frac{2}{\mathit{Pe}\left(x_E-x_W\right)}\left(\frac{{Ф}_E-{Ф}_P}{x_E-x_P}-\frac{Ф-{Ф}_W}{x_P-x_W}\right)+S\left(x_P\right).$$E94

In (94) we pass to the limit at $x_E\to x_P$ and, assuming the existence of the limit, we have

$$\frac{{Ф}_P-{Ф}_W}{x_P-x_W}=\frac{2}{\mathit{Pe}\left(x_P-x_W\right)}\left(\frac{{\mathit{dФ}}_P^{-}}{{\mathit{dx}}_P}-\frac{{Ф}_P-{Ф}_W}{x_P-x_W}\right)+S\left(x_P\right).$$

Here, ${\mathit{dФ}}_P^{-}/{\mathit{dx}}_P$ is the left-hand derivative of the unknown function at the point $x_P$. From here

$$\frac{{\mathit{dФ}}_P^{-}}{{\mathit{dx}}_P}=\frac{2+\mathit{Pe}\left(x_P-x_W\right)}{2}\cdot \frac{{Ф}_P-{Ф}_W}{x_P-x_W}-\frac{\mathit{Pe}\left(x_P-x_W\right)}{2}\cdot S\left(x_P\right),$$E95

Similarly, taking an arbitrary point $x\in \left(x_P,x_E\right)$ and passing to the limit $x\to x_P$, we find

$$\frac{{\mathit{dФ}}_P^{+}}{{\mathit{dx}}_P}=\frac{2}{2+\mathit{Pe}\left(x_E-x_P\right)}\cdot \frac{{Ф}_E-{Ф}_P}{x_E-x_P}+\frac{\mathit{Pe}\left(x_E-x_P\right)}{2+\mathit{Pe}\left(x_E-x_P\right)}\cdot S\left(x_P\right),$$

By equating ${\mathit{dФ}}^{+}/\mathit{dx}={\mathit{dФ}}^{-}/\mathit{dx}$ the flows, we get an improved scheme:

$$c_P{Ф}_P=a_P{Ф}_W+b_P{Ф}_E+d_{P}S\left(x_P\right)$$E96

where

$a_P=\frac{2+\mathit{Pe}\left(x_P-x_W\right)}{\left(x_P-x_W\right)},$$b_P=\frac{2}{\left[2+\mathit{Pe}\left(x_E-x_P\right)\right]\left(x_E-x_P\right)},c_P=a_P+b_P,$.

Figure 45 shows a comparison of the exact solution and the schemes according to (87) and (96) for Pe = 5, with one moving node (S(х) = 0). It can be seen from the graph that the solution is improving. Numerical diffusion decreases.

Figure 45.

Comparison of schemes. The solid curve is the exact solution, the dotted line according to (87), the dotted line according to (96).

3.6 Investigation of the scheme by the MNM

At this point, we are dealing with monotonicity and MMN approximation of the circuit. On the basis of the analytical form of the approximate solution of the problem between the nodes, which is obtained on the basis of the MMN, it is possible to investigate monotonicity and the type of approximation of the scheme.

3.6.1 Investagation of monotonicity

Scheme with central-difference approximation of the convective term. Consider Eq. (73). Take a segment $\left[x_{i-1},x_{i+1}\right]\subset \left[0,1\right]$ and any point $x\equiv x_i\in \left(x_{i-1},x_{i+1}\right)$. Consider the grid analog (73)

$$\frac{u_{i+1}-u_{i-1}}{x_{i+1}-x_{i-1}}=\frac{2}{\mathit{Pe}\left(x_{i+1}-x_{i-1}\right)}\left(\frac{u_{i+1}-u}{x_{i+1}-x}-\frac{u-u_{i-1}}{x-x_{i-1}}\right)+S\left(x\right)$$E97

If we set $x=\left(x_{i+1}+x_{i-1}\right)/2$, we have a central-difference approximation. Here, $u_{i+1}$ is the approximate value of the solution at the point $x_{i+1}$, $u$ is the approximate value of the solution at the point $x.$ To obtain a physically plausible solution in simple cases, we set S(x) = 0.

From (97) we find

$$u=\frac{\left(x-x_{i-1}\right)\left(2-\mathit{Pe}\left(x_{i+1}-x\right)\right)u_{i+1}+\left(x_{i+1}-x\right)\left(2+\mathit{Pe}\left(x-x_{i-1}\right)\right)u_{i-1}}{2\left(x_{i+1}-x_{i-1}\right)}.$$E98

By changing $x$ the values on the interval $\left(x_{i-1},x_{i+1}\right)$, we can determine the behavior of the solution. For given values $x_{i+1},x_{i-1},u_{i-1},u_{i+1}$ (98) is a parabola.

From (98) one can get

$$\frac{u-u_{i-1}}{u_{i+1}-u_{i-1}}=\frac{\left(x-x_{i-1}\right)\left(2-\mathit{Pe}\left(x_{i+1}-x\right)\right)}{2\left(x_{i+1}-x_{i-1}\right)}.$$E99

A physically plausible solution is obtained if $0\le \frac{u-u_{i-1}}{u_{i+1}-u_{i-1}}\le 1$. This condition imposes a restriction $2-\mathit{Pe}\left(x_{i+1}-x\right)\ge 0$. This condition is the condition of monotonicity of the central-difference scheme for a non-uniform grid. In the case of a uniform grid, we have $2\ge \mathit{Pe}\cdot h$. This condition is the well-known monotonicity condition [46]. For a coarse grid (N = 2, one movable node) at Pe = 5, the solution of exact and approximate solutions are shown in Figure 46.

Figure 46.

Comparison of solutions in a coarse grid. The dotted curve is approximate, the solid curve is exact, Pe = 5 (Ф₀ = 0, Ф₁ = 1).

In Figure 46, the solid curve represents the exact solution, while the dotted one represents the approximate solution obtained on the basis of (99). It can be seen from the graph that scheme (99) does not give a physically plausible analytical solution. That is why scheme (99) for large Peclet numbers gives an oscillatory numerical solution. A plausible solution should have the same qualitative character as the exact solution. When solving numerically, scheme (97) is implemented using a sweep, and for the stability of the sweep, the nodes are selected so that $2-\mathit{Pe}\left(x_{i+1}-x_i\right)\ge 0$. For example, for a coarse grid (one nodal point), the credibility condition gives $x_i\ge 0,6$. Indeed, for Pe = 5, it $2-\mathit{Pe}\left(1-х\right)\ge 0$ follows that $x_i\ge 0,6$ (see Figure 46).

For Pe = 2, comparisons of the solutions are shown in Figure 47, which gives a physically plausible solution.

Figure 47.

Comparison of the solution in a coarse grid. The dotted curve is approximate, the solid curve is exact, Pe = 2 (Ф₀ = 0, Ф₁ = 1).

Upwind scheme. Let us consider a difference analog of Eq. (73), in which the convective term is approximated by a one-sided difference relation (without a source)

$$\frac{u-u_{i-1}}{x-x_{i-1}}=\frac{2}{\mathit{Pe}\left(x_{i+1}-x_{i-1}\right)}\left(\frac{u_{i+1}-u}{x_{i+1}-x}-\frac{u-u_{i-1}}{x-x_{i-1}}\right).$$

From here we get

$$u=\frac{2\left(x-x_{i-1}\right)u_{i+1}+\left(x_{i+1}-x\right)\left(2+\mathit{Pe}\left(x_{i+1}-x_{i-1}\right)\right)u_{i-1}}{\left(x_{i+1}-x_{i-1}\right)\left(2+\mathit{Pe}\left(x_{i+1}-x\right)\right)}$$

or

$$\frac{u-u_{i-1}}{u_{i+1}-u_{i-1}}=\frac{2\left(x-x_{i-1}\right)}{\left(x_{i+1}-x_{i-1}\right)\left(2+\mathit{Pe}\left(x_{i+1}-x\right)\right)}.$$E100

Since, the right side of relation (100) into segments is a hyperbola and therefore we have $0\le \frac{u-u_{i-1}}{u_{i+1}-u_{i-1}}\le 1$. Those the upstream circuit is always monotonic. Figure 48 shows a comparison of the exact and approximate analytical solutions (Pe = 5). However, numerical diffusion occurs.

Figure 48.

Comparison of the solution in a coarse grid. The dotted curve is approximate, and the solid curve is exact (Ф₀ = 0, Ф₁ = 1).

3.7 An explicit expression of the approximation error of ordinary differential equations based on the moved node method

Here discusses the issue of the possibility of calculating the approximation error. When replacing differential equations with discrete ones, one of the key issues is the closeness of the discrete solution to the exact solution. For the difference solution to the problem, a grid area is formed. The discrete solution is determined at the nodal points. Traditionally, in questions of replacing a differential equation with a descriptive one, one usually indicates the degree of approximation of the O(h^p) type. Here h is the grid step.

However, it is possible to calculate the approximation error at nodal points based on the method of moving nodes. The method of moving nodes allows for obtaining an approximate analytical expression. On the basis of the approximate form, it is possible to calculate the approximation error. On the other hand, at each node one can construct a differential analog of the difference equation. Using simple examples, the calculation of approximation errors is demonstrated and schemes of the collocation type are constructed.

3.7.1 Introduction

Here describes the application of the moving nodes method to the calculation of the approximation error. When a two-point boundary value problem is solved by different methods, the question of the degree of approximation usually appears. The closeness of the exact and approximation of the solution, and the quality of the difference scheme are evaluated based on the degree of this parameter. With such an analysis, other parameters (the coefficients of the differential equation) are not explicitly involved in the approximation error expression. Obtaining an explicit expression for the approximation error makes it possible to analyze it.

Consider the simplest ordinary differential equation with boundary conditions

$$\frac{d^{2}u}{{\mathit{dx}}^2}=C,u\left(0\right)=0,u\left(1\right)=1$$E101

where С—const.

Create a uniform grid on segments [0,1] with step h. A uniform grid on a segment x ∈ [0, 1] with step h has the form:

$${\overline{\omega }}_h=\left\{x_k=\mathit{hk},k=0,1,\dots ,N,h\cdot N=1\right\}$$

Let us replace the second-order derivative with the difference relation:

$$\frac{U_{i+1}-2U_i+U_{i-1}}{h^2}=C,1\le i\le N-1,U_0=0,U_N=1$$E102

Difference scheme (102) traditionally has order O(h²). However, if we solve system (102) by the Tomas algorithm, we obtain a numerical solution that coincides with the exact analytical solution for any grid steps h at the grid nodes. Those. scheme (102) approximates (101) exactly.

3.7.2 Methodology

Let us have a differential equation

$$\mathit{Lu}=f,$$E103

where L is a differential operator, f is a known function, and u is an unknown function. (103) the equation is considered in some domain D with appropriate boundary conditions. The differential Eq. (103) is replaced by the difference equation:

$$L_h{u}_h=f_h,$$E104

where L_h is the difference operator, u_h is the unknown grid function, and f_h is the approximation of the function f at the grid nodes.

Usually, the approximation error is given as [2, 3]:

$$Q_h=L_h{\left[u\right]}_h-f_h,$$E105

where [u]_h is the exact solution of (103) at the grid nodes. Using the Taylor series, from (105) one obtains that, Q_h = O(h^m), where h is the grid step and m is the degree of approximation.

You can determine an explicit approximation error if you use the method of a moving node, which allows you to extend the definition to the entire area D. This allows you to introduce an approximation error like this:

$$R_h=L_h{\left\{u\right\}}_h-f_h.$$E106

Here {u}_h is a predefined continuous function by means of a moveable node. The approximate calculation of the approximation error of type (106) is demonstrated using simple examples.

3.7.3 Results and discussion

As an application of the above approach, consider examples.

1. Consider a simple boundary value problem:

$$\frac{d^{2}u}{{\mathit{dx}}^2}=f\left(x\right),u\left(0\right)=u_a,u\left(1\right)=u_b$$E107

Let us build a non-uniform grid on segments [0; 1]:

$${\overline{\omega }}_h=\left\{0=x_0,<x_1<\dots <x_{N-1}<x_N=1,k=0,1,\dots ,N\right\}$$

In the non-uniform grid, we replace (107) with the difference problem:

$$\frac{2}{x_{i+1}-x_{i-1}}\left(\frac{U_{i+1}-U_i}{x_{i+1}-x_i}-\frac{U_i-U_{i-1}}{x_i-x_{i-1}}\right)=f\left(x_i\right),i=1,2,\dots ,N-1.$$E108

Here $U_i$ is the grid solution of the problem. From here

$$U_i=\frac{U_{i+1}\left(x_i-x_{i-1}\right)+U_{i-1}\left(x_{i+1}-x_i\right)}{x_{i+1}-x_{i-1}}-\frac{1}{2}f\left(x_i\right)\left(x_i-x_{i-1}\right)\left(x_{i+1}-x_i\right),i=1,2,\dots ,N-1.$$E109

We redefine the value of the function at non-nodal points as follows. To do this, we consider in (109) $x_{i+1},x_{i-1},U_{i-1},U_{i+1}$, to be fixed, and $x_i$ to be moved, and the function f(x) to be smooth. Thus, we will complete the grid function on each segment $\left(x_{i-1},x_{i+1}\right)$. From (109) we get

$$U_i^{″}\left(x_i\right)=-\frac{1}{2}{f}^{″}\left(x_i\right)\left(x_{i+1}-x_i\right)\left(x_i-x_{i-1}\right)-f^{\prime }\left(x_i\right)\left(x_{i+1}+x_{i-1}-2x_i\right)+f\left(x_i\right)$$E110

Then the approximation error for the nodal points looks like this:

$$R_h\left(x_i\right)=-\frac{1}{2}{f}^{″}\left(x_i\right)\left(x_{i+1}-x_i\right)\left(x_i-x_{i-1}\right)-f^{\prime }\left(x_i\right)\left(x_{i+1}+x_{i-1}-2x_i\right)$$E111

If the grid is uniform for the approximation error, we obtain the expression

$$R_h\left(x_i\right)=-\frac{1}{2}{f}^{″}\left(x_i\right)h^2,i=1,2,\dots ,N-1.$$E112

If on the segments $\left(x_{i-1},x_{i+1}\right)$ the function constant approximation error is identically equal to zero and we get the exact solution.

Based on expression (110), the following conclusion can be drawn.

Given a two-point boundary value problem

$$\frac{d^{2}u}{{\mathit{dx}}^2}=f^{\ast }\left(x\right),u\left(0\right)=u_a,u\left(1\right)=u_b$$

and $f^{\ast }\left(x\right)$ can be represented as

$$f^{\ast }\left(x_i\right)=-\frac{1}{2}{f}^{″}\left(x_i\right)\left(x_{i+1}-x_i\right)\left(x_i-x_{i-1}\right)-f^{\prime }\left(x_i\right)\left(x_{i+1}+x_{i-1}-2x_i\right)+f\left(x_i\right)$$

then the difference scheme

$$\frac{2}{x_{i+1}-x_{i-1}}\left(\frac{U_{i+1}-U_i}{x_{i+1}-x_i}-\frac{U_i-U_{i-1}}{x_i-x_{i-1}}\right)=f\left(x_i\right),i=1,2,\dots ,N-1,$$

gives a grid solution coinciding with the exact solution at the nodal points.

If there is only one internal node point (the node being moved is one), then an approximate analytical solution can be obtained. Indeed, if we rewrite scheme (108) for one moving node, we have

$$2\left(\frac{U_i-U_{i-1}}{1-x}-\frac{U\left(x\right)-U_a}{x}\right)=f\left(x_i\right).$$E113

From here we obtain an approximate analytical solution:

$$U\left(x\right)=\frac{U_{b}x+U_a\left(1-x\right)}{x_{i+1}-x_{i-1}}-\frac{1}{2}f\left(x_i\right)\left(1-x\right)x.$$E114

In this case, (114) represents the exact solution to the problem (107).

if we put

$$f^{\ast }\left(x\right)=-\frac{1}{2}{f}^{″}\left(x\right)\left(1-x\right)x-f^{\prime }\left(x\right)\left(1-2x\right)+f\left(x\right).$$

The form of the approximation error (111) allows the construction of new schemes of the collocation type. Indeed, if in problem (108) we replace the right side with the expression

$$f\left(x_i\right)+A\left(x_i-x_{i-1}\right)\left(x_{i+1}-x_i\right),$$

Here A is still an unknown constant. Parameter A is determined so that the approximation error (111) for a uniform step at node $x_i$ is equal to zero, i.e. collocation type scheme. Then we have

$$A=\frac{1}{4}{f}^{″}\left(x_i\right)$$

2. Consider a stationary equation in which only convection and diffusion are present without a source.

$$\epsilon v^{″}+v^{\prime }=0,$$E115

with boundary conditions $v\left(0\right)=0,v\left(1\right)=1.$

There are various schemes for the difference solution (115). Based on the moving node technique, it is possible to explicitly express local errors in the approximation of differential equations. Using the moving node method, we will show the efficient calculation of local approximation errors for the model problem (115).

Scheme with central-difference approximation of the convective term. Take a segment $\left(x_{i-1},x_{i+1}\right)$ and any point $x\in \left(x_{i-1},x_{i+1}\right)$. Consider the different analog (115).

$$\frac{2\epsilon }{x_{i+1}-x_{i-1}}\left(\frac{u_{i+1}-u}{x_{i+1}-x}-\frac{u-u_{i-1}}{x-x_{i-1}}\right)+\frac{u_{i+1}-u_{i-1}}{x_{i+1}-x_{i-1}}=0$$E116

At $x=\left(x_{i-1}+x_{i+1}\right)/2$, we have a central difference approximation. Here, $u$ is the approximate value of the solution at point x.

From (116) we find.

$$u=\frac{\left(x-x_{i-1}\right)\left(2\epsilon +x_{i+1}-x\right)u_{i+1}+\left(x_{i+1}-x\right)\left(2\epsilon -x+x_{i-1}\right)u_{i-1}}{2\epsilon \left(x_{i+1}-x_{i-1}\right)}.$$E117

From here we get,

$$u^{\prime }=\frac{2\epsilon +x_{i+1}+x_{i-1}-2x}{2\epsilon }\frac{u_{i+1}-u_{i-1}}{x_{i+1}-x_{i-1}},$$E118

$$u^{″}=-\frac{1}{\epsilon }\frac{u_{i+1}-u_{i-1}}{x_{i+1}-x_{i-1}}.$$E119

If the difference solution at nodal points is known, then formula (117) makes it possible to determine the unknown at points that are not nodal.

Using formulas (118) and (119), the derivatives are restored at any point of the segment. Multiplying (119) by and adding with (118), we obtain.

$$\epsilon u^{″}+u^{\prime }={\mathrm{\Psi }}_1,$$E120

where

$${\mathrm{\Psi }}_1=\frac{x_{i+1}+x_{i-1}-2x}{2\epsilon }\frac{u_{i+1}-u_{i-1}}{x_{i+1}-x_{i-1}}.$$

Eq. (120) can be called a differential analog of the difference Eq. (16); difference Eq. (116) is a collocation-type scheme.

Using (119), the approximation error can be written as.

$${\mathrm{\Psi }}_1=-\frac{x_{i+1}+x_{i-1}-2x}{2}{u}^{″}.$$

Then Eq. (120) takes the form

$$\left(\epsilon +\frac{x_{i+1}+x_{i-1}-2x}{2}\right)u^{″}+u^{\prime }=0.$$E121

Thus, difference Eq. (116) exactly approximates differential Eq. (121) on the segment $\left[x_{i-1},x_{i+1}\right].$

Comparison of Eqs. (115) and (121) shows that when Eq. (115) is approximated by scheme (116), scheme diffusion appears with a variable coefficient $\left(x_{i+1}+x_{i-1}-2x\right)/2.$

Upwind Scheme. Let us consider the difference analog of Eq. (115), in which the convective term is approximated by the one-sided difference relation.

$$\frac{2\epsilon }{x_{i+1}-x_{i-1}}\left(\frac{u_{i+1}-u}{x_{i+1}-x}-\frac{u-u_{i-1}}{x-x_{i-1}}\right)+\frac{u_{i+1}-u}{x_{i+1}-x}=0.$$E122

From here we get

$$u=\frac{\left(x-x_{i-1}\right)\left(2\epsilon +x_{i+1}-x_{i-1}\right)}{\left(x_{i+1}-x_{i-1}\right)\left(2\epsilon +x-x_{i-1}\right)}\frac{u_{i+1}+2\epsilon \left(x_{i+1}-x\right)u_{i-1}}{}$$E123

Determine the first and second derivatives:

$$u^{\prime }=\frac{2\epsilon \left(2\epsilon +x_{i+1}-x_{i-1}\right)}{{\left(2\epsilon +x-x_{i-1}\right)}^2}\frac{u_{i+1}-u_{i-1}}{x_{i+1}-x_{i-1}},$$E124

$$u^{″}=\frac{-4\epsilon \left(2\epsilon +x_{i+1}-x_{i-1}\right)}{{\left(2\epsilon +x-x_{i-1}\right)}^3}\frac{u_{i+1}-u_{i-1}}{x_{i+1}-x_{i-1}}$$E125

Let us calculate the approximation error.

$${\mathrm{\Psi }}_2=\frac{2\epsilon \left(x-x_{i-1}\right)\left(2\epsilon +x_{i+1}-x_{i-1}\right)}{{\left(2\epsilon +x-x_{i-1}\right)}^3}\frac{u_{i+1}-u_{i-1}}{x_{i+1}-x_{i-1}}$$

The differential analog of scheme (122) has the form.

$$\left(\epsilon +\frac{x-x_{i-1}}{2}\right)u^{″}+u^{\prime }=0,$$E126

those with a scheme against the flow, we have a scheme diffusion with a coefficient $\left(x_{i+1}-x\right)/2$. Based on (123)—is a hyperbola, which is monotone on the segment, i.e. scheme (122) is monotonic.

Based on the form of the differential analog (126), we can conclude that the differential equation

$$\left(\epsilon +\frac{x}{2}\right)u^{″}+u^{\prime }=0$$E127

is exactly approximated by the scheme

$$2\epsilon \left(\frac{u_b-u}{1-x}+\frac{u-u_a}{x}\right)+\frac{u_b-u}{1-x}=0$$E128

Thus solving (128) with respect to u, we obtain the exact solution of differential Eq. (127).

3.8 On convergence of MNM

Let us show the convergence of MNM on model problems.

1. Consider the Cauchy problem

$$\frac{\mathit{du}}{\mathit{dx}}=-u,u\left(0\right)=1.$$E129

Let us replace the derivative with the forward difference,

$$\frac{\mathit{du}}{\mathit{dx}}\approx \frac{U_1\left(x\right)-U_1\left(0\right)}{x-0}=\frac{U_1\left(x\right)-1}{x},$$E130

In (130) $U_1\left(x\right)$ the approximate value of the unknown function at the moving point is if there is only one moving node.

Using (130) we write the difference Eq. (129)

$$\frac{U_1\left(x\right)-1}{x}=-U_1\left(x\right),$$E131

Take, now, two moving $x$ nodes and $x/2$. For these points, we write difference equations of the type (131)

$$\frac{U_2\left(x/2\right)-1}{x/2}=-U_2\left(x/2\right),\frac{U_2\left(x\right)-U_2\left(x/2\right)}{x-x/2}=-U_2\left(x\right),$$E132

Eliminating these equations $U_2\left(x/2\right)$, we get

$$U_2\left(x\right)=\frac{1}{{\left(1+x/2\right)}^2}.$$

For three moved nodes $x/3$, $2x/3$ and $x$ we get

$$U_3\left(x\right)=\frac{1}{{\left(1+x/3\right)}^3}.$$

If the number of nodes $n$, we get

$$U_n\left(x\right)=\frac{1}{{\left(1+x/n\right)}^n}.$$E133

If we strive for the number of nodes to infinity, we get

$$\underset{n\to \infty }{lim}{U}_n\left(x\right)=\underset{n\to \infty }{lim}\frac{1}{{\left(1+x/n\right)}^n}=e^{-x}.$$

Thus, we obtain the exact solution to problem (129).

2. Consider the problem

$$\frac{\mathit{dФ}}{\mathit{dx}}=\frac{1}{\mathit{Pe}}\frac{d^2Ф}{{\mathit{dx}}^2},Ф\left(0\right)=0,Ф\left(1\right)=1.$$E134

For this problem, the difference scheme with $2^k-1$ moving nodes has the form (29):

$$a_P^{\left(2^k-1\right)}{U}^{\left(2^k-1\right)}=a_E^{\left(2^k-1\right)}{U}_E^{\left(2^k-1\right)}+a_W^{\left(2^k-1\right)}{U}_W^{\left(2^k-1\right)}$$E135

where

$a_E^{\left(2^k-1\right)}=\frac{2^{2k+1}\left(1-{\gamma }_k\right)}{\left(1-x\right)\left(1-{\gamma }_k^{2^k}\right)},a_W^{\left(2^k-1\right)}=\frac{2^{2k+1}\mathit{Pe}\left(1-{\tau }_k\right)}{x\left(1-{\tau }_k^{2^k}\right)}+\frac{2^{2k+1}\left(1-{\tau }_k\right)}{x\left(1-{\tau }_k^{2^k}\right)},$$a_P^{\left(2^{л}-1\right)}=a_W^{\left(2^{л}-1\right)}+a_E^{\left(2^{л}-1\right)}.$ ${\tau }_k=2^k/\left(2^k+\sigma \right),{\gamma }_k=\left(2^k+\theta \right)/2^k$, $\theta =\mathit{Pe}\left(1-x\right).$

If we find from (135) $U^{\left(2^k-1\right)}$ and pass to the limit at $k\to \infty$, we have

$$\underset{k\to \infty }{lim}{U}^{\left(2^k-1\right)}\left(x\right)=\frac{e^{\mathit{Pe}x}-1}{e^{\mathit{Pe}}-1}.$$

The obtained limit coincides with the exact solution.

4. Conclusions

Summary: Derivation of approximate analytical solutions of differential equations by the moving nodes method

• The method of moving nodes allows one to obtain approximate analytical solutions for boundary value problems of mathematical physics.

• This is especially true in engineering applications, to obtain a rough analytical representation of the solution. The analytical method has its advantages over the numerical ones for its subsequent use and analysis of the structure of the solution.

• To refine the solution of differential equations, you can achieve this by adding the number of nodes to be moved.

• The examples given show the possibilities of applying and using the method of moving nodes for applied problems.

• Using the method of moving nodes based on the upwind scheme, compact schemes with high resolution for convective-diffusion problems are constructed.

Summary: Application of the moving node method

• Using the method of moving nodes based on the control volume method, compact schemes with high resolution for convective-diffusion problems are constructed.

• Using a combination of moving node methods and Richardson’s extrapolation, compact, high-resolution schemes for convective-diffusion problems are constructed.

• Choices of the influence of the profile on the faces of the control volume are studied.

• The possibilities of using movable nodes for the analysis of schemes are shown.

• Based on the method of moving nodes, the possibilities of finding errors in the approximation of differential equations are shown.

• For simple problems, the convergence of the moving nodes method is given.