On some numerical integration curves for PDE in neighborhood of "butterfly" catastrophe point Текст научной статьи по специальности «Математика»

URAL MATHEMATICAL JOURNAL, Vol. 2, No. 2, 2016

ON SOME NUMERICAL INTEGRATION CURVES FOR PDE IN NEIGHBORHOOD OF "BUTTERFLY" CATASTROPHE POINT1

Oleg Y. Khachay

Ural Federal University, Ekaterinburg, Russia, khachay@yandex.ru

Pavel A. Nosov

Ural Federal University, Ekaterinburg, Russia, pavel.nosov@urfu.ru

Abstract: We consider a three-dimensional nonlinear wave equation with the source term smoothly changing over time and space due to a small parameter. The behavior of solutions of this PDE near the typical "butterfly" catastrophe point is studied. In the framework of matched asymptotic expansions method we derive a nonlinear ODE of the second order depending on three parameters to search for the special solution describing the rapid restructuring of the solution of the PDE in a small neighborhood of the catastrophe point, matching with expansion in a more outer layer. Numerical integration curves of the equation for the leading term of the inner asymptotic expansion are obtained.

Key words: Matched asymptotic expansions, Numerical integration, Butterfly catastrophe, Nonlinear ODE and PDE.

Introduction

This paper is devoted to the study of specific behavior of a solution of the nonlinear wave PDE d2 U d2U d2U d2U

d2T + dT2+W2 + dz-2= /(eT'u) (0-1}

as a smoothed shock wave, the so called step-like contrast structure. Since this equation (0.1) contains 4 independent variables in every open domain of the arguments of function in the right-hand side of this equation there typically exists a point of catastrophe related to degeneration of f up to 5-th order with respect to the unknown function, that corresponds to "butterfly" type catastrophe [2].

The purpose of this paper is by applying the matched asymptotic expansions method [3] and catastrophe theory [2] to deduce the nonlinear ODE of the second order, which depends on three parameters, the ODE, which would be satisfied by the special solution related to a step-like contrast structure. We explore also the variants of such special solutions behavior depending on the parameter settings.

The similar equation with two independent variables and the corresponding typical point of "cusp" catastrophe was considered in [12]. The detailed study of special solutions with obtaining a uniform asymptotic expansion was carried out in [5], [6], this paper mainly follows the approach taken in these works.

1This work was supported by RFBR, research project No 16-31-00222.

1. Preliminary constructions

We consider first a more general than in equation (0.1) form of the differential operator in the left-hand side of PDE the resultant ODE in some cases will have the first order while in others it will be the second order. We suppose to study the constructed ODE in our following works in the framework of RFBR Research Project mentioned above. Consider a nonlinear PDE of such form

2

Y^ Aa(eR)daU(e, R) = f (eR, U(e, R)). (1.1)

|a| = 1

Following [13, p. 15] we denote by da the operator of differentiation with respect to independent variables and assume that the multi-index a = (a1,... , a4) corresponds to the independent variables, in particular, (R1,..., R4) = R in equality (1.1).

The relation (0.1) is a particular case of the equation (1.1), which corresponds to the following values of coefficients: Aa = 0 for 0 < an < 2, n = 1,..., 4 and 1 = —A2;o,o,o = Ao;2;o,o = Ao;o,2,o =

A-0,0,0,2.

In order to reduce the equation (1.1) to the standard form of singular equations with a small parameter multiplying the derivative, we make the change of variables s = eR, V(e, s) = U(e, eR), resulting equation (1.1) takes the form

2

e|a| Aa(s)daV(e, s) = f (s, V(e, s)). (1.2)

|a| = 1

Let the function f (s, V) be smooth and satisfy the inequality fV(s, V) = 0, (s, V) € x Qv. Let some conditions relevant to the equation (1.2) specify a solution on the domain asymptotically approximated by a series of the following form

ro

V(e, s) = Vo(s) + J] enVn(s), e ^ 0, (1.3)

