Asymptotic behaviour of a solution to a nonlinear equation modelling capillary rise

TL;DR

This paper analyzes the asymptotic behavior of solutions to a nonlinear ODE modeling capillary rise, proving convergence and developing approximation methods despite non-Lipschitz nonlinearities.

Contribution

It introduces a novel analysis of a singular nonlinear ODE for capillary rise, including convergence proofs and asymptotic approximation techniques for non-Lipschitz nonlinearities.

Findings

Proves convergence of solutions as a parameter approaches zero

Develops an accurate asymptotic approximation method for large times

Addresses nonlinear ODEs with non-Lipschitz components

Abstract

We are concerned with the asymptotics and perturbation analysis of a singular second-order nonlinear ODE that models capillary rise of a fluid inside a narrow vertical tube. We prove the convergence of the exact solution to a unperturbed solution when a nondimensional parameter decrease to zero. Furthermore, we provide an accurate method of approximating the asymptotic solution for large times. Due to the a fact that the nonlinear component in the main equation does not satisfy the Lipschitz continuity condition the methods used to prove the main theorems are nonstandard, require careful analysis, and can be useful in dealing with similar nonlinear ODEs.

Equations

$$\frac{8\nu }{r^2}h{h}^{\prime }+\rho gh+\rho \frac{d}{dt}\left(h{h}^{\prime }\right)=\frac{2\gamma \cos \theta }{r},$$

$$u^{\prime \prime }(s)+\frac{1}{\sqrt{\omega }}{u}^{\prime }(s)+\sqrt{2u(s)}=1,$$

$$u(0)=0,u^{\prime }(0)=0,$$

$$\omega =\frac{{\rho }^3{g}^2{r}^5}{128{\mu }^2\gamma \cos \theta }.$$

$$u^{\prime \prime }(s)+ϵu^{\prime }(s)+\sqrt{2u(s)}=1,u(0)=0,u^{\prime }(0)=0.$$

$$u_0(s)=\frac{1}{72}\sum _{i=1}^{\infty }(s-6(i-1){)}^2(6i-s{)}^2\text{mathbbm}{1}_{\{s\in [6(i-1),6i]\}}.$$

$$u_0^{\prime \prime }+\sqrt{2u_0}=1,u_0(0)=0,u_0^{\prime }(0)=0.$$

$$g\frac{dg}{d{u}_0}+\sqrt{2u_0}=1,g(0)=0,$$

$$u_{0}^{\prime}(x)=\pm\sqrt{\frac{2}{3}}\sqrt{-\biggl{(}2\sqrt{2}\sqrt{u(x)}-3\biggr{)}u_{0}(x)},\quad u_{0}(0)=0.$$

$$u_0(s)=\frac{s^2(6-s{)}^2}{72}\text{mathbbm}{1}_{\{s\in [0,6]\}}.$$

$$\frac{u^{\prime 2}}{2}+\frac{2\sqrt{2}}{3}{u}^{3/2}-u=-ϵ{\int }_0^s{u}^{\prime }(x{)}^{2}dx$$

$$\frac{1}{3}\frac{s^2}{2}\le u(s)\le \frac{s^2}{2}.$$

$$u^{\prime 2}(s)\le s^2-\frac{2\sqrt{2}}{9}{s}^3.$$

$$u(s)=\mathcal{O}(s^2),\text{as}s\to 0^{+}.$$

$$u^{\prime }(s)=\mathcal{O}(s),\text{as}s\to 0^{+}.$$

$$\left(1-\frac{1}{3}(1+ϵ)s\right)\frac{s^2}{2}\le u(s)\le \frac{s^2}{2}\text{for}s\in \left[0,\frac{3}{1+ϵ}\right].$$

$$\frac{1}{3}=\left(1-\frac{1}{3}(1+ϵ)\left(\frac{2}{1+ϵ}\right)\right)\le \left(1-\frac{1}{3}(1+ϵ)s\right)\text{for}s\in \left[0,\frac{2}{(1+ϵ)}\right].$$

$$\frac{u^{\prime 2}}{2}+\frac{2\sqrt{2}}{3}{u}^{3/2}-u\le 0.$$

$$u^{\prime 2}(s)\le 2u-\frac{4\sqrt{2}}{3}{u}^{3/2}\le s^2-\frac{2\sqrt{2}}{9}{s}^3.$$

$$|u^{\prime }(s)|\le s,$$

