Asymptotic behaviour of a solution to a nonlinear equation modelling capillary rise
Łukasz Płociniczak Mateusz Świtała
math.CA

Experimental investigation of superdiffusion via coherent disordered Quantum Walks
    Paolo Mataloni^1
    December 30, 2023
==================================================================================




	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. 


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




§ 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 <cit.> and <cit.>), 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 <cit.>). 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 <cit.> we provided brief explanation of each component in the main governing equation modelling the capillary rise process in a narrow vertical tube 
 
    8 ν/r^2 h h'+ρ g h+ρ d/ d t(h h')=2 γcosθ/r,
 
	where h=h(t) is the liquid column height at time t, r is a radius of the capillary tube, and μ, ρ denote viscosity and density of liquid respectively. As usual, g denotes the gravitational acceleration and γ denotes the surface tension. It is worth emphasis that in the current work the contact angle, here denoted by θ, preserves a constant value during the flow. Nevertheless, some results <cit.> indicates the dependence of contact angle on the velocity of flow. The meticulous investigation of the (<ref>) 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 <cit.> and <cit.>. 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'(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 <cit.>. Furthermore, the rigorous mathematical analysis of the governing equation (<ref>) for both cases of initial conditions Reader can find in <cit.>. Recently, a very interesting geometrical analysis of (<ref>) appeared in <cit.> 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 <cit.>. 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 <cit.>). 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. <cit.>). A nonlinear dynamical systems are also widely studied in two fundamental monographs <cit.> and <cit.>. 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 (<ref>) into a more convenient form. After reducing the governing equation into a dimensionless form (see <cit.>) and using a simple transformation for both independent variable and height function (for more details see <cit.>) we are able to write (<ref>) in following form
	 
    u”(s)+1/√(ω)u'(s)+√(2 u(s))=1,
 
		with 
		
    u(0)=0,    u'(0)=0,

		where, prime denotes a differentiation with respect to the independent variable. In subsequent analysis we will rather use  parameter ϵ=1/√(ω) instead of ω. Whereas, the parameter ω is defined as follows
		
    ω = ρ ^3 g^2 r^5/128 μ ^2 γcosθ.

	    Obviously, we are not able to solve (<ref>) 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 ϵ→ 0^+ and asymptotic behavior for large times for (<ref>). In the next section we use a perturbation analysis to investigate the solution of the main equation for small values of ϵ. 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 ϵ→ 0^+.  Moreover, an asymptotic behavior of the solution for small times is also presented.
	
	In Section <ref> we present an asymptotic analysis of the oscillatory behavior for (<ref>) for large times. For this reason we consider only the values of dimensionless parameter ϵ for which the solution oscillates around the Jurin's height (ϵ < 2). We introduce a convergent method of finding a successive approximations of the solution to (<ref>). The numerical analysis revealed that for a certain values of ϵ the error between the exact solution and the approximate one might be negligible even for the relatively small values of independent variable.
	


§ PERTURBATION ANALYSIS FOR SMALL VALUES OF Ε

Let us now recall the main equation with appropriate initial conditions

    u”(s)+ϵ u'(s)+√(2 u(s))=1,   u(0)=0,  u'(0)=0.

Using regular perturbation methods with respect to the small parameter ϵ we can obtain the zeroth order approximation to the exact solution of (<ref>)

    u_0(s)=1/72∑_i=1^∞(s-6(i-1))^2 (6 i-s)^21 _{s∈ [6(i-1),6 i]}.

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



To see how to obtain u_0 first we put ϵ=0 into (<ref>) and get a second order equation without the derivative term

    u_0”+√(2 u_0)=1,    u_0(0)=0,   u_0'(0)=0.

Since the above does not explicitly contain the independent variable we can introduce another function g(u_0):=u_0' and using u_0”=g   d g/d u_0 we are able to reduce the order what leads to

    g d g/d u_0+√(2 u_0)=1,    g(0)=0,
 
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

    u_0'(x)=±√(2/3)√(-(2 √(2)√(u(x))-3) u_0(x)),    u_0(0)=0.

