# A Method for Generating Approximate Similarity Solutions of Nonlinear Partial Differential Equations

Authors: Mazhar Iqbal; M. T. Mustafa; Azad A. Siddiqui
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

Standard application of similarity method to find solutions of PDEs mostly results in reduction to ODEs which are not easily integrable in terms of elementary or tabulated functions. Such situations usually demand solving reduced ODEs numerically. However, there are no systematic procedures available to utilize these numerical solutions of reduced ODE to obtain the solution of original PDE. A practical and tractable approach is proposed to deal with such situations and is applied to obtain approximate similarity solutions to different cases of an initial-boundary value problem of unsteady gas flow through a semi-infinite porous medium.


---

## Body

## 1. Introduction

The mathematical modeling of most of the physical processes in fields like diffusion, chemical kinetics, fluid mechanics, wave mechanics, and general transport problems is governed by such nonlinear PDEs whose analytic solutions are hard to find. Therefore, the approach of investigating nonlinear PDEs via reduction to ODEs becomes important and has been quite fruitful in analysis of many physical problems. The reader is referred to [1] for an introduction to different types of such reduction approaches and to have an idea about the advances made in the fields of nonlinear diffusion, fluid mechanics, and wave propagation from the utilization of reduction-to-ODE approach.Lie symmetry method provides a powerful general technique for analyzing nonlinear PDEs and reducing them to ODEs. PDEs modeling physical processes naturally inherit symmetries from the underlying physical system. Lie symmetry method takes advantage of these natural symmetries in a PDE and provides similarity variables that lead to the reduction to ODEs; compare, for example, [2–4]. A large amount of the literature about the Lie symmetry theory, its applications, and its extensions is available; see, for example, [1, 5–12].Although there have been some notable contributions in the applications of similarity method to boundary value problems of PDEs (cf. [1, 12–15]), in general the method has not been utilized in a great deal mainly because of a restriction explained below. The success of the method relies heavily on the success in solving the reduced ODE. Only in few cases, the reduced ODE is integrable in terms of elementary or tabulated functions, but in most cases it is not a simple matter and it is suggested to resort to numerical methods to solve the reduced ODE. This retreat from analytic calculations to numerical computations is not uncommon, but it has not proved practical due to the lack of practical and systematic procedures to utilize the numerical solution of reduced ODE to obtain the solution of original PDE. The aim of this work is to provide one practical way of dealing with such situations.As a test problem to demonstrate our method, we consider a physical problem arising in the transient flow of gas through a semi-infinite porous medium. In the investigation of the unsteady flow of gas through a semi-infinite porous medium [16–19] initially filled with gas at a uniform pressure p0>0 at time t=0, the pressure at the outflow face is suddenly reduced from p0 to p1≥0 (p1=0 is the case of diffusion into a vacuum) and is thereafter maintained at this lower pressure. The studies in [18, 20] show that the unsteady flow of gas in a porous medium is modeled by a nonlinear partial differential equation.The nonlinear partial differential equation that describes the unsteady flow of gas through a semi-infinite porous medium has been derived by Muskat [20] in the form
(1)∇2(p2)=2(ϕμk)∂p∂t,
where p is the pressure within porous medium, ϕ is the porosity, μ is the viscosity, and k is the permeability. In the one-dimensional medium extending from x=0 to x=∞, this reduces to
(2)∂∂x(p∂p∂x)=A∂p∂t,
with the boundary conditions
(3)p(x,0)=p0,0<x<∞p(0,t)=p1(<p0),0≤t<∞,
where the constant A=ϕμ/k is given by the properties of the medium.In the next section we provide details about the method and its application to obtain approximate similarity solutions of BVPs of the above form.


## 2. The Method and Its Implementation