$$|u^{\prime }(s)|\le \sqrt{s^2-\frac{2\sqrt{2}}{9}{s}^3}\le s,$$

$$\frac{1}{2}{u}^{\prime }(\delta {)}^2+\frac{2\sqrt{2}}{3}u(\delta {)}^{\frac{3}{2}}-u(\delta )=-ϵ\underset{0}{\overset{\delta }{\int }}{u}^{\prime }(t{)}^{2}dt.$$

$$-ϵ\underset{0}{\overset{\delta }{\int }}{u}^{\prime }(t{)}^{2}dt\ge -ϵ\delta u^{\prime }(\delta {)}^2.$$

$$u^{\prime }(\delta {)}^2\ge \frac{1}{3(\frac{1}{2}+ϵ\delta )}({\delta }^2-{\delta }^3),$$

$$u^{\prime }(\delta )\ge M_1\frac{1}{\sqrt{3(\frac{1}{2}+ϵ\delta )}}\delta ,$$

$$M_1\frac{1}{\sqrt{3(\frac{1}{2}+ϵ\delta )}}\delta \le u^{\prime }(\delta )\le \delta ,$$

$$\underset{ϵ\to 0^{+}}{\lim }|u(s)-u_0(s)|=0,$$

$$u(s)=u_0(s)+R(s),$$

$$R^{\prime \prime }(s)=\sqrt{2u_0(s)}-\sqrt{2(u_0(s)+R(s))}.$$

$$R(s)={\int }_0^{s}G(s-x)f(u_0(x),R(x))dx$$

Full text

Keywords: nonlinear oscillations, singular perturbations, asymptotic behaviour, non-Lipschitzian function, Washburn’s equation, capillary rise

1 Introduction

Capillary rise is an extraordinary physical phenomenon that is ubiquitous in nature and the effects of which we can notice in many everyday situations. This natural phenomenon is mainly responsible for water distribution in plants, in particular it causes the main driving force in the process of water supply to the upper branches of trees. Capillary rise can also be observed in porous media such as soils and rocks. In addition to its obvious natural applications, the discussed phenomenon appears also in industry (see [20] and [3]), e.g. in the printing processes and in a technique called thin layer chromatography to separate mixtures from substances. The appropriate mathematical models describing the behavior of capillary rise are gaining in importance in view of the applications of a capillary flow in industry.

One of the first equations modelling the dynamics of capillary rise was introduced by E. W. Washburn in 1921 (see [18]). However, after subsequent researches and conducted experiments by others scientists, many improved mathematical models were proposed in order to provide the more accurate description of that extraordinary physical phenomenon. In our previous work [13] we provided brief explanation of each component in the main governing equation modelling the capillary rise process in a narrow vertical tube

[TABLE]

where $h=h(t)$ is the liquid column height at time $t$, $r$ is a radius of the capillary tube, and $\mu,\,\rho$ denote viscosity and density of liquid respectively. As usual, $g$ denotes the gravitational acceleration and $\gamma$ denotes the surface tension. It is worth emphasis that in the current work the contact angle, here denoted by $\theta$, preserves a constant value during the flow. Nevertheless, some results [14, 2, 6] indicates the dependence of contact angle on the velocity of flow. The meticulous investigation of the (1) with various models of the velocity dependent contact angle will be an object of study in further works. Moreover, we assume that the narrow tube is initially filled with air and the walls of tube are smooth. For simplification and clarity of the main process, we do not take into account the influence of the container walls where the narrow tube is immersed. The comprehensive derivation of above governing equation can be found in [10, 4, 16] and [13]. In our analysis we assume that the initial height is equal to zero and to complete the description of the governing equation we have to choose a proper initial condition for velocity of the flow. In order to ensure the consistency of the problem we have to take $h^{\prime}(0)=0$. The discussion on the different choice of initial condition as well as the extensive analysis for more physical initial condition can be find in [20, 9, 17]. Furthermore, the rigorous mathematical analysis of the governing equation (1) for both cases of initial conditions Reader can find in [13]. Recently, a very interesting geometrical analysis of (1) appeared in [19] and the Reader is invited to compare these two different approaches.