n=1

on the domain Qs except for a neighborhood of its borders, and let the principal term Vo(s) of the series (1.3) satisfy the equation

f (s,Vo(s)) = 0. (1.4)

Examples of such statements of problems are presented, for instance, in the monograph [11] and in the paper [5].

At the boundary of could be located points of degeneracy of the function f (s, V) with respect to the unknown function, i.e., zeros of the function fV(s, V) on the boundary of the set which may contain the manifolds of points of "fold" catastrophe and "cusp" catastrophe, smooth lines of "swallowtail" catastrophe and isolated points of "butterfly" catastrophe [2].

2. Derivation of special equation

Let a point (s°, V°) on the boundary of the domain Qs x Qv be an isolated point of "butterfly" (A±5) type catastrophe [2, p. 5,11] of the function f (s, V). Then Taylor asymptotic expansion of f in the neighborhood of this point is

f (s, V) = £ (bm + Gm(s) + Gm(s)) (V — V°)m, m=o

where

4

Gm{s) — ^ bm,j (Sj — s°),

j=i

œ

Gm(s) y ^ bm,a(s - s°

|a|=2

0 — bo — bi — b2 — b3 — b4, (2.1)

b5 — K — 0. (2.2)

We introduce the following notation

B — {bn— 1,^n,j=1,...,4 .

The vanishing of the five coefficients (2.1) is achieved by the choice of the values of the five coordinates of a point (s°, V°); the existence of a solution of a system of five equations

dn

f(s,V)= 0, n = 0,..., 4

dV

with respect to five unknowns is a typical situation, while satisfying any additional independent

relations would not be the typical for functions f (s, V), and this, in particular, validates the assumption (2.2) and the inequality

det B = 0. (2.3)

Performing the change of variables

= -Gn-1(s), n = 1,..., 4; a = B (s - s°); V (e, s) = V (e, s° + B-1a) = W (e, a) (2.4)

and using the approach of catastrophe theory [2, p. 37,38,43], we find that there is a diffeomorphism of a neighborhood of (s°, V°) onto a neighborhood of the origin of variables (a, W) satisfying the asymptotic relation

V - V° — Ho(o-) + W(1 + Hi (a)) + J2 Wm(Cm + Hm(a)), (2.5)

m=2

œ

Hm(a) — ^ Cm,a&a |a| = 1

with a set of numerical coefficients Cm,Cm,a. Additionally, in the neighborhood of (ff°, V°) this diffeomorphism satisfies the identity

f (a, V) = -^i - - 2 - 3 + kW5, (2.6)

where the number k is the same as in (2.2), and the coefficients ^>n — ^>n(a) have the asymptotic representation

œ

^n(a) — an + ^ c,n,aaa. (2.7)

|a|=2

To study in details the behavior of the solution near the point of catastrophe s° corresponding to the value a — 0 we perform the coordinate stretching

an — Sn, Pn > 0, n — 1,..., 4; W (e, a)— eY W (e, S), y > 0 (2.8)

with some exponents (^1,..., p4) = p and y, the exact values of which will be defined below.

Since the leading term W of the series (2.5) has to satisfy the limiting equation (1.4) with expression (2.6) substituted into it instead of function f, then it turns out, that it depends on the all ^>n(ct), n = 1,..., 4, and hence, by virtue of the asymptotics (2.7), it is dependent on all an, n = 1,..., 4. To achieve this effect, when making the change of variables (2.8), we need to balance all the exponents of powers of e arising from the terms in equation (1.4). Setting

Pn = (n + 1)Yn, n = 1,..., 4, (2.9)

we obtain the asymptotic approximation

ro

f (s, V(e, s)) = e5Y (-S1 — S2W — S3W2 — S4W3 + kW5) + ^ emYPm(S, W), e ^ 0,

m=6

where W = W(e, S). Below, we will analyze the behavior of the principal term of the asymptotic expansion

ro

W(e, S) = w(S) + enwn(S), e ^ 0. (2.10)

n=1

To derive the principal term of the asymptotics of the left-hand side part of the equation (1.2) we transform the operators of differentiation with the use of changes of variables written in equalities (2.4), (2.8):

4

e^n dSn = d<rn = ^ bn-1 j dsj. j=1

To simplify the calculations we denote eYW(e, S) = W(e, S), then, as 4 d ~ ~ 4 d

S)) dSn = dW(e, S)) = dV(e, s) = s) ds3,

