The Analytical and a Higher-Accuracy Numerical Solution of a Free Boundary Problem in a Class of Discontinuous Functions

Abstract

A new method is suggested for obtaining the exact and numerical solutions of the initial-boundary value problem for a nonlinear parabolic type equation in the domain with the free boundary. With this aim, a special auxiliary problem having some advantages over the main problem and being equivalent to the main problem in a definite sense is introduced. The auxiliary problem allows us to obtain the weak solution in a class of discontinuous functions. Moreover, on the basis of the auxiliary problem a higher-resolution numerical method is developed so that the solution accurately describes all physical properties of the problem. In order to extract the significance of the numerical solutions obtained by using the suggested auxiliary problem, some computer experiments are carried out.

1. Introduction

It is known that many practical problems such as distribution of heat waves, melting glaciers, and filtration of a gas in a porous medium, and so forth, are described by nonlinear equations of the parabolic type.

In [1], at first the effect of localization of the solution of the equation describing the motion of perfect gas in a porous medium is observed and the solution in the traveling wave form is structured. Then, the mentioned properties of the solution for the nonlinear parabolic type equation are studied in [2, 3], and so forth.

These problems are also called free boundary problems. Therefore, it is necessary to obtain the moving unknown boundary together with the solution of a differential problem. Its nature raises several difficulties for finding analytical as well as numerical solutions of this problem.

The questions of the existence and uniqueness of the solutions of the free boundary problems are studied in [4, 5]. In [5], Oleĭnik introduced the notion of a generalized solution of the Stefan problem whose uniqueness and existence were guaranteed in the class of measurable bounded functions. In [4], Kamenomostskaya considered the classical quasilinear heat conduction equation and constructed the generalized solution by the use of an explicit difference scheme.

In the literature there are some numerical algorithms (homogeneous schemes) which are approximated by finite differences of the differential problem without taking into account the properties occurring in the exact solution [6–8].

2. Traveling Wave Solution of the Main Problem

We consider the equation 𝜕𝑢=𝜕𝜕𝑡²𝜑(𝑢)𝜕𝑥²,in𝑅²₊ (2.1) with following initial 𝑢(𝑥,0)=𝑢₀(𝑥)=0,in𝐼=[0,∞) (2.2) and boundary 𝑢(0,𝑡)=𝑢₁(𝑡)=𝜇₀𝑡^𝑛,𝑡>0 (2.3) conditions, where 𝑅²₊=𝐼×[0,𝑇). Here, 𝜇₀ and 𝑛 are real known constants. In order to study the properties of the exact solution of the problem (2.1)–(2.3) for the sake of simplicity the case 𝜑(𝑢)=𝑢^𝜎 is considered. It is clear that the function 𝜑(𝑢) satisfies the following conditions:(i)𝜑(𝑢)∈𝐶²(𝑅²₊),(ii)𝜑^′(𝑢)≥0, for 𝑢≥0 and 𝜎≥2,(iii)for 𝜎≥2, 𝜑^′′(𝑢) have alternative signs on the domain when 𝑢(𝑥,𝑡)≠0.

It is easily shown that the problem (2.1)–(2.3) has the solution in the traveling wave form as 𝑢()=𝐷𝑥,𝑡𝜎−1𝜎^1/(𝜎−1)(𝐷𝑡−𝑥)^1/(𝜎−1),0<𝑥<𝐷𝑡,0,𝑥≥𝐷𝑡. (2.4) Via simple calculation we get the following:(1)the function 𝑢(𝑥,𝑡) and 𝑤(𝑥,𝑡)=−𝜕𝑢^𝜎(𝑥,𝑡)/𝜕𝑥=𝐷𝑢(𝑥,𝑡) are continuous in 𝐷_ℓ(𝑡)={(𝑥,𝑡)∣0≤𝑥≤𝐷𝑡,0≤𝑡≤𝑇}, but 𝑢_𝑡 and 𝑢_𝑥 do not exist when 𝜎>2;(2)when 𝜎=2, 𝑢_𝑡 and 𝑢_𝑥 are finite;(3)for 1<𝜎<3/2, all derivatives 𝑢_𝑡, 𝑢_𝑥, and 𝑢_𝑥𝑥 exist;(4)for 3/2<𝜎<2 the 𝑢_𝑡 and 𝑢_𝑥 exist, but 𝑢_𝑥𝑥 does not exist.