As was shown in various experiments, the flow in a narrow tube may not be monotone [9, 15]. Thus a two regimes are distinguished where the fluid column height either monotonically increase or oscillates near the Jurin’s height (the point, where the capillary pressure balances the weight of a fluid). The change of the behavior of flow from monotone to oscillatory for certain values of nondimensional parameter was also proved analytically (see [13]). Nonlinear oscillations are very often difficult to analyze due to its complexity and unique character, hence, a non standard methods have to be used to analyze such a phenomenon efficiently. Nonlinear oscillations modelled by differential equations has been extensively analyzed in many papers (see for ex. [7, 12, 8]). A nonlinear dynamical systems are also widely studied in two fundamental monographs [5] and [11]. That monographs concern mostly the theoretical and applicable mathematical results related to the dynamical systems and nonlinear oscillations.

To proceed further with the asymptotic analysis we have to transform (1) into a more convenient form. After reducing the governing equation into a dimensionless form (see [10, 13]) and using a simple transformation for both independent variable and height function (for more details see [13]) we are able to write (1) in following form

[TABLE]

with

[TABLE]

where, prime denotes a differentiation with respect to the independent variable. In subsequent analysis we will rather use parameter $\epsilon=1/\sqrt{\omega}$ instead of $\omega$. Whereas, the parameter $\omega$ is defined as follows

[TABLE]

Obviously, we are not able to solve (2) explicitly. Hence, a advanced approximation methods need to be used to analyzed the above equation and its solution in most accurate manner.

The following article is composed of the two main sections concerning the perturbation with respect to $\epsilon\rightarrow 0^{+}$ and asymptotic behavior for large times for (2). In the next section we use a perturbation analysis to investigate the solution of the main equation for small values of $\epsilon$. Due to appearance of the non-Lipschitz function in the main equation we can not use any standard theorem to prove the convergence. Therefore, a new approach was proposed to prove that the exact solution approaches to the zeroth order approximation as $\epsilon\to 0^{+}$. Moreover, an asymptotic behavior of the solution for small times is also presented.

In Section 3 we present an asymptotic analysis of the oscillatory behavior for (2) for large times. For this reason we consider only the values of dimensionless parameter $\epsilon$ for which the solution oscillates around the Jurin’s height ($\epsilon<2$). We introduce a convergent method of finding a successive approximations of the solution to (2). The numerical analysis revealed that for a certain values of $\epsilon$ the error between the exact solution and the approximate one might be negligible even for the relatively small values of independent variable.

2 Perturbation analysis for small values of $\epsilon$

Let us now recall the main equation with appropriate initial conditions

[TABLE]

Using regular perturbation methods with respect to the small parameter $\epsilon$ we can obtain the zeroth order approximation to the exact solution of (5)

[TABLE]

The above function is depicted on the Figure 1. We can see that it represents an oscillation with a period equal to $6$.

To see how to obtain $u_{0}$ first we put $\epsilon=0$ into (5) and get a second order equation without the derivative term

[TABLE]

Since the above does not explicitly contain the independent variable we can introduce another function $g(u_{0}):=u_{0}^{\prime}$ and using $u_{0}^{\prime\prime}=g\,dg/du_{0}$ we are able to reduce the order what leads to

[TABLE]

where $g$ is now the sought function while $u$ is treated as an independent variable. We can easily integrate the above equation to get $g(u_{0})$. Then, using the definition of $g$ we have

[TABLE]

If we integrate above equation we get two solutions but only one of them is oscillatory for $\epsilon=0$ with $u_{0}(s)\geq 0$ and $s\geq 0$. Thus, we have

[TABLE]

To verify the periodicity of (6) we multiply (5) by $u^{\prime}$ and integrate to get the energy equation

[TABLE]

On the left hand side we have a sum of the kinetic and potential energy. It is easy to notice that the potential energy has a minimum. Moreover, for the boundary points of subintervals given in (6) for all $i$, the minimum value of potential energy is not achieved. Hence, the solution will oscillate with decaying amplitude (for a different treatment see [19]).

The above considerations were only formal and we do not know whether the exact solution of (5) converges to $u_{0}(s)$ as $\epsilon\to 0^{+}$. The standard convergence theorems from perturbation theory cannot be applied here because the Lipschitz condition for nonlinearity is not satisfied. As we will see with a suitable regularization we are able to avoid this requirement.

Before we proceed to the main result we will derive certain inequalities that will be useful in subsequent analysis. The following Lemma provides an asymptotic behavior of function $u(s)$ and its derivative $u^{\prime}(s)$ for $s\to 0^{+}$.

Lemma 1**.**

If $u(s)$ is a solution of (5) then for $s\in[0,2/(1+\epsilon)]$ we have the following estimates.