we have

whence

dSn ^ dsj

n=1 J=1 J

4 d ~ 4 4 d E £~^dS~nm£>s)) s ^ W ^ = E s)

n=1 n j=1 j=1 J

^ —/3 9

J n=1

2 2 4 / 4 „ \ aj

d

|a| = 1 |a| = 1 j=1 \n=1 n/

where

Aa(S) = Aa(s° + B-1CT(S)), CT(S) = (e^nSn)n=1,...,4 .

Consequently, under the condition of smoothness of functions Aa(s) and taking into account (2.8)-(2.10) at a fixed value of S and e ^ 0 the equation (1.2) can be written as

s) + °(£7)) + e2"2/3l+7 2W(£> s) + °(£7)) =

e5Y (—S1 — S2W — S3W2 — S4W3 + kw5 + O(eY)) , (2.11)

if we collect on the left-hand side in one bracket only the terms with derivatives of the first order and in another bracket - only with second-order derivatives; here

M = £ Q(a), N = £ Q(a),

|a| = 1 |a|=2

Q(a) = Aa(s°)ba, b = (60,1 ,...,60,4).

In particular, the constants M = 0 and N = — (60>1)2 + (60,2)2 + (60,3)2 + (60,4)2 correspond to the equation (0.1) and the constants M = — 60)1 and N = (60,2)2 + (60,3)2 + (60,4)2 are matched with the diffusion equation

dU d2U d2U d2U ^ ^ ^ ^ rn

~dT + ^ + ^ = /(eT'^

In accordance with the practice of matched asymptotic expansions method [3], when transiting to an internal scale, we need to choose 7 such, that the exponents of powers of e at the main terms of the asymptotics of the left-hand and right-hand sides of the equation coincide with each other after the transition to the new variables.

Let us consider the following two situations separately:

1) when the constant M = 0;

2) when in the original equation (1.1) all the coefficients of the first derivatives are identical to zero and the constant N = 0.

Note that implementation of these inequalities for the constant M and N, the value of which, except for special cases (for example, the lack of the first order in the original equation (1.1)), depends on the choice of f on the right-hand part of the original equation, is a typical situation similar to inequality (2.3), as noted above.

In the first case, by virtue of (2.9) and (2.11) we come to the conclusion, that 7 has to satisfy the relation

min{1 — 47, 2 — 97} = 57.

Solving it and taking into account the inequality 7 > 0, we obtain the value 7 = 1/9. Substituting it into the equality (2.11), we receive the estimate

d

M—w(S) + Si + S2w(S) + S3w2(S) + S4w3(S) - xw5(S) = Oie1/9). (2.12)

ds1

The limiting equation to 2.14 is an ODE of the first order with respect to w(S) as a function of one variable S1 and three parameters S2, S3, S4. By making in the ODE the linear change of variables

S1 = (M5k-1)1/9 x, S2 = (M4K)1/9 y, S3 = (MK)1/3 z, S4 = (M2K5)1/9 t, u(x) = (M k_2)1/9 w(S1 (x), S2(y), S3(z), S4(t)),

we obtain the first order nonlinear ODE

ux = u5 — tu3 — zu2 — yu — x, (2.13)

which depends on three parameters y, z, t. We plan the study of the behavior of solutions of ODE (2.13) to be hold in subsequent papers in the framework of the above mentioned RFBR research project.