If we integrate above equation we get two solutions but only one of them is oscillatory for ϵ=0 with u_0(s)≥ 0 and s≥ 0. Thus, we have

    u_0(s)=s^2 (6-s)^2/721 _{s∈ [0,6 ]}.


To verify the periodicity of (<ref>) we multiply (<ref>) by u' and integrate to get the energy equation

    u'^2/2+2 √(2)/3u^3/2-u=-ϵ∫_0^su'(x)^2 d  x

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  (<ref>) 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 <cit.>).

The above considerations were only formal and we do not know whether the exact solution of (<ref>) converges to u_0(s) as ϵ→ 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'(s) for s→ 0^+.

	If u(s) is a solution of (<ref>) then for  s∈ [0,2/(1+ϵ)] we have the following estimates.
	
		
  *  Function u:
	
    1/3s^2/2≤ u(s)≤s^2/2.

	
  *  The derivative u':
	
    u'^2(s)≤ s^2-2 √(2)/9s^3.

	
	Furthermore, for s→ 0^+ the following asymptotic behavior holds

    
  i) Function u:
    
    u(s)=𝒪(s^2),  as   s→ 0^+.

    
  ii) The derivative u':
    
    u'(s)=𝒪(s),  as   s → 0^+.




	From <cit.> we have 
	
    (1-1/3(1+ϵ)s)s^2/2≤ u(s)≤s^2/2  for   s∈[0,3/1+ϵ].

	We see that the function of s on the leftmost side attains its maximal value at s=2/(1+ϵ). Moreover, 1-1/3(1+ϵ)s is a decreasing function of s hence we have 
	
    1/3=(1-1/3(1+ϵ)(2/1+ϵ))≤(1-1/3(1+ϵ)s)  for   s∈[0,2/(1+ϵ)].

	Thus (<ref>) follows from the above. From (<ref>) we write
	
    u'^2/2+2 √(2)/3u^3/2-u≤ 0.

	and using (<ref>) we arrive at
	
    u'^2(s)≤ 2 u-4 √(2)/3u^3/2≤ s^2-2 √(2)/9s^3.

	Specifically, we have
	
    |u'(s)|≤ s,

	what completes the first part of proof.
	
	Now, it immediate to see that (<ref>) is a consequence of (<ref>). Next, let us consider the derivative u'(s) for small argument s. From (<ref>) we have
    
    |u'(s)|≤√(s^2-2√(2)/9s^3)≤ s,

    for s∈ [0,2/(1+ϵ)]. We see that the absolute value of derivative of u is bounded by s with u'(0)=0. Now, looking on (<ref>) we notice that for a very small argument the second derivative is positive which implies that the function u' must be a  increasing function on the interval [0,δ] for certain δ>0. Let us rewrite (<ref>) in the energy form for s=δ
    
    1/2u'(δ)^2+2√(2)/3u(δ)^3/2-u(δ)=-ϵ∫_0^δu'(t)^2  d t.

    Taking advantage of the fact that function u'(t) is increasing on the interval [0,δ] we can write
    
    -ϵ∫_0^δu'(t)^2  d t ≥ -ϵδ u'(δ)^2.

    Hence, combining the above with (<ref>) along with using (<ref>) we obtain
    
    u'(δ)^2≥1/3(1/2+ϵδ)(δ^2-δ^3),

    and finally 
    
    u'(δ)≥ M_1 1/√(3(1/2+ϵδ))δ,

    where M_1 is a positive constant, δ≪ 1 and δ≤ 1- M_1^2. From (<ref>) and the above inequality it follows that
    
    M_1 1/√(3(1/2+ϵδ))δ≤ u'(δ)≤δ,

    for δ→ 0^+ what completed the proof.