The main idea of the method rests on finding an approximation of solution of reduced ODE in the form of a function. This has clear advantages over the option of working only with numerical solution of reduced ODE. In the first place, it can be used to generate approximate solution of original PDE via the similarity variables. Secondly, having the approximate solution of PDE in function form can be more useful than numerical solution as it displays the parameters and variables of the problem, so it requires less processing time and can be used for applications in real time. A detailed description of the method is provided below.To explain and illustrate our method we consider the IBVP(4)∂∂x(p∂p∂x)=A∂p∂t,
with
(5)p(x,0)=p0,p(0,t)=p1(<p0),p(∞,t)=p0.
Without loss of generality we can assume p0=1 because the change of variable p¯(x,t)=p(x,t)/p0 leads to a similar problem with p¯(x,0)=1. So the IBVP under study here is
(6)∂∂x(p∂p∂x)=A∂p∂t,
with initial and boundary conditions
(7)p(x,0)=1,0<x<∞,p(0,t)=p1(<1),0≤t<∞,p(∞,t)=1,0<t<∞.Step 1 (Reduction of IVBP to a BVP of ODE).
The similarity transformations [16, 17]
(8)z=xt(A4)1/2,w(z)=α-1(1-(p(x,t))2),
with α=1-p12, reduce the above IBVP (6)-(7) to BVP of ODE as follows:
(9)w′′+2z1-αww′=0,w(z=0)=1,w(z⟶∞)=0.
The aim of the remaining steps is to find an approximate solution WApprox(z) of BVP (9) in function form and then use the similarity transformations (8) to obtain approximate p(x,t) of PDE problem (6)-(7).Step 2.
Find numerical solutionWNum of BVP (9) and use this as a benchmark for obtaining function form WApprox of the solution.Step 3 (Obtain an initial guess forWApprox).
This is a crucial step and in our case we use the lower solution of BVP (9), obtained in [17], as our initial guess. As shown in numerical simulations below, lower solution provides a good initial guess that leads to an accurate enough approximate solution in few iterations. So the initial guess for approximate solution of BVP (9) for all the cases below is taken as follows, see [17, Example 2.1]:
(10)initialapproximation=WLower=1-erf(zp1),
where erf denotes the error function.
As an alternate, a solution obtained by homotopy analysis method can also be used as initial guess.Step 4.
Improve the initial approximation to getWApprox up to the desired level of accuracy.
Here we adopt the following procedure for improving the level of accuracy, starting from initial approximation. The lower solution is of the form(11)1-erf(kz),
with k=k0=1/p1 giving the initial approximation Wk0, that is, the lower solution. Numerical simulations suggest that as the values of k decrease from k0 by a small decrement, the lower solution moves uniformly towards the numerical solution.
Given a functionN(x) and a number ɛ>0, we say that f(x) lies within ɛ-band of N(x) on an interval I if
(12)|f(x)-N(x)|<ɛ,∀x∈I.
For a suitable value n and numbers ɛ>0, δi>0, using the sequences of values
(13)k=ki=k0-δi,(i=1,2,…,n),
in (11) generates a sequence of curves Wki that uniformly approach numerical solution, finally resulting in the curve
(14)WApprox=Wkn,
which lies in an ɛ-band around the graph of numerical solution WNum. The number ɛ is chosen according to the desired level of accuracy and the value of kn is approximated via numerical simulations.Step 5.
Use the similarity variables (8) to get the approximate solution of original IBVP (6)-(7).
In the subsequent subsections we illustrate implementation of the above procedure, and we provide simulation results and approximate similarity solutions for different cases of values of parameterp1.


### 2.1. Approximate Similarity Solution of the IBVP (6)-(7) for p1=0.9