In the second case, in the relation (2.11) the first term of the equation is completely missing and the equation for 7 becomes simpler 2 — 97 = 57, whence 7 = 1/7 and therefore equality (2.11) takes the form:

d2

Nt^t^w(S) + Si + S2w(S) + S3w2(S) + S4w3(S) - xw\S) = O^1/7). (2.14) d (S1)2

We turn to the equation obtained from (2.14) as a result of passing to the limit as e ^ 0 and at the same time producing a linear change of variables

S1 = — sgn(x) \N5x-1|1/14 x, S2 = — sgn(x) \N2x|1/7 y, S3 = — sgn(x) \N3k5\1/14 z, S4 = — sgn(x) \Nk4\1/7 t, u(x) = \N k-3\1/14 w(S1(x),S2(y),S3(z),S4 (t)).

Thus, we obtain the desired nonlinear second-order ODE

sgn(N) uxx = u5 — tu3 — zu2 — yu — x, (2.15)

depending on three parameters y, z, t. This ODE is held for a special solution describing the rapid reconstruction of the original PDE solution in a small neighborhood of the catastrophe point s°.

3. Matching with more outer layer condition

Consider the multi-valued relation

u = h(x; y, z, t), (3.1)

each value of which is a root of the equation

0 = u5 — tu3 — zu2 — yu — x.

In order to apply the matched asymptotic expansions method, it is necessary to construct a function u(x) = u(x; y,z,t), which is a solution of the equation (2.15), i.e., it is the principal term of the asymptotic expansion of the solution of the original problem for PDE in the inner layer (the layer, projection of which onto the axis Ox is the set |x| < e-2/7+<51) in such way that it would be matched with the solution of the original problem in the more outer layer (projection of which onto the axis Ox is the set |x| > e-2/7+^2), where 0 < ¿1 < ¿2 < 2/7 are some numbers. Therefore it is necessary to match u(x; y,z,t) and Vo(s), the leading term of the series (1.3), which satisfy the equation (1.4). As a consequence of that we will look for those solutions u(x) = u(x; y,z,t) of equation (2.15), which satisfies the limiting relation

lim |u(x; y,z,t)+ H (x; y,z,t)| = 0, (3.2)

where u = H(x) = H(x; y, z, t) stands for the maximal extension over axis Ox in the right-hand direction of a smooth branch of the root (3.1) defined in a neighborhood of x = —to. Thus, we will configure the function u(x; y, z, t) to match with the Vo(s) on the left-hand side along axis Ox and then we will study the behavior of such a curve with an increase in the variable x.

4. Numerical search for special integral curves

Integral curves of solutions given in this paper shown in bold and green on the figures below were calculated using explicit Runge-Kutta-Felberg (4,5) method [1] with variable step and accuracy control. In all figures a red thick line represents the graph of the root (3.1).

Note that the behavior of solutions of equations (2.15) significantly depends on the sgn(x) in the left-hand side. For minus sign it is typical to appear rapid fluctuations of solutions after passing the branching point of a multi-valued relation (3.1). Figure 1 illustrates the occurrence of such fluctuations of the integral curve of the equation

—uxx = u5 + 9u3 - 8u - x, (4.1)

computed under condition (4.6). In this work we do not perform a detailed study of such solutions

6 5 4 3 2 1

u 0 -1 -2 -3 -4 -5 -6

Figure 1. Integral curve for equation (4.1) under condition (4.6) and graph of root line (3.1)

with fluctuations; we focus on so-called step-like contrast structures. Therefore, we will consider only the equation with plus sign:

uxx = u5 — tu3 — zu2 — yu — x. (4.2)