Now, we are able to introduce the main result in which we show that the solution of (<ref>) converges uniformly to (<ref>) as ϵ→ 0^+. 

	Let u(t) be the solution of (<ref>). Then, the following asymptotic behaviour holds
	
    lim_ϵ→ 0^+|u(s)-u_0(s)|=0,

	uniformly for s∈ [0,T] for any fixed T>0.  


	Let us consider the equation (<ref>) on the interval [0, T], where T is any fixed value greater than zero. Write 
	
    u(s)=u_0(s)+R(s),

	where function u_0 is defined as in (<ref>) and function R(s) is a unknown remainder of a zeroth order approximation to the exact solution of (<ref>). We will find an upper bound on R(s). We substitute (<ref>) into (<ref>) and obtain following equation
	
    R”(s)=√(2 u_0(s))-√(2(u_0(s)+R(s))).

	We can easily transform above to the form of integral equation
	
    R(s)=∫_0^sG(s-x)f(u_0(x),R(x)) d x
 
	where 
	
	G(z)    = 1/ϵ(1-e^-ϵz),

	f(u_0(s),R(s))    =√(2 u_0(s))-√(2(u_0(s)+R(s)))-ϵu_0'(s).
	 
	Thereafter we divide the interval [0,T] into [0,δ] and [δ, T]. We assume that the number δ is small. Its exact value will be found later.
	
	Firstly, let us consider the interval [0,δ]. We can estimate the
	integral in (<ref>) from the above as follows
	
    |R(s)|≤∫_0^δ|G(s-x)||f(u_0(x),R(x))| d x≤δ∫_0^δ|f(u_0(x),R(x))| d x.

	It is a matter of simple calculus exercise to see that u_0 is bounded by 1. Moreover, from <cit.> we know that the exact solution of (<ref>) is no greater than 9/8, hence the function f can be easily bounded from above by 3/2+ϵ. Whereas, the function G(z) is concave which together with G(s-x)≤ G(δ) gives G(s-x)≤δ. Therefore, 
 
    |R(s)|≤(3/2+ϵ)δ^2,    s∈[0,δ].

 Next, let us consider subinterval [δ, T]. To this end, define a new function w_δ which satisfies the following initial value problem
 
    w_δ”+√(2 w_δ)=1,   w_δ(δ)=u_δ, w_δ'(δ)=u'_δ,
 
 where u_δ = u(δ), u'_δ=u'(δ) and u is a solution of (<ref>). Then we can estimate an error on the interval [δ,T] as follows
 
    |R(s)|=|u(s)-u_0(s)|≤ |u(s)-w_δ(s)|+|u_0(s)-w_δ(s)|=|R_1|+|R_2|.

Now, we will analyze the error function R_1 and R_2 separately. First, let us insert u=R_1+w_δ into eq. (<ref>) to get

    R_1”+w_δ”+ϵ (R_1+w_δ)+√(2(R_1+w_δ))=1,

where R_1 is an unknown function. Next, we are using (<ref>) to introduce above equation in a more familiar form 

    R_1”+ϵ R_1'=1-√(2(R_1+w_δ))-√(2 w_δ)-ϵ w'_δ,

with the initial conditions R_1(δ)=u(δ)-u_δ=0 and R_1'(δ)=u'(δ)-u'_(δ)=0.
Subsequently, we integrate (<ref>) twice and use a initial conditions of function R_1 to obtain

    |R_1|=|∫_δ^s G(s-x)f(w_δ(x), R_1(x))  d x|,