In this case,k=1/0.9 and the lower solution of the ODE problem (9) becomes
(15)WLower=1-erf(kz)=1-erf(z0.9)≈1-erf(1.05409z),
which serves as our initial approximation Wk0 for approximating the solution of ODE problem. Solving the BVP (9) numerically to get WNum and uniformly improving the approximations Wki by simulating the procedure explained above we obtain an approximate solution of the ODE problem as
(16)WApprox(z)=1-erf(0.7281z),
with
(17)Max⁡|WApprox-WNum|=0.0019642487.The plots of the initial approximationWLower(z), the numerical solution WNum(z), the approximate solution WApprox(z), and the Error(z) are given in Figure 1 where
(18)Error(z)=WApprox(z)-WNum(z).
The dotted curves demonstrate some intermediary curves involved in simulations of the uniform approximation process from WLower to WApprox, in a manner that as k decreases from 1/0.9 by a small decrement the curves 1-erf(kz) move uniformly towards numerical solution WNum.Figure 1Finally the similarity variables (8) provide the approximate solution of IBVP (6)-(7) for p1=0.9 as
(19)p(x,t)=0.81+0.19erf(0.7281xt(A4)1/2).


### 2.2. Approximate Similarity Solution of the IBVP (6)-(7) for p1=0.3

In this case,k=1/0.3 and hence the initial approximation determined by the lower solution of the ODE problem (9) becomes
(20)WLower=1-erf(kz)=1-erf(z0.3)≈1-erf(1.82574z).
Proceeding as in Section 2.1 we obtain an approximate solution of the ODE problem, for the case p1=0.3, as
(21)WApprox(z)=1-erf(0.77362z),
with
(22)Max⁡|WApprox-WNum|=0.0096829849.
The plots of the initial approximation WLower(z), the numerical solution WNum(z), the approximate solution WApprox(z), and the Error(z), for the case p1=0.3, are given in Figure 2.Figure 2As in Section2.1 the dotted curves demonstrate some intermediary curves involved in simulations of the uniform approximation process from WLower to WApprox, in a manner that as k decreases from 1/0.3 by a small decrement the curves 1-erf(kz) move uniformly towards numerical solution WNum.The approximate solution of IBVP (6)-(7) for p1=0.3 can be found using the similarity variables (8) and is given by
(23)p(x,t)=0.09+0.91erf(0.77362xt(A4)1/2).


### 2.3. Approximate Similarity Solution of the IBVP (6)-(7) for p1=0.1

Usingk=1/0.1 gives the lower solution of the ODE problem (9) as
(24)WLower=1-erf(kz)=1-erf(z0.1)≈1-erf(3.162277z),
which serves as our initial approximation Wk0 for approximating the solution of ODE problem for this case. As above, solving the BVP (9) numerically to get WNum and uniformly improving the approximations Wki by simulating the procedure explained above, we obtain an approximate solution of the ODE problem as
(25)WApprox(z)=1-erf(0.782743z),
with
(26)Max⁡|WApprox-WNum|=0.0116428185.The plots of the initial approximationWLower(z), the numerical solution WNum(z), the approximate solution WApprox(z), and the Error(z), for the case p1=0.1, are given in Figure 3. As in the above sections, the dotted curves demonstrate some intermediary curves involved in simulations of the uniform approximation process from WLower to WApprox, in a manner that as k decreases from 1/0.1 by a small decrement the curves 1-erf(kz) move uniformly towards numerical solution WNum.Figure 3Using the similarity variables, the approximate solution of IBVP (6)-(7) for p1=0.1 is found as
(27)p(x,t)=0.01+0.99erf(0.782743xt(A4)1/2).


### 2.4. Discussion