The purpose of construction and further study of the behavior of the integral curves was to test the hypothesis about the possibility of proving the existence of such special solutions within the framework of the approach of [5], [6]. Therefore, at this stage, we have considered the case when the parameter z in the right-hand side of equation (4.2) is zero, since this option, on the one hand, includes the very point of "butterfly" catastrophe (at the origin of all variables and parameters). On the other hand, it is very close to "cusp" catastrophe type of discussed in these papers, since in this case, the root line (3.1) is odd (symmetric with respect to the origin), and therefore the solution can be constructed as the extension based on oddness onto the negative half-axis (x < 0) of a solution u(x) of the Cauchy problem for the equation

uxx = u5 — tu3 — yu — x (4.3)

with the initial conditions

u(0) = 0, ux(0) = a.

(4.4)

(4.5)

In this case, the limiting relation (3.2) takes the form

lim |u(x; y, 0,t)+ H(x; y, 0, t) | = 0. (4.6)

x^+ro

Thus, below, we consider integral curves only for the Cauchy problem (4.3)-(4.5). As a problem of obtaining the principal term of the internal expansion of the original problem for PDE it presents the extension of study conducted by one of the authors of this paper of bisingular initial problems with one small parameter for the systems of one, two or more of ODEs, which also has the property of degeneration of high-order of right-hand part of the equation in respect to the unknown function [4, 7-10].

According to the theorem proved in [6], one can find the initial value of ao of the derivative ux(x) corresponding to the solution that satisfies the asymptotics (4.6) by calculating the exact lower boundary for a certain set of M. The numerical set M depending on the parameters of the differential equation was determined in [6] as follows: the number of a belong to the set M, if and only if the solution of the Cauchy ua(x) corresponding to initial condition (4.5) with this value of a has the property, that there exists a point xa such that

• for 0 function u«(x) is less than the function H(x);

• at the point xa functions ua and H(x) coincide.

It is not difficult to see that the most top branch of the root (3.1) defined by the formula u = H(—x), to which the desired solution has to approach (see, eg., Fig. 4.a), consists by virtue of the equation (4.3) (except, perhaps, for its leftmost point) of points of repulsion, which damps on approaching the line u = H(—x) and rapidly increases with increasing the distance from this line in any vertical direction. Therefore, most of the curves produced with the initial conditions (4.4), (4.5), either pass through a line of u = H(—x), and quickly grow up, tending to +to, or, if the initial velocity (4.5) is not sufficient, can not come close to this line and are broken down to —to. Search for the fine line of the balance between these two states of computed integral curves is consistent with the idea of the works [5], [6] of the initial value of ao = inf M in (4.5).

We have implemented a simple binary search algorithm of this balance: in the transition from the range [am,6m] of possible values ux(0) at the current step to the next step range [am+1,6m+1] we always choose one of the intervals [am,cm], [cm,6m] (where cm is the middle of the segment [am,bm]), for which the integral curve started with an initial rate of am+1 breaks down to —to, and one with initial rate of bm+1 grows up to +to. Clearly, the value pretend to be an element a of the set M described above and the limit of the decreasing sequence could play the role of ao = inf M.

5. Illustrations of section at z = 0 of separatrix of "butterfly" catastrophe (A±5)

This section of the paper is complementary, although it still has an indirect relationship to the purpose of this paper. Here, for the first time, as far as we know, the three-dimensional illustration of the cross section for z = 0 separatrix [2, p. 51] "butterfly" catastrophe (A±5) are given. Separatrix corresponding to such a point of catastrophe is the set of points (x°,y°, z°, t°) included into the parameter space of the parametric family of functions

■0(u; x, y, z, t) = u5 — tu3 — zu2 — yu — x, (5.1)

with the following property: when such a point (x°, y°, z°, t°) is substituted into (5.1) the resultant function ^°(u) has at least one non-Morse point, i.e., the following relation holds

3u°: (ip°(u°) = 0, = 0)-

The separatrix [2, p. 57-58] divides the parameter space into open domains corresponding to which subfamilies of the general family of functions (5.1) are structurally stable.

To construct a three-dimensional surface data we implemented the approach applied in [2, p. 6263] to the function of "swallowtail" catastrophe (A4). Knowledge about the configuration of the