1. a)

Function $u$:

[TABLE] 2. b)

The derivative $u^{\prime}$:

[TABLE]

Furthermore, for $s\to 0^{+}$ the following asymptotic behavior holds

1. i)

Function $u$:

[TABLE] 2. ii)

The derivative $u^{\prime}$:

[TABLE]

Proof.

From [13] we have

[TABLE]

We see that the function of $s$ on the leftmost side attains its maximal value at $s=2/(1+\epsilon)$. Moreover, $1-\frac{1}{3}(1+\epsilon)s$ is a decreasing function of $s$ hence we have

[TABLE]

Thus (12) follows from the above. From (11) we write

[TABLE]

and using (12) we arrive at

[TABLE]

Specifically, we have

[TABLE]

what completes the first part of proof.

Now, it immediate to see that (14) is a consequence of (12). Next, let us consider the derivative $u^{\prime}(s)$ for small argument $s$. From (13) we have

[TABLE]

for $s\in[0,2/(1+\epsilon)]$. We see that the absolute value of derivative of $u$ is bounded by $s$ with $u^{\prime}(0)=0$. Now, looking on (5) we notice that for a very small argument the second derivative is positive which implies that the function $u^{\prime}$ must be a increasing function on the interval $[0,\delta]$ for certain $\delta>0$. Let us rewrite (5) in the energy form for $s=\delta$

[TABLE]

Taking advantage of the fact that function $u^{\prime}(t)$ is increasing on the interval $[0,\delta]$ we can write

[TABLE]

Hence, combining the above with (22) along with using (12) we obtain

[TABLE]

and finally

[TABLE]

where $M_{1}$ is a positive constant, $\delta\ll 1$ and $\delta\leq 1-M_{1}^{2}$. From (21) and the above inequality it follows that

[TABLE]

for $\delta\to 0^{+}$ what completed the proof. ∎

Now, we are able to introduce the main result in which we show that the solution of (5) converges uniformly to (6) as $\epsilon\to 0^{+}$.

Theorem 1**.**

Let $u(t)$ be the solution of (5). Then, the following asymptotic behaviour holds

[TABLE]

uniformly for $s\in[0,T]$ for any fixed $T>0$.

Proof.

Let us consider the equation (5) on the interval $[0,T]$, where $T$ is any fixed value greater than zero. Write

[TABLE]

where function $u_{0}$ is defined as in (6) and function $R(s)$ is a unknown remainder of a zeroth order approximation to the exact solution of (5). We will find an upper bound on $R(s)$. We substitute (28) into (5) and obtain following equation

[TABLE]

We can easily transform above to the form of integral equation

[TABLE]

where

[TABLE]

Thereafter we divide the interval $[0,T]$ into $[0,\delta]$ and $[\delta,T]$. We assume that the number $\delta$ is small. Its exact value will be found later.

Firstly, let us consider the interval $[0,\delta]$. We can estimate the integral in (30) from the above as follows

[TABLE]

It is a matter of simple calculus exercise to see that $u_{0}$ is bounded by $1$. Moreover, from [13] we know that the exact solution of (5) is no greater than $9/8$, hence the function $f$ can be easily bounded from above by $3/2+\epsilon$. Whereas, the function $G(z)$ is concave which together with $G(s-x)\leq G(\delta)$ gives $G(s-x)\leq\delta$. Therefore,

[TABLE]

Next, let us consider subinterval $[\delta,T]$. To this end, define a new function $w_{\delta}$ which satisfies the following initial value problem

[TABLE]

where $u_{\delta}=u(\delta),\ u^{\prime}_{\delta}=u^{\prime}(\delta)$ and $u$ is a solution of (5). Then we can estimate an error on the interval $[\delta,T]$ as follows

[TABLE]

Now, we will analyze the error function $R_{1}$ and $R_{2}$ separately. First, let us insert $u=R_{1}+w_{\delta}$ into eq. (5) to get

[TABLE]

where $R_{1}$ is an unknown function. Next, we are using (33) to introduce above equation in a more familiar form

[TABLE]

with the initial conditions $R_{1}(\delta)=u(\delta)-u_{\delta}=0$ and $R_{1}^{\prime}(\delta)=u^{\prime}(\delta)-u^{\prime}_{(}\delta)=0$. Subsequently, we integrate (36) twice and use a initial conditions of function $R_{1}$ to obtain