Therefore when 𝜎>2, we must only consider the weak solution for the problem (2.1)–(2.3). As it is seen from the formula (2.4), when 𝑛=1/(𝜎−1), the solution is equal to zero at 𝑥=ℓ(𝑡)=𝐷𝑡. Here, 𝐷=±(𝜎/(𝜎−1))𝜇₀^𝜎−1 is the speed of the front of the traveling wave and 𝑢(ℓ(𝑡),𝑡)=0. (2.5) Taking into account (2.5), we get 𝑑𝑥𝑑𝑡=−𝜕𝑢/𝜕𝑡=()𝜕𝑢/𝜕𝑥𝐷𝜕𝑢𝑥,𝑡/𝜕𝑥𝜕𝑢/𝜕𝑥=𝐷. (2.6) It is clear that when 𝑢=0, (2.1) degenerates to the first-order equation. The weak solution is defined as follows.

Definition 2.1. A nonnegative function 𝑢(𝑥,𝑡) which is satisfying the initial condition (2.2) and boundary conditions (2.3) and (2.5) is called a weak solution of the problem (2.1)–(2.3), and (2.5), if the following integral relation _{𝐷ℓ(𝑡)}𝑢()()𝑥,𝑡𝜕𝑓𝑥,𝑡−𝜕𝑡𝜕𝑢^𝜎(𝑥,𝑡)()𝜕𝑥𝜕𝑓𝑥,𝑡𝜕𝑥𝑑𝑥𝑑𝑡=0 (2.7) holds for every test functions 𝑓(𝑥,𝑡) from 𝐶^𝑜_1,1(𝐷_ℓ(𝑡)) and 𝑓(𝑥,𝑇)=0.

Because the function 𝑢(𝑥,𝑡) and 𝑤(𝑥,𝑡)=−𝜕𝑢^𝜎(𝑥,𝑡)/𝜕𝑥 are continuous in the domain 𝐷_ℓ(𝑡), the integral involving (2.7) exists in the Riemann sense.

Now, we show that the solution defined by the formula (2.4) satisfies the integral equality (2.7); that is, 𝑢(𝑥,𝑡) is the weak solution of the problem (2.1)–(2.3).

According to the definition of the weak solution, we can write 0=^𝑇₀₀^𝐷𝑡𝑢()()𝑥,𝑡𝜕𝑓𝑥,𝑡()()𝜕𝑡+𝐷𝑢𝑥,𝑡𝜕𝑓𝑥,𝑡=𝜕𝑥𝑑𝑥𝑑𝑡_{𝑇00𝐷𝑡}𝑢()()𝑥,𝑡𝜕𝑓𝑥,𝑡𝜕𝑡𝑑𝑥𝑑𝑡+𝐷^𝑇₀₀^𝐷𝑡𝑢()()𝑥,𝑡𝜕𝑓𝑥,𝑡𝜕𝑥𝑑𝑥𝑑𝑡. (2.8) Changing the order of integration in the first integral and then applying the integration by parts to the inner integrals of the first and second terms in the last expression with respect to 𝑡 and 𝑥, respectively, we have ₀^𝐷𝑇𝑇_𝑥/𝐷𝑢()()𝑥,𝑡𝜕𝑓𝑥,𝑡𝜕𝑡𝑑𝑡𝑑𝑥+𝐷^𝑇₀₀^𝐷𝑡𝑢()()𝑥,𝑡𝜕𝑓𝑥,𝑡𝜕𝑥𝑑𝑥𝑑𝑡=0. (2.9) If we reapply integration by parts to the inner integrals in the first term and second term by 𝑡 and 𝑥, respectively, we prove that the solution in the form (2.4) satisfies the integral relation (2.7).