Figure 2. Section at z = 0 of separatrix of family of functions (5.1)

section of the separatrix allows us to study not all possible options for the behavior of curves at z = 0 through taking only one representative from each of the open domains of the parameters (x, y, t) corresponding to Morse functions from the family (5.1) and additional representatives from the border separatrix itself, detailed analysis of which can be done by ranking the degree of degeneracy of the function (5.1) with respect to u. Some of the chosen within the framework of this concept values of the parameters (x,y,t) are given in Table 1. The figures with the corresponding integral curves of the Cauchy problem (4.3)-(4.5) illustrating the results of the described above algorithm of the binary search for the solution of the asymptotic problem (4.3), (4.4), (4.6) are placed in the next section of the paper.

Table 1. Parameter values and corresponding numbers of figures

No y t Corresponding figures

1 -8 9 Fig. 3, 4

2 -10 5 Fig. 5

3 8 -9 Fig. 6

4 8 8 Fig. 7

5 -8 -8 Fig. 8

6 8 0 Fig. 9

7 0 0 Fig. 10

6. Illustrations of numerical calculations

This section is devoted to illustrations of the integral curves for the Cauchy problem (4.3)-(4.5) obtained using the binary search for the initial speed as described above. We remind that in the figures bold green line represents the integral curves and the red bold line corresponds to the graph of the root function (3.1). In figures with contour distributions the increase in the intensity of green shade in the color of lines corresponds to the increase of positive values of the function, and increase in the intensity of blue shade in the color of lines corresponds to the increase in the absolute value of the negative values. The contour distributions are given in the paper to show the rate of the function (5.1) changing with increasing distance from the line of the root (3.1) in different parts of the line.

(4.3)-(4.5) for y = —8, t = 9

The root line (3.1) for y = —8, t = 9 and z = 0 shown in Fig. 3 and 4 has an explicit double system of bends, in virtue of which the integral curves of the Cauchy problem (4.3)-(4.5) obtained for different values of the initial velocity demonstrate a variety of processes including fluctuations. Partial figures of the Fig. 3 marked with letters from a) to f) correspond to selected stages of the

binary search, which embodied such variants of the curve behavior as the gradual disappearance of the interval of oscillation and the slow emergence and expand of the interval, on which the curve passes very close to the line of the root (3.1).

3 2 1 u 0 -2 -3 3 2.5 2 1.5 1 0.5 u 0 -0.5 -1.5 -2 -2.5 -3 )

-30 -20 -10 0 10 20 30 a) -35 -30 -25 -20 -15 -10 -5 b) 5 10 15 20 25 30 35

Figure 4. Results obtained for y = -8, t = 9 z = 0: a) contour distributions of function (5.1); b) integral curves for the Cauchy problem (4.3)-(4.5) for ux(0) = 7.38348579007552236

2 1- u 0 -2 2.521.510.5-u0--0.5-1-1.5-2-2.5-

-40 -20 0 20 40 a) -40 -30 -20 -10 0 10 20 30 40 b)

Figure 5. Results obtained for y = -10, t = 5 z = 0: a) contour distributions of function (5.1); b) integral curves for the Cauchy problem (4.3)-(4.5) for ux(0) = 0.10031199880359876075

The situation that occurs when y = -10, t = 5 and z = 0 is shown in Fig. 5. It is characterized by the presence of two deflections of green curve, which is the approximation of the desired solution of the asymptotic problem (4.3), (4.4), (4.6), from the root line (3.1) in various ways: in the beginning integral curve follows the line of the root, then exceeds it, then crosses and passes below the root line and later again starts to follow it.

Situation in Fig. 6 that occurs when y = 8, t = —9 and z = 0, is very similar to the pattern for case of "cusp" catastrophe studied in the works [5], [6], the resemblance in the shape of the root lines (3.1) and of integral curves is explicit.