where function G(s-x) and f(w_δ, R_1) are defined like in eq. (<ref>).  If we already have the integral equation, we can proceed to the further analysis of the eq. (<ref>) to get a appropriate integral inequality
 
    |R_1(s)|   ≤∫_δ^s|G(s-x)||f(w_δ(x),R_1(x))| d x≤ T ∫_δ^s|f(w_δ(x),R_1(x))| d x
       ≤ T ∫_δ^s(|√(2 w_δ)-√(2(w_δ+R_1))|+ϵ|w'_δ|)  d x≤ T ∫_δ^s|R_1|/√(2 w_δ)+√(2(w_δ+R_1))+ϵ  dx
       ≤  T ∫_δ^s(|R_1|/√(2 u(δ))+ϵ)   dx

 
 
 Notice that for s∈[δ, T] we have u(s)>0. Using (<ref>) we have
 
    1/3δ≤√(2 u(δ)),

 and, therefore, we can estimate the denominator in (<ref>). Now, we define a new function
 
    S(s):=T ∫_δ^s(3 |R_1|/δ+ϵ)   dx.

 Then, of course
 
    |R_1(s)|≤ S(s).

 Differentiating (<ref>) with respect to s leads to
 
    d/d sS(s)=T(3 |R_1|/δ+ϵ) ≤ T(3 S/δ+ϵ)

 Multiplying by exp-(3 T s/δ ) and integrating implies that
 	
    S(s)≤ -ϵδ/3 T(e^-3 T/δs-e^-3 T)e^3 T/δs≤ϵδ/3 T(e^-3 T+3 T/δs-1)≤ϵδ/3 T e^-3 T+3 T^2/δ,

 where we have used the fact that S(δ)=0. Now, if we choose δ(ϵ):=-3 T^2/logϵ we see that δ(ϵ)→ 0 and max_s∈[δ,T]|S(s)|→ 0 as ϵ→ 0. Hence, the remainder R_1 vanishes. 
 	
 Now, we take a look at the second component on the right hand side in (<ref>), namely the remainder denoted by R_2. Let us notice that the functions w_δ and u_0 satisfy the same nonlinear differential equation but with different initial conditions at point s=δ. Using the results concerning function u_0 presented at beginning of the section we obtain that w_δ is also a periodic function. Multiplying (<ref>) by w_δ' and integrating we get
 	
    w_δ'(s) ^2/2+2√(2)/3w_δ(s)^3/2-w_δ(s)-u'_δ ^2/2-2√(2)/3u_δ ^3/2+u_δ=0,

 	for s∈ [δ, T]. Using above equation we can show that the function w_δ is bounded. Specifically, following inequalities are satisfied
 	
    min_s∈ [δ, T]{u_0(s)}=0≤ w_δ,min≤ w_δ(s)≤ w_δ, max≤9/8=max_s∈ [δ, T]{u_0(s)},

 	for all s∈ [δ, T] and all δ≥ 0. Next, performing some technical and tedious calculations it can be shown that lim _δ→ 0^+ w_δ, min=0 and lim _δ→ 0^+ w_δ, max=9/8. We expect that the function w_δ converges to u_0 as δ→ 0^+. In our previous work <cit.> we proved that (<ref>) posses a unique solution with zero initial conditions and with initial conditions of the form: u(0)=β^2/2 and u'(0)=0, where β is some arbitrary number. The uniqueness of the solution can be easily extended to all initial conditions belonging to the set Ω:={(x_0, y_0): y_0^2/2+2√(2)/3 x_0^3/2-x_0≤ 0}. Furthermore, the uniqueness theorem is true regardless of the value of the parameter ϵ. Therefore, setting ϵ=0 in (<ref>) implies that the solution of
 	
    u”+√(2 u)=1,

 	exists and is unique for all initial conditions from the set Ω. Notice that the initial conditions of function w_δ belongs to the set Ω. Setting x:=u and y:=u' we can transform eq. (<ref>) into a system of differential equations
 	
    {[      x'  =        y ,;       y' = 1-√(2 x) . ].

 	The uniqueness of the solution of (<ref>) for all initial condition belonging to the set Ω implies that all trajectories of the system (<ref>) located in the set Ω are uniquely determined. Let Γ_w_δ :={(w_δ(s),w'_δ(s)), δ≤ s≤ T} and Γ_u_0:={(u_0(t),u'_0(t)), δ≤ s≤ T} denote the trajectories associated with w_δ and u_0 respectively. After a detailed analysis of w_δ and u_0, we get  Γ_u_0∈∂Ω and Γ_w_δ∈Ω for all δ≥ 0. If a value of parameter δ will decrease to zero we anticipate that the elements of Γ_w_δ will be approaching the proper elements in Γ_u_0. Using the fact that the minimal and maximal value of w_δ converge to 0 and 9/8 respectively we get
 	
 	
 	
    lim_δ→ 0^+(w_δ,min,0)    = (0,0)=(u_0,min,0),
    lim_δ→ 0^+(w_δ,max,0)    = (9/8,0)=(u_0,max,0).

 	Now, applying a non-intersection property for uniquely defined phase plane of system (<ref>) we get lim_δ→ 0^+Γ_w_δ =Γ_u_0. Therefore, combining that with the fact that the difference between function w_δ and u_0 at initial point converge to zero as δ→ 0^+ we get lim_δ→ 0^+max_s∈ [δ, T]||(w_δ(s),w'_δ(s))-(u_0(s),u'_0(s))||=0 what immediately  implies lim_δ→ 0^+|R_2|=lim_δ→ 0^+|w_δ -u_0|=0. Here the ||·|| is a Euclidean norm in ℝ^2
 	
 	
 	Finally, we take togather the inegualities (<ref>) and (<ref>) concerning the bound for remainder on the intervals [0, δ] and [δ,  T] respectively and using lim_δ→ 0^+|w_δ -u_0|=0  we conclude the final convergence result (<ref>).
 	


The exact order of convergence in Theorem <ref> 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_δ and u_0 should not be worse than δ as δ→ 0^+. In the following remark we used the MATHEMATICA software to find a numerical solution to eq. (<ref>) with the non zero initial conditions. 

The numerical analysis revealed that the following convergence error holds

    max_s∈[0,T]|u(s)-u_0(s)|=𝒪(-1/logϵ)  as  ϵ→ 0^+,

where u(s) is an exact solution of (<ref>), function u_0 is defined as before and T is an any fixed constant greater than zero.

From the proof of Theorem <ref> we have that on the interval [0,δ] the error function between u(s) and u_0 can be bounded from above by

    |R(s)|≤(3/2+ϵ)δ^2.

Whereas on the interval [δ, T] we get

    |R(s)|≤δ/3 T e^-3 T+max_s∈ [δ, T]|w_δ(s)-u_0(s)|.

We numerically compute the values of the second term in above inequality  for various values of δ. The result are presented on Fig. <ref>. In calculation we use the asymptotic results from Lemma <ref>, specifically we use following initial conditions for function w_δ(s) on the interval [δ, T]

    w_δ(δ)=D_1 δ^2,   w_δ(δ)=D_2δ,

where D_1 and D_2 are some positive constants.
The numerical results showed that the error between function w_δ and u_0 might of order of δ. It is worth to emphasize that for some choice of initial conditions w_δ the quadratic order of convergence is observed. 

Combing numerical simulations with the above results we conclude that (<ref>) holds.



§ ASYMPTOTIC BEHAVIOR FOR LARGE TIMES


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

    v(s):=(u(s)-1/2)ϕ(s)^-1,

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

    ϕ v”+(2 ϕ' +ϵϕ) v'+(ϕ”+ϵϕ')v+√(1+2 ϕ v)=1.

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

    ϕ(s):=e^ϵ s/2

Now we can rewrite (<ref>) as

    v”+(2/1+√(1+2 ϕ v)-ϵ^2/4)v=0,   v(0)=-1/2, v'(0)=ϵ/4.

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 (<ref>) into the integral form

    v(s)=A_0cosτ s+ B_0sinτ s+∫_T^sG(s-t) f(ϕ(t),v(t))  d t=:L v,

where

    G(z)    = 1/τsinτz,    τ=√(1-ϵ^2/4),

    f(ϕ(t),v(t))     = (1-2/1+√(1+2 ϕ(t) v(t)))v(t),
 
and 

    A_0 =     v'(T
        ) cos (T τ )/τ+v(T) sin (T τ ),
    
        B_0 =      v(T) cos (T τ )-v'(T) sin (T τ )/τ,

The initial values v_T:=v(T), v'_T:=v'(T) are well defined because of our preceding results stating that (<ref>) has a unique global solution (see <cit.>). 

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

	Let v be a solution of (<ref>) and assume that ϵ > 0. Then, for sufficiently large s>0 we have
	
    -1/2<C_1≤ϕ(s)v(s)≤ C_2< 1/2,

	where C_1,2 are some fixed constants. 
  

From our previous result <cit.> we know that 0≤ u≤ 9/8 where u is a solution of (<ref>). We will show that the there exist such T that the inequality 0<1/2+C_1≤ u is satisfied on the interval [T,∞). To do this let rewrite the energy form of the governing equation (<ref>)
	
    u'^2/2+2√(2)/3u^3/2-u=-ϵ∫_0^su'^2 d t.

	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 <cit.> we know that the solution of (<ref>) oscillates around 1/2 and approaches it as s→∞. Hence, we are able to choose such T_1 that the inequality 0<1/2+C_1≤ u is satisfied for s∈[T_1,∞). Using the transformation (<ref>) yields the estimate from below. Similarly, on the grounds of the asymptotic convergence of function ϕ v there must exist T_2 such that the inequality ϕ v≤ C_2<1/2 is satisfied for s∈[T_2,∞). Finally, we put T:=max{T_1,T_2}.

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

Solution of (<ref>) exists and is uniquely defined on [T,∞), where T is a fixed number defined as in the proof of Corollary (<ref>) .




	Let us consider the existence problem of the solution of (<ref>) on the interval [T,T^*], where T^* is an arbitrary fixed number greater than T. Moreover we are looking for the solution of (<ref>) in the Banach space (C([T,T^*]),·_∞). Further, introduce the Bielecki's norm defined as follows <cit.>
	
    u_B:=sup_s∈[T,T^*]|e^-α su(s)|,  α>0,

	where α 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
	
    |f(ϕ(t), u(t)) -f(ϕ(t),v(t) )|≤   (1+2/√(1+2 ϕ u)+√(1+2 ϕ v))|u-v| 
    ≤   (1+1/√(1+2 C_1))|u-v|,

	where the last inequality is a consequence of Corollary <ref> with appropriate chosen value of T. Finally, using (<ref>) we can make the estimate
	
    |L(u)(s)-L(v)(s)|≤   ∫_T^s|G(s-t)| |f(ϕ(t),u(t))-f(ϕ(t),v(t))|  d t≤1/τ∫_T^s|f(ϕ(t),u(t))-f(ϕ(t),v(t))|  d t
    ≤   1/τ∫_T^se^α te^-α t|f(ϕ(t),u(t))-f(ϕ(t),v(t))|  d t≤e^α s-e^α T/ταf(ϕ,u)-f(ϕ,v)_B .
 
	Multiplying the above by exp(-α s), taking supremum and using (<ref>) leads to
	
    L(u)-L(v)_B≤1/τα(1+1/√(1+2 C_1))u-v_B.

	Choosing a sufficiently large α 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^*]),·_∞) asserts that (<ref>) posses a unique fixed point. From the arbitrariness of T^* we conclude that the solution can be continued to the unbounded domain [T,∞).