The plots of the numerical solutionWNum(z), the lower solution WLower(z), and the approximate solution WApprox(z) for the ODE problem (9) are given in Figures 1, 2, and 3, respectively, for p1=0.9, p1=0.3, and p1=0.1. It can be seen that the lower solution WLower(z) provides a better initial approximation for the values of p1 closer to 1. Yet a comparison of the plots of initial approximation WLower(z) and the accuracy benchmark curve WNum(z), in all cases, clearly emphasizes the need of improving the initial approximation to get WLower(z) to get accurate enough analytic approximation WApprox(z) of the benchmark curve WNum(z). The plots in Figures 1, 2, and 3 clearly show that WApprox(z) and WNum(z) curves overlap in each case and the corresponding plots of the error Error(z) demonstrate that, in each case, WApprox(z) provides an analytic approximation of the benchmark numerical solution curve WNum(z). In all cases, the solution WApprox(z) decreases with the increase in z but descends faster for the case p1=0.1 as compared to case p1=0.9.A comparison of the approximate solutionWApprox(z) of BVP (9) obtained by our method with solutions obtained by other methods [21–23], for α=0.5, is presented in Table 1.Table 1
Solution values forα=0.5.
z
W
Approx
(
z
)
Corresponding solution
Corresponding solution
W
kidder [18]
W
[
3
/
3
]  [21–23]
0.1
0.899431
0.881659
0.897917
0.2
0.800454
0.766308
0.798523
0.3
0.704584
0.656538
0.704113
0.4
0.613194
0.554402
0.616504
0.5
0.527452
0.461365
0.537053
0.6
0.448284
0.378311
0.466563
0.7
0.37634
0.305598
0.406243
0.8
0.311998
0.243133
0.35608
0.9
0.255363
0.190462
0.317997
1.0
0.206302
0.158769
0.290026


### 2.5. Comments on the Conservation Laws

The conservation laws of(28)∂∂x(p∂p∂x)=∂p∂t
can be found using the direct multiplier method [24, 25]. Since the procedure is standard, only the results are stated. Implementing the direct multiplier procedure for multipliers of the form Λ(t,x,p,px,pxx) generates the determining equations
(29)Λxx=0,Λt=0,Λp=0,Λpx=0,Λpxx=0,
which yield the solution Λ=C1+C2x. Hence there are two linearly independent conservation laws arising from the multipliers Λ1=1 and Λ2=x. Next, elementary direct calculations [26] give the corresponding conservation laws as listed below.ForΛ1=1, the conservation law is
(30)Dt(p)+Dx(-ppx)=0,
and, for Λ2=x, the conservation law is
(31)Dt(xp)+Dx(p22-xppx)=0,
where
(32)Dt=∂∂t+pt∂∂p+ptt∂∂pt+ptx∂∂px+⋯,Dx=∂∂x+px∂∂p+pxt∂∂pt+pxx∂∂px+⋯.
For a given conservation law
(33)Dt(T1)+Dx(T2)=0,
of (28), if the spatial flux T2 vanishes on the boundary x=0 and x=∞ of the semi-infinite medium, then integration from x=0 to x=∞ provides conserved quantity of the boundary value problem.As an example, the conservation law (30) is applied to derive conserved quantity of the BVP
(34)∂∂x(p∂p∂x)=∂p∂t,∂p∂x(0,t)=0,∂p∂x⟶ke-x2/t,asx⟶∞.
The conservation law (30) can be written as
(35)∂∂t(p)-∂∂x(ppx)=0.
Integrating with respect to x from x=0 to x=∞ and using the boundary conditions imply that (36)∂∂t∫0∞p(x,t)dx=0,
which gives the time independent conserved quantity ∫0∞‍p(x,t)dx of the BVP.In general, the boundary conditions will determine which conservation law is to be applied to obtain conserved quantities of the BVP of (28).


## 2.1. Approximate Similarity Solution of the IBVP (6)-(7) for p1=0.9

In this case,k=1/0.9 and the lower solution of the ODE problem (9) becomes
(15)WLower=1-erf(kz)=1-erf(z0.9)≈1-erf(1.05409z),
which serves as our initial approximation Wk0 for approximating the solution of ODE problem. Solving the BVP (9) numerically to get WNum and uniformly improving the approximations Wki by simulating the procedure explained above we obtain an approximate solution of the ODE problem as
(16)WApprox(z)=1-erf(0.7281z),
with
(17)Max⁡|WApprox-WNum|=0.0019642487.The plots of the initial approximationWLower(z), the numerical solution WNum(z), the approximate solution WApprox(z), and the Error(z) are given in Figure 1 where
(18)Error(z)=WApprox(z)-WNum(z).
The dotted curves demonstrate some intermediary curves involved in simulations of the uniform approximation process from WLower to WApprox, in a manner that as k decreases from 1/0.9 by a small decrement the curves 1-erf(kz) move uniformly towards numerical solution WNum.Figure 1Finally the similarity variables (8) provide the approximate solution of IBVP (6)-(7) for p1=0.9 as
(19)p(x,t)=0.81+0.19erf(0.7281xt(A4)1/2).