If y = 8, t = 8 and z = 0, as shown in Fig. 7, the root line (3.1) has a larger number inflection points, than one in the case of "cusp" catastrophe or than in Fig. 6, but the behavior of the integral curve remains broadly similar: a sharp increase is replaced by following the root line.

The graph the function (3.1) in Fig. 8 for y = —8, t = —8 and z = 0 increases gradually and remains single-valued for all points of positive half-axis (x > 0). By virtue of this fact, immediately from the very point x = u = 0, the root line (3.1) becomes an attractor for the approximation curve for the solution of asymptotic problem (4.3), (4.4), (4.6) and the integral curve following it until the final rapid movement to the vertical infinity.

1.5 0.5 0 -0.5 -1.5 1.61.41.210.80.60.40.2-u0--0.2-0.4-0.6-0.8-1-1.2-1.4-1.6-

-30 -20 -10 0 10 20 30 a) -35 -30 -25 -20 -15 -10 -5 0 5 b) 10 15 20 25 30 35

Figure 6. Results obtained for y = 8, t = —9 z = 0: a) contour distributions of function (5.1); b) integral curves for the Cauchy problem (4.3)-(4.5) for ux(0) = 1.982684405999750261

3 2 1 u 0 -2 -3 s 32.521.510.5-u0--0.5-1-1.5-2-2.5-3-

-60 -40 -20 0 20 40 60 a) -60 -50 ^10 -30 -20 -10 0 10 20 30 40 50 60 b)

Figure 7. Results obtained for y = 8, t = 8 z = 0: a) contour distributions of function (5.1); b) integral curves for the Cauchy problem (4.3)-(4.5) for ux(0) = 12.406010731733428955

1.5 0.5 u 0 -0.5 -1.5 1.61.41.210.80.60.40.2-u0--0.2 --0.4 --0.6 --0.8-1.2 --1.4 --1.6 -

-60 -40 -20 0 20 40 60 a) 60 -50 -10 -30 -20 -10 0 10 20 30 40 50 60 b)

Figure 8. Results obtained for y = —8, t = —8 z = 0: a) contour distributions of function (5.1); b) integral curves for the Cauchy problem (4.3)-(4.5) for ux(0) = 0.1236970813110499762

The situation in Fig. 9, corresponding to the values y = 8, t = 0 and z = 0 of the parameters and therefore to the absence of not only a quadratic but also a cubic term in the formula (5.1), the shape of the root line (3.1) has a more extended than in Fig. 6, the interval of domination of the linear term near the origin, but the behavior of the integral curve remains broadly the same as in the Fig. 6.

The values of the parameters for which the curves in Fig. 10 are constructed, correspond to the

2 1- u 0 -2 2 1.5 1 0.5 u 0 -0.5 -1.5 -2

-30 -20 -10 0 10 20 30 a) 5 -30 -25 -20 -15 -10 -5 0 5 10 15 20 25 30 35 b)

Figure 9. Results obtained for y = 8, t = 0 z = 0: a) contour distributions of function (5.1); b) integral curves for the Cauchy problem (4.3)-(4.5) for ux(0) = 4.0017096188039309541127989

2 1- u 0 -2 21.510.5- u0--0.5-1-1.5-2-

-40 -20 0 20 40 a) -40 -30 -20 -10 -5 0 5 10 15 20 25 30 35 40 b)

Figure 10. Results obtained for y = 0, t = 0 z = 0: a) contour distributions of function (5.1); b) integral curves for the Cauchy problem (4.3)-(4.5) for ux(0) = 0.9725146994563902

origin of the parameter space of the function (5.1). In Fig. 10 the point x = u = 0 is the very point of "butterfly" catastrophe. One can note that in a neighborhood of this point the integral curve slightly deviates from almost vertically growing root line.