The above result suggests that we may be able to explicitly construct the exact solution of the problem (<ref>) 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 (<ref>) can be expanded in a Taylor series with respect to the product ϕ(s) v(s). The function ϕ(s) decreases monotonically to zero while the function v is bounded. Therefore, for s→∞ we have ϕ(s)v(s)→ 0. Hence, we can rewrite (<ref>) in the following form

    v(s)=A_0cosτ s+B_0sinτ s-∑_m=1^∞σ_m∫_T^sG(s-t) ϕ(t)^m v^m+1  d t=:P v,

where A_0 and B_0 are defined as before and the coefficients are

    σ_m=1/2m+12^m+1.

Equation (<ref>) can now be used to construct a iteration

    v_n+1=P v_n,     v_0=A_0 cosτ s+B_0 sinτ s.


Let v^* be the unique solution of (<ref>). Then {v_n} is convergent to unique v^*. Moreover, 

    v_n+1(s)∼ A_n sinτ s+B_n cosτ s,  as s→∞,

with

    A_n+1=   v'(T
        ) cos (T τ )/τ+v(T) sin (T τ )-1/τ∑_m=1^∞ ∑_j,k
     j+k=m+1m+1jσ_m A_n^jB_n^m+1-jJ_1,
    
            B_n+1=    v(T) cos (T τ )-v'(T) sin (T τ )/τ+1/τ∑_m=1^∞∑_j,k
     j+k=m+1m+1jσ_m A_n^jB_n^m+1-jJ_2,