3. Auxiliary Problem and Its Exact Solution

In order to find the weak solution of the problem (2.1)–(2.3), according to [9, 10] the special auxiliary problem 𝜕𝑣=𝜕𝜕𝑡𝜕𝑥𝜕𝑣𝜕𝑥^𝜎𝑣(),𝑥,0=𝑣₀(𝑥)𝑢()=(),0,𝑡𝜕𝑣0,𝑡𝜕𝑥=𝜇₀𝑡^𝑛 (3.1)(3.2)(3.3) is introduced. Here, the function 𝑣₀(𝑥) is any solution of the equation 𝑑𝑣₀(𝑥)𝑑𝑥=0. (3.4) The problem (3.1)–(3.3) has the solution 𝑣(𝑥,𝑡)=−𝐷^1/(𝜎−1)𝜎−1𝜎^{𝜎/(𝜎−1)}(𝐷𝑡−𝑥)^{𝜎/(𝜎−1)},𝑥<𝐷𝑡,0,𝑥≥𝐷𝑡 (3.5) in the traveling wave form [9]. As is seen from (3.5), the differentiability property of the function 𝑣(𝑥,𝑡) is more than the differentiability property of the solution 𝑢(𝑥,𝑡). In addition to this, from (3.5) we get𝑢()=()𝑥,𝑡𝜕𝑣𝑥,𝑡.𝜕𝑥 (3.6)

Theorem 3.1. If the function 𝑣(𝑥,𝑡) is a soft solution of the problem (3.1)–(3.3), then the function 𝑢(𝑥,𝑡) obtained by (3.6) is a weak solution of the main problem (2.1)–(2.3) in sense of (2.7).

The auxiliary problem has the following advantages.(i)The function 𝑣(𝑥,𝑡) is smoother than 𝑢(𝑥,𝑡).(ii)The function 𝑣(𝑥,𝑡) is an absolutely continuous function.(iii)In the process of finding the solution 𝑢(𝑥,𝑡), one does not need to use the first and second derivatives of 𝑢(𝑥,𝑡) with respect to 𝑥.

The graphs of the function structured by the formulas (2.4), (3.5), and (3.6) are shown in the Figures 1(a), 1(b), and 2(a), respectively.

(a)

(b)

Figure 1 

(a) The exact solution of the main problem. (b) The exact solution of the auxiliary problem.

(a)

(b)

Figure 2 

(a) The function 𝑢(𝑥,𝑡)=𝜕𝑣(𝑥,𝑡)/𝜕𝑥. (b) Numerical solution obtained by using the classical algorithm (4.6)–(4.7).

4. Developing a Numerical Algorithm in a Class of Discontinuous Functions

In this section we investigate an algorithm for finding a numerical solution of the problem (2.1)–(2.3). At first, we cover the region 𝐷_ℓ(𝑡) by a special grid: 𝜔_𝜏=𝑥_𝑖,𝑡_𝑘∣𝑥_𝑖=𝐷𝑡_𝑖,𝑡_𝑘=𝑘𝜏,𝑖,𝑘=0,1,2,…, (4.1) where 𝜏 is the step of the grid with respect to 𝑡 variable.

Now, we will develop a numerical algorithm as follows. Since the function 𝜕𝑢^𝜎/𝜕𝑥 is continuous, we can approximate the problem (3.1)–(3.3) by the following finite differences schemes: 𝑉_𝑖,𝑘+1=𝑉_𝑖,𝑘+𝜏ℎ^𝜎+1𝑉_{𝑖+1,𝑘+1}−𝑉_𝑖,𝑘+1^𝜎−𝑉_𝑖,𝑘+1−𝑉_{𝑖−1,𝑘+1}^𝜎,𝑉_𝑖,0=𝑣₀𝑥_𝑖𝑉_0,𝑘+1,=𝑉_1,𝑘+1−ℎ𝜇₀𝑡^𝑛_𝑘+1 (4.2)(4.3)(4.4)(𝑖=0,1,2,…;𝑘=0,1,2,…). Here, ℎ=𝐷𝜏 and 𝑣₀(𝑥_𝑖) is any grid function of (3.4).