Thus, the results of calculations for the values of the parameters in Table 1 show that, gradually increasing the accuracy of the curve fitting method [1] and using a binary search algorithm, it is possible to obtain curves, which for more and more long time remain close to the line of balance noted above, that is consistent with hypothesis about the construction of the desired curve satisfying the asymptotic condition (4.6) as the solution of the Cauchy problem (4.3)-(4.5) with the value a0 = inf M of the initial rate.

7. Conclusion

We considered a PDE with an arbitrary linear combination of differentiation operators with smooth coefficients of the first and second orders on the left-hand side of the equation and a nonlinear function in the right-hand side, which depends on the desired function and contains a small parameter. We consider and analyze also the three-dimensional nonlinear wave equation with the source term smoothly changing over time and space.

The primary study of the behavior of solutions of this PDE near the typical point of "butterfly" catastrophe was held. We deduce two nonlinear ODE of the first and the second orders, respectively, depending on three parameters. The order of equation depends on the configuration

of the coefficients of the linear combination on the left-hand side of the original PDE and of the properties of the nonlinear function on its right-hand side.

We use the resulting second-order ODE to search for a special solution describing the rapid reconstruction of the solution of the wave PDE in a small neighborhood of the catastrophe point matched with expansion in a more outer layer. We have done a primary study: we produced integral curves that allow one to analyze the behavior of such a special solution. The results revealed no contradictions with the possibility to prove the existence of special solution within the framework of the approach given in [5,6].

Acknowledgements

This work was supported by RFBR, research project No. 16-31-00222.

REFERENCES

1. Fehlberg E. Low-order classical Runge-Kutta formulas with stepsize control// NASA Technical Report R-315, 1969.

2. Gilmore R. Catastrophe theory for scientists and engineers// New York: Dover Publications, 1993. 666 p.

3. Il'in A.M. Matching of asymptotic expansions of solutions of boundary value problems// Transl. Math. Monogr., Vol.102, Providence: Amer. Math. Soc., 1992. 281 p.

4. Il'in A.M., Leonychev Y.A., Khachay, O.Y. The asymptotic behaviour of the solution to a system of differential equations with a small parameter and singular initial point// Sbornik Mathematics, 2010. Vol. 201, no. 1, P. 79-101.

5. Il'in A.M., Suleimanov B.I. On two special functions related to fold singularities// Doklady Mathematics, 2002. Vol. 66, no. 3, P. 327-329.

6. Il'in A.M., Suleimanov B.I. Birth of step-like contrast structures connected with a cusp catastrophe// Sbornik Mathematics, 2004. Vol. 195, no 12. P. 1727-1746.

7. Khachay O.Y. Asymptotic expansion of the solution of the initial value problem for a singularly perturbed ordinary differential equation// Differential Equations, 2008. Vol. 44, no. 2, P. 282-285.

8. Khachay O.Y. Asymptotics of the solution of a system of nonlinear differential equations with a small parameter and with a higher-order extinction effect// Differential Equations, 2011. Vol. 47, no. 4, P. 604-607.

9. Khachay O.Y. On the matching of powerlogarithmic asymptotic expansions of a solution to a singular Cauchy problem for a system of ordinary differential equations// Trudy Instituta Matematiki i Mekhaniki, 2013. Vol. 19, no. 1, P. 300-315. [in Russian]

10. Khachai O.Y. On the application of the method of matching asymptotic expansions to a singular system of ordinary differential equations with a small parameter// Differential Equations, 2014. Vol. 50, no. 5, P. 608-622.

11. Maslov V.P., Danilov V.G., Volosov K.A. Mathematical modeling of heat and mass transfer processes. Evolution of dissipative structures// Nauka. Moscow, 1987. 352 p. [in Russian]

12. Suleimanov B.I. Cusp catastrophe in slowly changing equilibrium positions// JETPh, 2002. Vol. 122, no. 5 (11). P. 1093-1106. [in Russian]

13. Vladimirov V.S., Jarinov V.V. Equations of mathematical physics// Fizmatlit, Moscow, 2004. 400 p., ISBN: 5-9221-0310-5. [in Russian]