## 2.2. Approximate Similarity Solution of the IBVP (6)-(7) for p1=0.3

In this case,k=1/0.3 and hence the initial approximation determined by the lower solution of the ODE problem (9) becomes
(20)WLower=1-erf(kz)=1-erf(z0.3)≈1-erf(1.82574z).
Proceeding as in Section 2.1 we obtain an approximate solution of the ODE problem, for the case p1=0.3, as
(21)WApprox(z)=1-erf(0.77362z),
with
(22)Max⁡|WApprox-WNum|=0.0096829849.
The plots of the initial approximation WLower(z), the numerical solution WNum(z), the approximate solution WApprox(z), and the Error(z), for the case p1=0.3, are given in Figure 2.Figure 2As in Section2.1 the dotted curves demonstrate some intermediary curves involved in simulations of the uniform approximation process from WLower to WApprox, in a manner that as k decreases from 1/0.3 by a small decrement the curves 1-erf(kz) move uniformly towards numerical solution WNum.The approximate solution of IBVP (6)-(7) for p1=0.3 can be found using the similarity variables (8) and is given by
(23)p(x,t)=0.09+0.91erf(0.77362xt(A4)1/2).


## 2.3. Approximate Similarity Solution of the IBVP (6)-(7) for p1=0.1

Usingk=1/0.1 gives the lower solution of the ODE problem (9) as
(24)WLower=1-erf(kz)=1-erf(z0.1)≈1-erf(3.162277z),
which serves as our initial approximation Wk0 for approximating the solution of ODE problem for this case. As above, solving the BVP (9) numerically to get WNum and uniformly improving the approximations Wki by simulating the procedure explained above, we obtain an approximate solution of the ODE problem as
(25)WApprox(z)=1-erf(0.782743z),
with
(26)Max⁡|WApprox-WNum|=0.0116428185.The plots of the initial approximationWLower(z), the numerical solution WNum(z), the approximate solution WApprox(z), and the Error(z), for the case p1=0.1, are given in Figure 3. As in the above sections, the dotted curves demonstrate some intermediary curves involved in simulations of the uniform approximation process from WLower to WApprox, in a manner that as k decreases from 1/0.1 by a small decrement the curves 1-erf(kz) move uniformly towards numerical solution WNum.Figure 3Using the similarity variables, the approximate solution of IBVP (6)-(7) for p1=0.1 is found as
(27)p(x,t)=0.01+0.99erf(0.782743xt(A4)1/2).


## 2.4. Discussion

The plots of the numerical solutionWNum(z), the lower solution WLower(z), and the approximate solution WApprox(z) for the ODE problem (9) are given in Figures 1, 2, and 3, respectively, for p1=0.9, p1=0.3, and p1=0.1. It can be seen that the lower solution WLower(z) provides a better initial approximation for the values of p1 closer to 1. Yet a comparison of the plots of initial approximation WLower(z) and the accuracy benchmark curve WNum(z), in all cases, clearly emphasizes the need of improving the initial approximation to get WLower(z) to get accurate enough analytic approximation WApprox(z) of the benchmark curve WNum(z). The plots in Figures 1, 2, and 3 clearly show that WApprox(z) and WNum(z) curves overlap in each case and the corresponding plots of the error Error(z) demonstrate that, in each case, WApprox(z) provides an analytic approximation of the benchmark numerical solution curve WNum(z). In all cases, the solution WApprox(z) decreases with the increase in z but descends faster for the case p1=0.1 as compared to case p1=0.9.A comparison of the approximate solutionWApprox(z) of BVP (9) obtained by our method with solutions obtained by other methods [21–23], for α=0.5, is presented in Table 1.Table 1
Solution values forα=0.5.
z
W
Approx
(
z
)
Corresponding solution
Corresponding solution
W
kidder [18]
W
[
3
/
3
]  [21–23]
0.1
0.899431
0.881659
0.897917
0.2
0.800454
0.766308
0.798523
0.3
0.704584
0.656538
0.704113
0.4
0.613194
0.554402
0.616504
0.5
0.527452
0.461365
0.537053
0.6
0.448284
0.378311
0.466563
0.7
0.37634
0.305598
0.406243
0.8
0.311998
0.243133
0.35608
0.9
0.255363
0.190462
0.317997
1.0
0.206302
0.158769
0.290026