It can be easily shown that 𝑈_𝑖,𝑘+1=𝑉_{𝑖+1,𝑘+1}−𝑉_𝑖,𝑘+1ℎ, (4.5) and the grid function 𝑈_𝑖,𝑘+1 defined by (4.5) is a solution of the nonlinear system of algebraic equations: 𝑈_𝑖,𝑘+1=𝑈_𝑖,𝑘+𝜏ℎ²𝑈^𝜎_{𝑖+1,𝑘+1}−2𝑈^𝜎_𝑖,𝑘+1+𝑈^𝜎_{𝑖−1,𝑘+1} (4.6)(𝑖=1,2,…;𝑘=0,1,2,…). The initial and boundary conditions for (4.6) are 𝑈_𝑖,0(),𝑈=0𝑖=0,1,2,…_0,𝑘+1=𝜇₀𝑡^𝑛_𝑘+1(𝑘=0,1,2,…). (4.7) Here, 𝑉_𝑖,𝑘 and 𝑈_𝑖,𝑘 are the approximation values of 𝑣(𝑥,𝑡) and 𝑢(𝑥,𝑡) at any point (𝑥_𝑖,𝑡_𝑘) of the grid 𝜔_𝜏. For the sake of simplicity we introduce the notations 𝑉_𝑖,𝑘=𝑉, 𝑉_𝑖,𝑘+1=𝑉, and 𝑉_{𝑖±1,𝑘+1}=𝑉_±.

4.1. Convergence

In this section we will investigate some properties of the numerical solution and of the question of convergence of the numerical solution to the weak exact solution. Suppose that 𝜀_𝑖,𝑘,𝛿_𝑖,𝑘, and 𝜂_𝑖,𝑘 are the errors of approximations of the functions 𝜕𝑣/𝜕𝑥,𝜕𝑣/𝜕𝑡, and 𝜕𝜑/𝜕𝑥 by finite differences, respectively. Then, we can write (3.1) in the following form: 𝑣_𝑡+𝛿_𝑖,𝑘+1=𝜑_𝑥𝑣_𝑥+𝜀_𝑖,𝑘+1+𝜂_𝑖,𝑘+1 (4.8) or 𝑣_𝑡=𝜑_𝑥𝑣_𝑥+𝛾_𝑖,𝑘+1, (4.9) where 𝛾_𝑖,𝑘+1=𝛿_𝑖,𝑘+1+𝜑_𝑥𝑥𝑦_𝑥𝜀_𝑖,𝑘+1+𝜂_𝑖,𝑘+1. (4.10)

