(Hokusai, The Great Wave off Kanagawa, 1823)

Let's consider a two-dimensional medium consisting of particles moving by inertia, without intersection and in the absence of external forces. Denote $u(t,x)$ is the speed of a particle located at a point in time $t$ and space point $x$. If $x=q(t)$ is the trajectory of motion of some fixed particle, then its speed $q'(t)=u(t, q(t))$, and acceleration $q''(t)$ is equal to zero.

The resulting equation describing the evolution of the velocity field of non - interacting particles is called Hopf equation

The analytical solution of Hopf equation

The analytical solution of Hopf equation can be determined by method of characteristics:

where, $g$ is the continuous function.

It is necessary to invert the map $x = x_{0} + u_{0}(x_{0})t$ concerning $x_{0}$, and then substitute it into the expression for $u$ to find the explicit dependence $u = u(x,t)$ (the parametric solution), which is most often found in literature:

The solution of the Hopf equation can be easily analyzed. According to the solution, each point of the velocity profile moves with its own (constant) velocity. In this case, the velocity profile will be deformed since faster particles will overtake the slower ones. As a result, the shape is steeper, and at time $t = t^* $, the gradient of $x$ becomes endless. This phenomenon is called overturning or gradient catastrophe [2]. To find the overturning time $t^* $, it is necessary to calculate derivative $\frac{\partial{u}}{\partial{x}}$:

where is the Jacobian of $x = x_{0} + u_{0}(x_{0})t$. It can be seen from this expression that the derivative becomes infinite for the first time when the Jacobian, $J=0$ and corresponding time can be expressed as the following formula:

This minimum corresponds to the point x* (the breaking point). From this, it is clear that the overturning occurs for the initial profiles with a negative derivative $u_x < 0$. For $t > t^* $, the equation $J = 0$ has two roots, respectively. The derivative goes to infinity at two points. The area of ambiguity is between these points. From the mapping point of view, this means that the mapping $x = x(x_{0}, t)$ is single-valued — the equation $x = x(x_{0}, t)$ has only one (real) root $x_{0} = x_{0}(x, t)$ for all $x$. For $t > t^* $, the equation $x = x(x_{0} , t)$ has three roots $x_{0}$ for any $x$ from the ambiguity domain. A change in the number of roots of a mapping is called a bifurcation. In this case, one root converts to three roots at the bifurcation point . The bifurcation point is also called the cusp point. In physics, overturning is the cause of the formation of shock waves. When approaching the overturning point, the equations of ideal hydrodynamics lose their applicability in real gas.

Model description

We study the Hopf equation with following limits: 0 < t < 1 , -5 < x < 10, and the initial condition:

Using the the analytical solution of this problem is represented by the following formula:

The overturning effect at which $u_{x}$ turns to infinity and the function ceases to be single-valued arises at the moment, $t^{*}=1.649$:

We use Newton's method for two different numerical schemes (Lax-Friedrichs and MacCormack schemes) to solve the Hopf equation with the following limits and the initial condition.

Newton's method (Newton-Raphson method) is an iterative numerical method for finding the root of a given function [3, 4]. The search for a solution is carried out by constructing successive approximations based on simple iteration principles.

Description

To numerically solve the Eq. 10, by simple iteration it must be reduced to the equivalent expression Eq. 11.

where, $\phi$ is the contraction display. The following condition:

is used for the best convergence of the method at the next approximation point $y^*$. The solution to this equation is sought in the form:

that implies:

Assuming that the point of approximation is "close enough" to the root $\bar {y}$ and that the given function is continuous:

The final formula for $\alpha(y)$ is:

The function $\phi(y)$ is defined:

Algorithm

The algorithm for Newtons method can be described as follows:

1. Set the initial approximation $y_0$

2. Until the stop condition is met, which can be taken as $|y_{n+1}-y_{n}|&lt;\epsilon$, the error is within the required limits, calculate a new approximation:

We utilize Lax - Friedrichs and MacCormack numerical schemes using Newton's method to determine the numerical solution.

Lax - Friedrichs

The Lax-Friedrichs scheme is defined in numerical analysis as a method for the numerical solution of hyperbolic partial differential equations based on the finite difference method [5].

Extensions to Nonlinear Problems

The formula defines the hyperbolic system of one-dimensional conservation laws:

A generalization of the Lax - Friedrichs scheme to nonlinear systems can be expressed in the following form:

Lax - Friedrichs numerical scheme can be used to construct higher-order schemes for solving systems of hyperbolic partial differential equations, in the same way that Euler's method is used to construct more accurate Runge-Kutta method for solving ordinary differential equations [6].

This scheme can be written in a conservative form:

The Lax - Friedrichs scheme is explicit and has approximation error is $O(\Delta t, \Delta x)$. For this reason, it is not used as often in practice. This scheme is stable under the following condition:

MacCormack scheme

Difference schemes "predictor-corrector"

Predictor-corrector schemes are a family of methods related to algorithms designed to integrate ordinary differential equations [7, 8]. All such techniques involve two steps [9]:

1. At the first step (predictor), some function is determined from values calculated in the previous step to get the approximated value of the desired function in the following point.

2. At the second step (corrector), the received approximation using the predicted value function and another operator to interpolate the value desired position at the same point.

The following steps can represented as the following expressions:

where, $u^{n+1/2}$ is the value obtained at the "predictor" step, and $u^{n+1}$ is the refined desired value. It can be seen that the first step is implemented using explicit methods, and the second step is based on the application of formulas of implicit methods. On the right side, instead of the unknown value $u^{n+1}$, the result is substituted predictions $u^{n+1/2}$.

Methods using the predictor-corrector scheme:

1. Milne method for ODE.

2. Heun method (predictor - Euler's method, corrector - Trapezium method).

3. Adams-Bashforth method for solutions non-rigid boundary value problems (the Adams-Bashforth-Moulton corrector is used).

4. MacCormack method.

MacCormack method

The MacCormack method is a modified two-step scheme Lax - Wendroff, but it is much easier to use [10]. Consider the following first-order hyperbolic equation:

Predictor: At this step, the predicted value of $w$ at the moment time $k+1$: $\bar{w}^{k+1}_{i}$ is evaluated as follows:

Corrector: At this step, the predicted value $\bar{w}^{k+1}_{i}$ corrected according to the equation:

Note that the intermediate value $\bar{w}^{k+1}_{i}$ doesn't make sense physically. The scheme has the second order of accuracy with an approximation error $O(\Delta t^2, \Delta x^2)$, while it is stable at $|a\frac{\Delta t}{\Delta x}|\leq1$. The differential approximation of the scheme has the following form:

The MacCormack scheme is often used due to a number of its advantages [11]. In particular, it operates only with quantities in the primary grid nodes and can be easily generalized to multidimensional problems. Also, this is the schematic second-order accuracy.

Lax - Friedrichs scheme

We convert the Hopf equation into the divergent form:

and using Eq. 21 we can get the Lax - Friedrichs numerical scheme:

where $dt$ and $dx$ are discretization steps in the time and space domain correspondingly. The upper index defines the time location, and the lower represents the spatial location of the function. The scheme requires a boundary condition on the left side (x=0). We use $u^n_{0}=0$ because of the exponential damping of the initial function.

We studied the solution of the Hopf equation for different values of time: $t = 0.5, 1, 1.649, 1.8$. The Fig. 1 compares analytical and numerical solutions for a given time. In addition, the maximum discrepancy $r$ between the analytical and numerical solutions was calculated for the Lax - Friedrichs scheme.

The Fig. 2 shows the 3D numerical solution of the Hopf equation for $t = 1$ and $t = 1.65$ As can be seen, for the time $t=1.65$ overturning effect occurs, the velocity profile becomes steep.

MacCormack scheme

Similarly, as in the previous numerical scheme, we used the divergent form of the differential equation. We calculate the predictor using the following formula:

Later, to determine the corrector, the predicted value of the function we approximated using the following formula:

The Fig. 3 shows the 3D numerical solution of Hopf equation using MacCormack scheme. In the same way as in the Lax - Friedrichs scheme, numerical solutions were constructed for different values of time $t$ (Fig. 4) for MacCormack scheme. We use the number of time discretization points $N_{t} = 100$ and the spatial discretization, $N_{x} =1500 $, for both numerical schemes. Fig. 1 and Fig. 4 show that the analytical and numerical solutions coincide very well. However, the maximum discrepancy $r$ when using the MacCormack scheme is less than for the Lax - Friedrichs scheme (Fig. 5). It can be explained by the fact that the MacCormack scheme has the second order of accuracy $O(\Delta t^2, \Delta x^2)$ in contrast to the Lax - Friedrichs scheme $O(\Delta t, \Delta x)$. Also, it can be seen that the maximum discrepancy $r$ is increasing with increasing time $t$. Up to the overturning time, $t^* $, the difference rises very slightly. After this time, the maximum value of discrepancy increases significantly. It is due to the overturning effect. The solution becomes ambiguous after reaching the overturning time $t^* $.

References

[1] Kuznetsov, E., and Shapiro, D., 2011. “Methods of mathematical physics: a course of lectures”. Novosibirsk State University.

[2] Kuznetsov, E., and Mikhailov, E., 2022. “Slipping flows and their breaking”. Annals of Physics,p. 169088.

[3] Mor ́e, J. J., and Sorensen, D. C., 1982. Newton’s method. Tech. rep., Argonne National Lab., IL (USA).

[4] Polyak, B. T., 2007. “Newton’s method and its use in optimization”. European Journal of Operational Research, 181(3), pp. 1086–1096.

[5] DuChateau, P., and Zachmann, D. W., 2002. Applied partial differential equations. Courier Corporation.

[6] CHU, C., 1978. Numerical methods in fluid dynamics, in “advances in applied mechanics”(cs. yih, ed.), vol. 18.

[7] Zhang, P.-G., and Wang, J.-P., 2012. “A predictor - corrector compact finite difference scheme for burgers’ equation”. Applied Mathematics and Computation, 219(3), pp. 892–898.

[8] Butcher, J. C., 2016. Numerical methods for ordinary differential equations. John Wiley & Sons.

[9] Press, W., Teukolsky, S., Vetterling, W., and Flannery, B., 2007. “Section 17.6. multistep, multivalue, and predictor-corrector methods”. Numerical Recipes: The Art of Scientific Computing.

[10] Anderson, J. D., and Wendt, J., 1995. Computational fluid dynamics, Vol. 206. Springer.

[11] Hixon, R., and Hixon, R., 1997. “On increasing the accuracy of maccormack schemes for aeroacoustic applications”. In 3rd AIAA/CEAS Aeroacoustics Conference, p. 1586.