[TABLE]

where function $G(s-x)$ and $f(w_{\delta},R_{1})$ are defined like in eq. (30). If we already have the integral equation, we can proceed to the further analysis of the eq. (37) to get a appropriate integral inequality

[TABLE]

Notice that for $s\in[\delta,T]$ we have $u(s)>0$. Using (12) we have

[TABLE]

and, therefore, we can estimate the denominator in (38). Now, we define a new function

[TABLE]

Then, of course

[TABLE]

Differentiating (40) with respect to $s$ leads to

[TABLE]

Multiplying by $\exp{-(3Ts/\delta)}$ and integrating implies that

[TABLE]

where we have used the fact that $S(\delta)=0$. Now, if we choose $\delta(\epsilon):=-3T^{2}/\log{\epsilon}$ we see that $\delta(\epsilon)\to 0$ and $\max_{s\in[\delta,T]}|S(s)|\to 0$ as $\epsilon\to 0$. Hence, the remainder $R_{1}$ vanishes.

Now, we take a look at the second component on the right hand side in (34), namely the remainder denoted by $R_{2}$. Let us notice that the functions $w_{\delta}$ and $u_{0}$ satisfy the same nonlinear differential equation but with different initial conditions at point $s=\delta$. Using the results concerning function $u_{0}$ presented at beginning of the section we obtain that $w_{\delta}$ is also a periodic function. Multiplying (33) by $w_{\delta}^{\prime}$ and integrating we get

[TABLE]

for $s\in[\delta,T]$. Using above equation we can show that the function $w_{\delta}$ is bounded. Specifically, following inequalities are satisfied

[TABLE]

for all $s\in[\delta,T]$ and all $\delta\geq 0$. Next, performing some technical and tedious calculations it can be shown that $\lim_{\delta\to 0^{+}}w_{\delta,min}=0$ and $\lim_{\delta\to 0^{+}}w_{\delta,max}=9/8$. We expect that the function $w_{\delta}$ converges to $u_{0}$ as $\delta\to 0^{+}$. In our previous work [13] we proved that (5) posses a unique solution with zero initial conditions and with initial conditions of the form: $u(0)=\beta^{2}/2$ and $u^{\prime}(0)=0$, where $\beta$ is some arbitrary number. The uniqueness of the solution can be easily extended to all initial conditions belonging to the set $\Omega:=\{(x_{0},y_{0}):y_{0}^{2}/2+2\sqrt{2}/3x_{0}^{3/2}-x_{0}\leq 0\}$. Furthermore, the uniqueness theorem is true regardless of the value of the parameter $\epsilon$. Therefore, setting $\epsilon=0$ in (5) implies that the solution of

[TABLE]

exists and is unique for all initial conditions from the set $\Omega$. Notice that the initial conditions of function $w_{\delta}$ belongs to the set $\Omega$. Setting $x:=u$ and $y:=u^{\prime}$ we can transform eq. (46) into a system of differential equations

[TABLE]

The uniqueness of the solution of (46) for all initial condition belonging to the set $\Omega$ implies that all trajectories of the system (47) located in the set $\Omega$ are uniquely determined. Let $\Gamma_{w_{\delta}}:=\{(w_{\delta}(s),w^{\prime}_{\delta}(s)),\ \delta\leq s\leq T\}$ and $\Gamma_{u_{0}}:=\{(u_{0}(t),u^{\prime}_{0}(t)),\ \delta\leq s\leq T\}$ denote the trajectories associated with $w_{\delta}$ and $u_{0}$ respectively. After a detailed analysis of $w_{\delta}$ and $u_{0}$, we get $\Gamma_{u_{0}}\in\partial\Omega$ and $\Gamma_{w_{\delta}}\in\Omega$ for all $\delta\geq 0$. If a value of parameter $\delta$ will decrease to zero we anticipate that the elements of $\Gamma_{w_{\delta}}$ will be approaching the proper elements in $\Gamma_{u_{0}}$. Using the fact that the minimal and maximal value of $w_{\delta}$ converge to [math] and $9/8$ respectively we get

[TABLE]