At first, we show that the finite difference scheme (4.2) is consistent; that is, 𝛾_𝑖,𝑘+1 approaches zero when 𝜏→0. It is known that the suitable characteristic of continuity of any function 𝑓(𝑥) on the any interval [𝑎,𝑏] is its modulus continuity:𝜔𝛿,𝑓≡𝜋𝑓=sup_{||𝑡−𝑥<𝛿}𝑓(𝑡)(𝑥,𝜀−𝑓)_𝑖,𝑘+1=𝜕𝑣𝑥_𝑖,𝑡_𝑘+1𝜕𝑥−𝑣_𝑥=𝜕𝑣𝑥_𝑖,𝑡_𝑘+1−𝑥𝜕𝑥𝜕𝑣^∗_𝑖,𝑡_𝑘+1𝑥𝜕𝑥=𝑢_𝑖,𝑡_𝑘+1−𝑢𝑥^∗_𝑖,𝑡_𝑘+1=𝜋(𝑢)⟶0,𝑥^∗_𝑖∈𝑥_𝑖−ℎ,𝑥_𝑖,𝛿_𝑖,𝑘+1=𝜕𝑣𝑥_𝑖,𝑡_𝑘+1𝜕𝑡−𝑣_𝑡=𝜕𝑣𝑥_𝑖,𝑡_𝑘+1−𝑥𝜕𝑡𝜕𝑣_𝑖,𝑡^∗_𝑘+1=𝑢𝑥𝜕𝑡𝜕𝜑_𝑖,𝑡_𝑘+1−𝑢𝑥𝜕𝑥𝜕𝜑_𝑖,𝑡^∗_𝑘+1(𝑢)𝜕𝑥=𝜋𝜕𝜑𝜕𝑥⟶0,𝑡^∗_𝑘+1∈𝑡_𝑘+1,𝑡_𝑘+1+𝜏. (4.11) Finally, 𝜂_𝑖,𝑘+1=𝜕𝜑(𝑢)𝜕𝑥−𝜑_𝑥()=𝑢𝑥𝑢𝜕𝜑_𝑖,𝑡_𝑘+1−𝑢𝑥𝜕𝑥𝜕𝜑_𝑖^∗∗,𝑡_𝑘+1𝜕𝑥⟶0,𝑥_𝑖^∗∗∈𝑥_𝑖−ℎ,𝑥_𝑖. (4.12) Therefore, 𝛾_𝑖,𝑘+1→0, if 𝜏→0.

Theorem 4.1 (Maximum Principle). The solution of problem (4.2)–(4.4) takes its maximum (or minimum) value on the boundary of the domain of definition of the solution, that is, 0≤𝑉_𝑖,𝑘+1≤𝑀=max_𝑖,𝑘𝑢₀,𝑢₁. (4.13)

Proof. At first, let us write (4.2) in the following form: 𝑉_𝑖,𝑘+1−𝑉_𝑖,𝑘𝜏=1ℎ𝑉_{𝑖+1,𝑘+1}−𝑉_𝑖,𝑘+1ℎ^𝜎−𝑉_{𝑖,𝑘+1+1}−𝑉_{𝑖−1,𝑘+1}ℎ^𝜎=1ℎ𝑈^𝜎_𝑖−𝑈^𝜎_𝑖−1=−𝑤𝑥_𝑖,𝑡_𝑘. (4.14)
Assume that 𝑉_𝑖,𝑘+1 is not constant and 𝑉_𝑖,𝑘+1 takes the greatest value at some point of the grid 𝜔_𝜏 rather than at boundary nodes of 𝛾_𝜏. Then, there is such a point (𝑥₁,𝑡₁)∈𝜔_𝜏 that 𝑉 takes the maximal value and even some neighborhood points 𝑉(𝑥₁,𝑡₁) less than 𝑉(𝑥₁,𝑡₁).
If 𝑉_𝑥1,𝑡1>𝑉_𝑥1,𝑡1, because the function 𝜑(𝑢)=𝑢^𝜎 is monotone, the left part of the relation (4.14) is positive, but the right part is negative. Hence we arrive to inconsistency. We arrive to the same inconsistency if 𝑉_±(𝑥₁,𝑡₁)>𝑉(𝑥₁,𝑡₁). Similarly, we can prove that 𝑉 does not take a minimal value at the inner nodes of the grid 𝜔_𝜏.

From Theorem 4.1, it follows that the solution of the problem (4.2)–(4.4) converges to the solution of problem (3.1)–(3.3) pointwise, that is, max_𝑖𝑣_𝑖,𝑘−𝑉_𝑖,𝑘⟶0. (4.15) It can be easily seen that the solution of the mentioned problem is continuously dependent on initial data.