where

    J_1=    -1/2^j+k+1∑_l=0^j∑_p=0^k+1jlk+1pf(j,l,k,p,m,τ,ϵ)/(τ(j-2 l+k+1-2 p))^2+λ^2,
    
            f(j,l,k,p,m,τ,ϵ)=   {[ (-1)^r+l+1(λcosν T-νsinν T)e^-λ T,                        for  j=2 r,;   (-1)^r+l(νsinν T-λcosν T)e^-λ T,                      for  j=2 r+1, ].

and

    J_2=    -1/2^j+k+1∑_l=0^j+1∑_p=0^kj+1lkpg(j,l,k,p,m,τ,ϵ)/(τ(j-2 l+k+1-2 p))^2+λ^2,
    
            g(j,l,k,p,m,τ,ϵ)=   {[   (-1)^r+l(λcosν T-νsinν T)e^-λ T,                      for  j=2 r+1,; (-1)^r+1+l(νcosν T-λsinν T)e^-λ T,                        for  j=2 r, ].

and ν=τ (j-2 l+k+1-2 p).


    Let us prove that the iteration
    
    v_n+1=P v_n,   v_0=A_0 sinτ s+B_0 cosτ s,

    converges to the exact solution of (<ref>) on the interval [T,T^*], where T^* is an arbitrary fixed number greater than T. Take any two solutions of (<ref>), call them u and v and calculate
    
    |P v-P u|≤∑_m=1^∞σ_m∫_T^s|G(s-t)|ϕ(t)^-1|(ϕ v)^m+1-(ϕ u)^m+1| d t.

        Then, we can apply Lagrange's Mean Value Theorem for the function (ϕ v)^m+1 for all values of t. Thus, we have
        
    ∑_m=1^∞σ_m∫_T^s|G(s-t)|ϕ(t)^-1|(ϕ v)^m+1-(ϕ u)^m+1| d t≤   ∑_m=1^∞σ_m(m+1)∫_T^s|G(s-t)|w(t)^m|v-u| d t,

where w is a some function satisfying min{ϕ(t) v(t),ϕ(t) u(t)} <w(t) < max{ϕ(t) v(t),ϕ(t) u(t)} for all t∈ [T,s]. From the Cor. <ref> we have that ϕ(s) v(s)≤ C_2<1/2, where v is a solution of (<ref>). Hence we can finally write

    |P v-P u|≤   ∑_m=1^∞σ_m(m+1)C_2^m∫_T^se^α t|G(s-t)|e^-α   t|v-u| d t
    ≤   v-u_B∑_m=1^∞σ_m(m+1)C_2^m∫_T^se^α t|G(s-t)| d t
    ≤   v-u_B∑_m=1^∞σ_m(m+1)C_2^m1/τ∫_T^se^α t d t
    
            =   e^α sv-u_B∑_m=1^∞σ_m(m+1)C^m1/τα(1-e^-α (s- T))
    ≤    e^α sv-u_B∑_m=1^∞σ_m(m+1)C_2^m1/τα,

    where ||·||_B is a Bielecki's norm defined as in (<ref>) and C_2 is a fixed number smaller than 1/2 introduced earlier in (<ref>).
    
    
     From the standard convergence criterion we know that the series in above converge. Moreover, we have
    
    ∑_m=1^∞-σ_m(m+1)C_2^m /τα =-C_2/2ατ(  _2F_1(1,1.5;3;-2 C_2)-   _2F_1(1.5,2;3;-2 C_2))=:H(α, T),

    where _2F_1 is a hypergeometric function
    
    _2 F_1 (a,b,c,z)=∑_k=0^∞(a)_k (b)_k/(c)_kz^k/k!,

    with
    
    (a)_k={[               1       for  k=0,; a(a+1)⋯ (a+n-1)       for  k>0. ].

    Now we are able to rewrite the main inequality of (<ref>) in the following form
    
    P(u)-P(v)_B≤  |H(α,T)| u-v_B.

    Let us notice the inverse dependence of function H(α, T) on the parameter α. Hence, choosing sufficiently large  α 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,∞].

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 L v_0=0, with appropriate initial conditions. Hence,
    
    v_0(s)=B_0cosτ s+A_0sinτ s.

    Next, we can expand the kernel G(z) in the integral in (<ref>) and get
    
    v(s)=    (B_0+1/τ∑_m=1^∞σ_m∫_T^ssin(τ t) ϕ(t)^m v(t)^m+1  d t)cosτ s
       +(A_0-1/τ∑_m=1^∞σ_m∫_T^scos(τ t) ϕ(t)^m v(t)^m+1  d t)sinτ s.

     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τ s and cosτ s. We will use the mathematical induction to prove that the function v_n can be reduced to the form 
    
    v_n(s)∼ A_n sinτ s+B_ncosτ s,  as s→∞,

    for all 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 (<ref>) we have that the function v_0 is a linear sum of sine and cosine functions. If we substitute the v_0 to (<ref>) and use a binomial expansion we get
    
    v_1=   (A_0-1/τ∑_m=1^∞ ∑_j,k
     j+k=m+1m+1jσ_m A_0^jB_0^m+1-jJ_1)sinτ s
       + (B_0+1/τ∑_m=1^∞∑_j,k
     j+k=m+1m+1jσ_m A_0^jB_0^m+1-jJ_2)cosτ s,

    where 
    
    J_1    = ∫_T^se^-λ tsin(τ t)^j cos(τ s)^k+1  d t,
    
            J_2    =  ∫_T^se^-λ tsin(τ t)^j+1cos(τ s)^k  d t.

    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
    
    J_1    = ∫_T^se^-λ tsin(τ t)^j cos(τ s)^k+1  d t=∫_0^se^-λ t(e^τ t i-e^-τ t i/2 i)^j(e^τ t i+e^-τ t i/2 )^k+1  d t
        = ∫_T^se^-λ t∑_l=0^j∑_p=0^k+1jlk+1p1/2^j+k+1i^j(-1)^l e^τ s i(j-2 l+k+1-2 p) d t
       = ∑_l=0^j∑_p=0^k+1jlk+1p1/2^j+k+1i^j(-1)^l ∫_T^se^-λ t+τ t i(j-2 l+k+1-2 p) d t
       = ∑_l=0^j∑_p=0^k+1jlk+1p1/2^j+k+1i^j(-1)^l/τ i(j-2 l+k+1-2 p)-λe^-λ t+τ t i(j-2 l+k+1-2 p)|_t=T^t=s
       ∼ -1/2^j+k+1∑_l=0^j∑_p=0^k+1jlk+1p(-1)^l+j(-λ-iτ(j-2 l+k+1-2 p))i^j e^-λ T+μ T i/(τ (j-2 l+k+1-2 p))^2+λ^2,  as   s→∞.

    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 (<ref>) 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 (<ref>) then repeating the tedious calculations performed above for function v_1 we get 
    
    v_n+1∼ A_n+1sinτ s +B_n+1cosτ s.

    Since both the base case and the inductive step have been performed,  by mathematical induction we conclude the asymptotic form of function v_n (<ref>) for arbitrary  n∈ℕ.


The asymptotic solution (<ref>)  of eq. (<ref>) 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 (<ref>) with n=0 and n=6 is depicted in  Figs. <ref> and <ref>. For simplicity we used a finite sum in formulas (<ref>) 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

    A_n+1≈   v'(T
        ) cos (T τ )/τ+v(T) sin (T τ )-1/τ∑_m=1^N ∑_j,k
     j+k=m+1m+1jσ_m A_n^jB_n^m+1-jJ_1,
    
            B_n+1≈    v(T) cos (T τ )-v'(T) sin (T τ )/τ+1/τ∑_m=1^N∑_j,k
     j+k=m+1m+1jσ_m A_n^jB_n^m+1-jJ_2,

where J_1 and J_2 are defined as before and N=30 in our analysis.



§ 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 ϵ what allow us to estimate the remainder from above by a function decreasing to zero as ϵ 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. 


plain