Now, applying a non-intersection property for uniquely defined phase plane of system (47) we get $\lim_{\delta\to 0^{+}}\Gamma_{w_{\delta}}=\Gamma_{u_{0}}$. Therefore, combining that with the fact that the difference between function $w_{\delta}$ and $u_{0}$ at initial point converge to zero as $\delta\to 0^{+}$ we get $\lim_{\delta\to 0^{+}}\max_{s\in[\delta,T]}||(w_{\delta}(s),w^{\prime}_{\delta}(s))-(u_{0}(s),u^{\prime}_{0}(s))||=0$ what immediately implies $\lim_{\delta\to 0^{+}}|R_{2}|=\lim_{\delta\to 0^{+}}|w_{\delta}-u_{0}|=0$. Here the $||\cdot||$ is a Euclidean norm in $\mathbb{R}^{2}$

Finally, we take togather the inegualities (32) and (34) concerning the bound for remainder on the intervals $[0,\,\delta]$ and $[\delta,\,T]$ respectively and using $\lim_{\delta\to 0^{+}}|w_{\delta}-u_{0}|=0$ we conclude the final convergence result (27).

∎

The exact order of convergence in Theorem 1 is difficult to obtain because of the appearance of the non-Lipschitz function in the main equation. However, the numerical analysis shows that the error rate between $w_{\delta}$ and $u_{0}$ should not be worse than $\delta$ as $\delta\to 0^{+}$. In the following remark we used the MATHEMATICA software to find a numerical solution to eq. (7) with the non zero initial conditions.

Remark 1**.**

The numerical analysis revealed that the following convergence error holds

[TABLE]

where $u(s)$ is an exact solution of (5), function $u_{0}$ is defined as before and $T$ is an any fixed constant greater than zero.

From the proof of Theorem 1 we have that on the interval $[0,\delta]$ the error function between $u(s)$ and $u_{0}$ can be bounded from above by

[TABLE]

Whereas on the interval $[\delta,T]$ we get

[TABLE]

We numerically compute the values of the second term in above inequality for various values of $\delta$. The result are presented on Fig. 2. In calculation we use the asymptotic results from Lemma 1, specifically we use following initial conditions for function $w_{\delta}(s)$ on the interval $[\delta,T]$

[TABLE]

where $D_{1}$ and $D_{2}$ are some positive constants. The numerical results showed that the error between function $w_{\delta}$ and $u_{0}$ might of order of $\delta$. It is worth to emphasize that for some choice of initial conditions $w_{\delta}$ the quadratic order of convergence is observed.

Combing numerical simulations with the above results we conclude that (50) holds.

3 Asymptotic behavior for large times

In this section we consider a large time behaviour of the solution of (5) in the oscillatory regime. To facilitate the analysis it is useful to introduce an additional function $\phi$ that will transform (5) into a simpler form. We set

[TABLE]

where the exact formula for $\phi(s)$ will be found later. Let us notice that (54) is simply the Liouville transformation. It is easy to see that $u(s)=1/2+\phi(s)v(s)$ and after substituting the above into (5) we obtain

[TABLE]

We can see that in order to eliminate the first derivative we have to take

[TABLE]

Now we can rewrite (5) as

[TABLE]

We can think about the above equation as one describing nonlinear oscillations in which the frequency is modulated with the amplitude via the term in the parentheses.

We will use integral equations techniques to investigate the large time regime. Let us consider a time interval $[T,T^{*}]$, where $T$ and $T^{*}$ are fixed numbers. With a standard use of the Green’s function we can transform (57) into the integral form

[TABLE]

where

[TABLE]

and

[TABLE]

The initial values $v_{T}:=v(T),\ v^{\prime}_{T}:=v^{\prime}(T)$ are well defined because of our preceding results stating that (5) has a unique global solution (see [13]).

Immediately, we can make some initial observations about the function $v$.

Corollary 1**.**

Let $v$ be a solution of (57) and assume that $\epsilon>0$. Then, for sufficiently large $s>0$ we have

[TABLE]

where $C_{1,2}$ are some fixed constants.

Proof.

From our previous result [13] we know that $0\leq u\leq 9/8$ where $u$ is a solution of (5). We will show that the there exist such $T$ that the inequality $0<1/2+C_{1}\leq u$ is satisfied on the interval $[T,\infty)$. To do this let rewrite the energy form of the governing equation (5)

[TABLE]

Let assume that $u$ reaches zero at some point $s^{*}>0$. Then, the above immediately leads to a contradiction since the right-hand side is strictly negative. Hence the only point in which function $u$ can reach zero is $s=0$.