Now, we will prove convergence of the 𝑈_𝑖,𝑘 to the solution of the main problem. To this end, by subtracting (4.2) from (4.9) and taking (4.5) into account, we have 𝑅_𝑡=𝜑^′(𝑢)𝑣_𝑥−𝑉_𝑥𝑥+𝛾_𝑖,𝑘 (4.16) or 𝑅_𝑡=𝜑^′(𝑢)𝑢_𝑖,𝑘−𝑈_{𝑖,𝑘𝑥}+𝛾_𝑖,𝑘,𝑢∈𝑢_𝑖,𝑘,𝑈_𝑖,𝑘, (4.17) where 𝑅=𝑣_𝑖,𝑘+1−𝑉_𝑖,𝑘+1. By multiplying the last equation with 𝑅 and summing with respect to 𝑖 and 𝑘, we get 𝑅,𝑅_{𝑡𝐿2(𝜔𝜏)}=𝑅,𝜑^′(𝑢)𝑢_𝑖,𝑘−𝑈_{𝑖,𝑘𝑥𝐿2(𝜔𝜏)}+𝑅,𝑣_{𝑖,𝑘𝐿2(𝜔𝜏)}. (4.18) Here, the symbol (𝑓,𝑔) denotes the differences analogy of the scalar product of the functions of 𝑓 and 𝑔 in 𝐿₂(𝜔_𝜏) sense: 𝑓,𝑔_{𝐿2(𝜔𝜏)}=𝜏ℎ_𝑖=1𝑘=0𝑓_𝑖,𝑘𝑔_𝑖,𝑘. (4.19) With some algebra we have 12_𝑖𝑅²_𝑖𝑡𝑘=𝑇+𝑢_𝑖,𝑘−𝑈_𝑖,𝑘,𝜑^′(𝑢)𝑢_𝑖,𝑘−𝑈_{𝑖,𝑘𝐿2(𝜔𝜏)}≤12_𝑖𝑅²_{𝑖𝑡𝑘=0}+𝑇max_𝑖𝑅²_𝑖𝛾_{𝑖,𝑘𝐿2(𝜔𝜏)}. (4.20) It follows that the numerical solution 𝑈_𝑖,𝑘 converges in mean to 𝑢_𝑖,𝑘.

5. Numerical Experiments

In order to extract the significance of the suggested method, the numerical solution obtained using the proposed auxiliary problem is compared with the exact solution of the problem (2.1)–(2.3) on an equal footing. With this aim, firstly, computer experiments are performed on the algorithm (4.2)–(4.4) with the values 𝜎=3, 𝑇=2,4, and 6. The graphs of the obtained numerical solutions are presented in Figure 3(a). The numerical solutions of the main problem (2.1)–(2.3) obtained by using the algorithm (4.2)–(4.4) are shown in Figure 3(b). The numerical solution of the main problem obtained using the algorithm (4.6)–(4.7) is demonstrated in Figure 2(b).

(a)

(b)

Figure 3 

(a) Numerical solution of the problem (4.2)–(4.4). (b) Numerical solution of the main problem obtained by using the solution of the problem (4.2)–(4.4).

Comparing Figures 1(a) and 3(b) shows that the numerical and the exact solutions of the main problem (2.1)–(2.3) coincide. Moreover, the graphs of the functions 𝑣(𝑥,𝑡) and 𝑉_𝑖,𝑘 coincide, too.

Thus, the numerical experiments carried out show that the suggested numerical algorithms are efficient and economical from a computer point of view. The proposed algorithms permit us to develop the higher-resolution methods where the obtained solution correctly describes all physical features of the problem, even if the differentiability order of the solution of the problem is less than the order of differentiability which is required by the equation from the solution.

The finite differences scheme (4.2)–(4.4) has the first-order approximation with respect to 𝑡. But using the different versions of Runge-Kutta's methods we can increase the order of the algorithms mentioned previously.

6. Conclusions

The new method is suggested for obtaining the regular weak solution for the free boundary problem of the nonlinear parabolic type equation.

The auxiliary problem which has some advantages over the main problem is introduced and it permits us to find the exact solution with singular properties.

The auxiliary problem introduced previously allows us to develop the higher-resolution method where the obtained solution correctly describes all physical features of the problem, even if the differentiability order of the solution of the problem is less than the order of differentiability which is required by the equation from the solution.