## 2.5. Comments on the Conservation Laws

The conservation laws of(28)∂∂x(p∂p∂x)=∂p∂t
can be found using the direct multiplier method [24, 25]. Since the procedure is standard, only the results are stated. Implementing the direct multiplier procedure for multipliers of the form Λ(t,x,p,px,pxx) generates the determining equations
(29)Λxx=0,Λt=0,Λp=0,Λpx=0,Λpxx=0,
which yield the solution Λ=C1+C2x. Hence there are two linearly independent conservation laws arising from the multipliers Λ1=1 and Λ2=x. Next, elementary direct calculations [26] give the corresponding conservation laws as listed below.ForΛ1=1, the conservation law is
(30)Dt(p)+Dx(-ppx)=0,
and, for Λ2=x, the conservation law is
(31)Dt(xp)+Dx(p22-xppx)=0,
where
(32)Dt=∂∂t+pt∂∂p+ptt∂∂pt+ptx∂∂px+⋯,Dx=∂∂x+px∂∂p+pxt∂∂pt+pxx∂∂px+⋯.
For a given conservation law
(33)Dt(T1)+Dx(T2)=0,
of (28), if the spatial flux T2 vanishes on the boundary x=0 and x=∞ of the semi-infinite medium, then integration from x=0 to x=∞ provides conserved quantity of the boundary value problem.As an example, the conservation law (30) is applied to derive conserved quantity of the BVP
(34)∂∂x(p∂p∂x)=∂p∂t,∂p∂x(0,t)=0,∂p∂x⟶ke-x2/t,asx⟶∞.
The conservation law (30) can be written as
(35)∂∂t(p)-∂∂x(ppx)=0.
Integrating with respect to x from x=0 to x=∞ and using the boundary conditions imply that (36)∂∂t∫0∞p(x,t)dx=0,
which gives the time independent conserved quantity ∫0∞‍p(x,t)dx of the BVP.In general, the boundary conditions will determine which conservation law is to be applied to obtain conserved quantities of the BVP of (28).


## 3. Conclusion

We present a practical way of obtaining approximate solution, in function form, for the class of PDEs where the PDE can be reduced to an ODE through similarity variables but the reduced ODE is not easily integrable in terms of elementary or tabulated functions. The idea presented here to get approximate solution of PDE, that is, approximating the surfacep(x,t), practically involves approximating a curve W(z) which is a tractable problem in comparison to increasingly complex and intractable problem of approximating the surface p(x,t) itself. A combination of simulations, initial approximation, and numerical solution of reduced ODE is utilized to obtain approximation of solution curve W(z) which readily generates, via similarity variables, the approximate solution surface p(x,t) of the PDE. This makes it a promising approach especially when reasonably accurate initial approximation of the solution of ODE can be obtained, as was the case here in the form of lower solution. The approach is applied to obtain approximate solutions for some cases of an initial-boundary value problem of unsteady flow of gas through a semi-infinite porous medium. The approach can be adapted for obtaining approximate analytic solutions for the class of PDEs where the PDE can be reduced to an ODE through similarity variables. For instance, the approach can be directly applied to all the reduced-via-similarity BVPs of ODEs in [27]. For further application, the approach can be extended to obtain approximate solutions where the reduction is a system of ODEs like the reduced flow problems in [28, 29].



---

*Source: 105414-2014-07-22.xml*