Moreover, from our previous work [13] we know that the solution of (5) oscillates around $1/2$ and approaches it as $s\rightarrow\infty$. Hence, we are able to choose such $T_{1}$ that the inequality $0<1/2+C_{1}\leq u$ is satisfied for $s\in[T_{1},\infty)$. Using the transformation (54) yields the estimate from below. Similarly, on the grounds of the asymptotic convergence of function $\phi v$ there must exist $T_{2}$ such that the inequality $\phi v\leq C_{2}<1/2$ is satisfied for $s\in[T_{2},\infty)$. Finally, we put $T:=\max\{T_{1},T_{2}\}$. ∎

Now, we are able to introduce result concerning existence and uniqueness of solutions to (58). Later, we will use it to construct several large time approximations.

Theorem 2**.**

Solution of (58) exists and is uniquely defined on $[T,\infty)$, where $T$ is a fixed number defined as in the proof of Corollary (1) .

Proof.

Let us consider the existence problem of the solution of (58) on the interval $[T,T^{*}]$, where $T^{*}$ is an arbitrary fixed number greater than $T$. Moreover we are looking for the solution of (58) in the Banach space $(C([T,T^{*}]),\left\|\cdot\right\|_{\infty}$). Further, introduce the Bielecki’s norm defined as follows [1]

[TABLE]

where $\alpha$ is a certain fixed number. Because of the the fact that Bielecki’s norm is equivalent to the usual supremum norm on $[T,T^{*}]$ we can use it in the proof that the operator $A$ is a contraction with respect to the supremum norm.

For any continuous $u$ and $v$ we have

[TABLE]

where the last inequality is a consequence of Corollary 1 with appropriate chosen value of $T$. Finally, using (58) we can make the estimate

[TABLE]

Multiplying the above by $\exp(-\alpha s)$, taking supremum and using (64) leads to

[TABLE]

Choosing a sufficiently large $\alpha$ we can make the prefactor on the right-hand side of the above to be smaller than $1$ which shows that the operator $L$ is a contraction in the Bielecki’s norm. The equivalence of the supremum and Bielecki’s norms implies that $L$ is also a contraction in the supremum norm. Banach Contraction Mapping Theorem used on the space $(C([T,T^{*}]),\left\|\cdot\right\|_{\infty})$ asserts that (58) posses a unique fixed point. From the arbitrariness of $T^{*}$ we conclude that the solution can be continued to the unbounded domain $[T,\infty)$. ∎

The above result suggests that we may be able to explicitly construct the exact solution of the problem (58) as a limit of iterative sequence. This can be utilized in finding the leading order asymptotic behaviour for large times and provide some useful approximations.

To this end, let us notice that the function $f$ in (57) can be expanded in a Taylor series with respect to the product $\phi(s)v(s)$. The function $\phi(s)$ decreases monotonically to zero while the function $v$ is bounded. Therefore, for $s\to\infty$ we have $\phi(s)v(s)\to 0$. Hence, we can rewrite (58) in the following form

[TABLE]

where $A_{0}$ and $B_{0}$ are defined as before and the coefficients are

[TABLE]

Equation (67) can now be used to construct a iteration

[TABLE]

Theorem 3**.**

Let $v^{*}$ be the unique solution of (67). Then $\{v_{n}\}$ is convergent to unique $v^{*}$. Moreover,

[TABLE]

with

[TABLE]

where

[TABLE]

and

[TABLE]

and $\nu=\tau(j-2l+k+1-2p)$.

Proof.

Let us prove that the iteration

[TABLE]

converges to the exact solution of (58) on the interval $[T,T^{*}]$, where $T^{*}$ is an arbitrary fixed number greater than $T$. Take any two solutions of (58), call them $u$ and $v$ and calculate

[TABLE]

Then, we can apply Lagrange’s Mean Value Theorem for the function $(\phi v)^{m+1}$ for all values of $t$. Thus, we have

[TABLE]

where $w$ is a some function satisfying $\min\{\phi(t)v(t),\phi(t)u(t)\}<w(t)<\max\{\phi(t)v(t),\phi(t)u(t)\}$ for all $t\in[T,s]$. From the Cor. 1 we have that $\phi(s)v(s)\leq C_{2}<1/2$, where $v$ is a solution of (58). Hence we can finally write

[TABLE]

where $||\cdot||_{B}$ is a Bielecki’s norm defined as in (63) and $C_{2}$ is a fixed number smaller than $1/2$ introduced earlier in (1).

From the standard convergence criterion we know that the series in above converge. Moreover, we have

[TABLE]

where ${}_{2}F_{1}$ is a hypergeometric function

[TABLE]

with

[TABLE]

Now we are able to rewrite the main inequality of (75) in the following form

[TABLE]

Let us notice the inverse dependence of function $H(\alpha,T)$ on the parameter $\alpha$. Hence, choosing sufficiently large $\alpha$ we can make expression in brackets smaller than $1$ and whence show that $P$ is a contraction in Bielecki’s metric on considered interval.

From the arbitrariness of $T^{*}$ we conclude that the operator $P$ is a contraction on the whole interval $[0,\infty]$.

From the numerical analysis and observations for large arguments it is justified to approximate function $v$ by some trigonometric function. Therefore, the natural candidate for the initial guess is $v_{0}$ satisfying $Lv_{0}=0$, with appropriate initial conditions. Hence,

[TABLE]

Next, we can expand the kernel $G(z)$ in the integral in (67) and get

[TABLE]

It is easy notice that if a $v_{0}$ is a sum of trigonometric functions then the subsequent result of iteration will be also in the similar form. In the another words, function $v_{n}$ will contain only a linear combination of functions $\sin{\tau s}$ and $\cos{\tau s}$. We will use the mathematical induction to prove that the function $v_{n}$ can be reduced to the form

[TABLE]

for all $n\in\mathbb{N}$. As will be shown the $A_{n}$ and $B_{n}$ coefficients in fact dependent only on the values of $A_{n-1}$, $B_{n-1}$ and on initial coefficients $A_{0}$ and $B_{0}$.

The base step of mathematical induction of this issue is crucial. From (82) we have that the function $v_{0}$ is a linear sum of sine and cosine functions. If we substitute the $v_{0}$ to (83) and use a binomial expansion we get

[TABLE]

where

[TABLE]

Now, we will compute the values of above integrals to get the analytic form of the function $v_{n}(s)$ for large values of $s$. Let us first consider integral $J_{1}$

[TABLE]

The imaginary parts of above expressions have to sum up to zero. Using the Euler’s formula we expanded the real trigonometric functions and we removed only the vanishing terms at the upper limit of integration which does not affect on the realness of the final function. Hence we get the final expression (72) of $J_{1}$ for $v_{1}$. The above steps can be repeated as well for $J_{2}$ to get final asymptotic form of that integral.

As we have seen the asymptotic expression for function $v_{1}$ is a linear sum of sine and cosine functions, where the coefficients are dependent only on $A_{0}$ and $B_{0}$. Moreover, if we now assume that the asymptotic behavior of function $v_{n}$ for some $n$ is expressed by the formula (84) then repeating the tedious calculations performed above for function $v_{1}$ we get

[TABLE]

Since both the base case and the inductive step have been performed, by mathematical induction we conclude the asymptotic form of function $v_{n}$ (84) for arbitrary $n\in\mathbb{N}$. ∎

The asymptotic solution (70) of eq. (58) describes the exact solution very well, especially for large value of T. As a illustration of the above result we present some numerical simulations. The asymptotic solution (70) with $n=0$ and $n=6$ is depicted in Figs. 4 and 4. For simplicity we used a finite sum in formulas (71) to calculate an adequate coefficients. Hence, in numerical analysis it is more convenient to use the following approximate expression for $A_{n}$ and $B_{n}$

[TABLE]

where $J_{1}$ and $J_{2}$ are defined as before and $N=30$ in our analysis.

4 Conclusion

We have provided the rigorous asymptotic analysis of a governing equation modelling the capillary rise phenomenon in a vertical narrow tube. The perturbation analysis of the main equation was a nontrival task due the a appearance of a non-Lipschitz component in the equation. However, the new approach was proposed where the convergence of the exact solution to the unperturbed solution was considered separately on two intervals dependent on the decreasing parameter. For one interval the nonlinear component of the governing equation is not a Lipschitz continuous function. To deal with this problem we use the fact that the length of that interval decreases with $\epsilon$ what allow us to estimate the remainder from above by a function decreasing to zero as $\epsilon$ goes to zero. We conclude the final statement performing a comprehensive calculations.

In the subsequent section we investigated the asymptotic behavior of a solution for large values of independent variable. The method of finding successive approximation to the exact oscillatory solution was introduced. Finally, the numerical analysis was presented to compare the numerical solution of the governing equation with the approximate periodic asymptotic solutions.