Loading [MathJax]/jax/output/SVG/jax.js
Research article

Some results for a variation-inequality problem with fourth order p(x)-Kirchhoff operator arising from options on fresh agricultural products

  • Received: 26 October 2022 Revised: 20 December 2022 Accepted: 23 December 2022 Published: 09 January 2023
  • MSC : 35K99, 97M30

  • In this paper, we study variation-inequality initial-boundary value problems with fouth order p(x)-Kirchhoff operators. First, an operator is constructed based on the Leray Schauder principle, and the existence of solutions is obtained. Secondly, the stability and uniqueness of the solution are analyzed after the conditions are appropriately relaxed on the Kirchhoff operators.

    Citation: Tao Wu. Some results for a variation-inequality problem with fourth order p(x)-Kirchhoff operator arising from options on fresh agricultural products[J]. AIMS Mathematics, 2023, 8(3): 6749-6762. doi: 10.3934/math.2023343

    Related Papers:

    [1] Gunisetty Ramasekhar, Shaik Jakeer, Seethi Reddy Reddisekhar Reddy, Shalan Alkarni, Nehad Ali Shah . Biomedical importance of Casson nanofluid flow with silver and Fe2O3 nanoparticles delivered into a stenotic artery: Numerical study. AIMS Mathematics, 2024, 9(8): 23142-23157. doi: 10.3934/math.20241125
    [2] Kolade M. Owolabi, Kailash C. Patidar, Albert Shikongo . Numerical solution for a problem arising in angiogenic signalling. AIMS Mathematics, 2019, 4(1): 43-63. doi: 10.3934/Math.2019.1.43
    [3] Huiyan Peng, Xuemei Wei . Qualitative analysis of a time-delayed free boundary problem for tumor growth with Gibbs-Thomson relation in the presence of inhibitors. AIMS Mathematics, 2023, 8(9): 22354-22370. doi: 10.3934/math.20231140
    [4] Guillaume Herpe, Julien Dambrine, Inès Bennis, Clément Thomas, Stéphane Velasco, Rémy Guillevin . Separation of perfusion phases in angiographies. AIMS Mathematics, 2021, 6(1): 938-951. doi: 10.3934/math.2021056
    [5] Shabiha Naz, Tamizharasi Renganathan . An exact asymptotic solution for a non-Newtonian fluid in a generalized Couette flow subject to an inclined magnetic field and a first-order chemical reaction. AIMS Mathematics, 2024, 9(8): 20245-20270. doi: 10.3934/math.2024986
    [6] Qiongru Wu, Ling Yu, Xuezhi Li, Wei Li . Dynamic analysis of a Filippov blood glucose insulin model. AIMS Mathematics, 2024, 9(7): 18356-18373. doi: 10.3934/math.2024895
    [7] Yan Yan . Multiplicity of positive periodic solutions for a discrete impulsive blood cell production model. AIMS Mathematics, 2023, 8(11): 26515-26531. doi: 10.3934/math.20231354
    [8] M. Adel, M. M. Khader, Hijaz Ahmad, T. A. Assiri . Approximate analytical solutions for the blood ethanol concentration system and predator-prey equations by using variational iteration method. AIMS Mathematics, 2023, 8(8): 19083-19096. doi: 10.3934/math.2023974
    [9] Sumiya Nasir, Nadeem ul Hassan Awan, Fozia Bashir Farooq, Saima Parveen . Topological indices of novel drugs used in blood cancer treatment and its QSPR modeling. AIMS Mathematics, 2022, 7(7): 11829-11850. doi: 10.3934/math.2022660
    [10] J. Vanterler da C. Sousa, E. Capelas de Oliveira, L. A. Magna . Fractional calculus and the ESR test. AIMS Mathematics, 2017, 2(4): 692-705. doi: 10.3934/Math.2017.4.692
  • In this paper, we study variation-inequality initial-boundary value problems with fouth order p(x)-Kirchhoff operators. First, an operator is constructed based on the Leray Schauder principle, and the existence of solutions is obtained. Secondly, the stability and uniqueness of the solution are analyzed after the conditions are appropriately relaxed on the Kirchhoff operators.



    Arterial deformations arise in blood flow when surrounding tissue invades the space available for a blood vessel to maintain its circular cross section. This may occur in steady state when the invading tissue is pathological, or in oscillatory state when the invading tissue is driven by the effects of pulsatile blood flow.

    In the brain, the presence of a tumor may compress surrounding lymphatic and blood vessels, causing flow disruptions, especially within the restrictive environment of the rigid skull [30,34,40]. In the heart, coronary vasculature embedded within the ventricular walls undergo periodic compression and deformation with each contraction of the heart muscle [43]. Segments of the aorta near the heart have also been reported to undergo periodic deformations from circular to elliptic cross section with each heart beat [27]. Coronary arteries tethered to the surface of the heart undergo a different kind of deformation as they are laterally displaced with each heart beat, causing a lateral acceleration of fluid and a lateral force on the tube wall, resulting in a change in its shape from a circular to an elliptic cross section [7]. Flow in tubes of noncircular cross sections, both steady and pulsatile, have also been discussed in relation to the movement of spinal fluids under normal and pathological conditions [8,17,19,21,39].

    It is well known that the flow in a tube of circular cross section is singular in the sense that any departure from the circular geometry of the cross section causes a reduction in the flow rate as well as a redistribution of the shear stress along the circumference of the tube wall whereby the shear stress at some points will be higher than that in an equivalent tube of circular cross section [12]. Both of these changes are important in blood flow, the latter in particular in relation to atherosclerosis [5,14,22,36].

    While blood vessel deformation by surrounding tissue may lead to many different forms of deformation of the vessel cross section, in the present study, to keep the problem mathematically tractable, we consider the limited problem of deformations from circular to elliptic cross sections.

    Flow within a blood vessel is generally under neurovascular control whereby a change in flow rate is mediated by a change of vessel diameter. The latter in turn is mediated by a change in muscular tension within the vessel wall to the effect of changing the length of the wall circumference [35]. If the vessel is deformed by surrounding tissue such that its cross section is transformed from circular to elliptic form, two distinctly different scenarios may follow, which we shall refer to as "passive" and "active" scenarios (Figure 1). Under a passive scenario the neurovascular control is absent, and a change from circular to elliptic cross section occurs with the circumference of the vessel wall remaining constant. Under the active scenario the neurovascular control responds by changing the tension within the vessel wall in an attempt to maintain the flow rate by keeping the cross sectional area available to the flow constant. The aim of the present study is to outline the analyses associated with these two scenarios and to present results illustrating the hemodynamic consequences in the two cases.

    Figure 1.  A blood vessel of circular cross section is compressed by surrounding tissue such that its cross section becomes elliptic with semiminor axis β=fea where fe is a prescribed fraction of the circle radius a. Under a passive scenario (blue) regulatory control is absent and the length of circumference of the resulting ellipse is the same as that of the circle. Under an active scenario (red) the regulatory system intervenes in an attempt to keep the cross sectional area of the resulting ellipse the same as that of the circle.

    While from a geometrical perspective the change from circular to elliptic cross sections may seem to be a "smooth" change, from a mathematical perspective it presents a discontinuity in the character of the equations governing pulsatile flow as well as in their solutions. Specifically, in the case of circular cross sections the equations governing the flow are Bessel equations and the solutions involve Bessel functions, while in the case of elliptic cross sections the flow is governed by Mathieu equations and the solutions involve Mathieu functions [13].

    The study of pulsatile flow in tubes of elliptic cross sections has been hampered in the past because of difficulties involved in the solution of these equations and in the numerical evaluation of Mathieu functions with complex arguments [3,8,12,33,45]. In the present study we use a methodology described in [4] to overcome these difficulties and to extend the range of ellipticity at which flow properties can be evaluated. In particular, the effects of vessel deformation on flow rate and on shear stress distribution along the vessel wall are presented.

    Consider an ellipse with semi-major and semi-minor axes, α, β, respectively. Figure 2 shows ellipses in confocal elliptic ξ, η coordinates, which we will find useful. If the foci are at (±d,0) then the normal Cartesian coordinates are x=dcoshξcosη and y=dsinhξsinη. If the parameter of the outer ellipse is ξ0 then α=dcoshξ0 and β=dsinhξ0. We have d2+β2=α2 from elementary geometry. The eccentricity ε of the outermost ellipse, at ξ=ξ0, is defined by ε=d/α=sechξ0. Thus the eccentricity of the confocal ellipses changes as ξ0 changes.

    Figure 2.  Confocal elliptic coordinate system x=dcoshξcosη, y=dsinhξsinη used in the solution of the governing equations, where ξ=ξ0 is the outer circumference of a cross section of the tube of elliptic cross section. Foci at (±d,0) indicated by solid black dots. The length of the semimajor axis of the outermost ellipse is α=dcoshξ0 and the length of the semiminor axis is β=dsinhξ0.

    Using polar coordinates x=αsinθ and y=βcosθ, the circumference of the ellipse is given by

    4π/20(dxdθ)2+(dydθ)2 dθ (2.1)
    =4αE(ε) (2.2)

    where

    E(ε)=π/201ε2sin2θ dθ (2.3)

    is the complete elliptic integral of the second kind [20].

    Passive scenario  Under this scenario the change from circular to elliptic cross section occurs while keeping the length of the circumference constant. The relationship between the cross sectional areas of the circle and of the ellipse is shown in Figure 3.

    Figure 3.  Relationships between areas of the circular and the elliptic cross sections when the length of their circumferences are the same. Se and Sc are the areas of the elliptic and circular cross sections, respectively. In the limit, as the vessel cross section is flattened such that fe0, the area of the ellipse vanishes.

    For an ellipse of eccentricity ε and a circle of radius a to have the same length of circumference, we have

    2πa=4αE(ε) (2.4)

    therefore

    αa=π2E(ε) (2.5)
    βa=1ε2π2E(ε) (2.6)

    If the area of the ellipse is denoted by Se (=παβ) and the area of the circle is denoted by Sc (=πa2), then the ratio of the two is given by

    SeSc=(π2E(ε))21ε2 (2.7)

    In the passive scenario, where the circumference remains constant on deformation from a circle of radius a, the foci of the ellipse are located at (±d,0) where

    d=πε2E(ε)a. (2.8)

    This tends to πa/2 as ε tends to 1.

    Compressing a circle of original radius a to an ellipse with semi-minor axis β=fea with fe<1 while keeping the circumference 2πa=4E(ε)α constant requires that α=gea where ge>1 is given by the following implicit formulae from Eqs (2.5)–(2.6):

    fe=π1ε22E(ε) (2.9)
    ge=π2E(ε). (2.10)

    To find ge for a given fe one must solve the transcendental Eq (2.9) for ε, and then use that in the equation for ge. This is straightforward in Maple, by use of the command fsolve. For convenience, we tabulate some fractions in Table 1.

    Table 1.  Table of the major and minor axes, α, β, of an ellipse having the same length of circumference of a circle of radius a, where ε is the eccentricity of the ellipse.
    fe=β/a ε ge=α/a
    0.4 0.9611 1.448
    0.5 0.9334 1.392
    0.6 0.8925 1.330
    0.7 0.8314 1.260
    0.8 0.7359 1.182
    0.9 0.5698 1.096

     | Show Table
    DownLoad: CSV

    Active scenario  Under this scenario the change from circular to elliptic cross section occurs while keeping the cross sectional area constant. In the active scenario, compressing a circle of radius a so that its semi-minor axis β=fea is a given fraction fe of the original radius, then since the area is παβ we must have α=a/fe.

    The ratio of the circumference of an ellipse to the circumference of a circle with the same area is

    CeCc=2E(ε)π(1ε2)1/4. (2.11)

    This is plotted in Figure 4.

    Figure 4.  The ratio of the circumference of an ellipse to the circumference of a circle with the same area. The formula is Ce/Cc=2E(ε)/(π(1ε2)1/4), with E(ε) being the complete elliptic integral. As the fraction fe0, the eccentricity ε1. The ratio of circumferences is singular as the fraction goes to 0, as one would expect. The figure thus points to an intrinsic limitation of the regulatory system to increase the length of circumference of a blood vessel of elliptic cross section in an attempt to maintain its cross sectional area.

    In the active scenario, where on deformation from a circle of radius a the circumference is stretched by the regulatory system in order to keep the area constant, the foci of the ellipse are located at (±d,0) where

    d=ε(1ε2)1/4a. (2.12)

    This is unbounded as ε1.

    The properties of steady flow in a tube of elliptic cross section will be used as reference for the corresponding properties in pulsatile flow. The function governing the axial velocity u0,e is given by [44]:

    u0,e=k0α2β22μ(α2+β2)(1x2α2y2β2) (2.13)

    where μ is viscosity, and k0 is the constant pressure gradient driving the flow. The maximum velocity occurs at x=0,y=0, the center of the ellipse:

    ˆu0,e=k0α2β22μ(α2+β2) (2.14)

    and volumetric flow rate is given by [44]

    q0,e=ˆu0,eSe2 (2.15)

    Shear stress on the tube wall is given by [31]

    τ0,e(x,y)=k0α2β2α2+β2(x2α4+y2β4)1/2 (2.16)

    Maximum shear occurs at the ends of the minor axis

    ˆτ0,e=k0α2βα2+β2 (2.17)

    Minimum shear occurs at the ends of the major axis

    ˇτ0,e=k0αβ2α2+β2 (2.18)

    Maximum velocity and maximum shear are related by

    ˆu0,e=β2μˆτ0,e (2.19)

    The corresponding quantities for a circle are obtained by letting βα or equivalently fe1. Then, for instance, the maximum and minimum shear both become k0a/2.

    In Figures 57 we use these formulas to compare the difference in steady flow in the two different scenarios, active (depicted with red curves in the figures) and passive (depicted with black curves). All of the quantities above involve the semimajor and semiminor axes, α and β. While for both scenarios β=fea is the same, the value of α will be different in the active scenario (a/fe) to the passive scenario gea where ge is computed by solving a transcendental equation. We see that there is indeed some difference in the flow quantities that arises in the two scenarios.

    Figure 5.  Maximum velocity in steady flow in a tube of elliptic cross section from Eq (2.14) compared to that in a tube of circular cross section in the two scenarios, active (red) and passive (black).
    Figure 6.  Maximum wall shear stress in steady flow in a tube of elliptic cross section from Eq (2.17) scaled by the constant shear stress on boundary of a tube of circular cross section and compared under the active (red) and passive (black) regulatory scenarios.
    Figure 7.  Flow rate in steady flow in a tube of elliptic cross section from Eq (2.15) scaled by the flow rate in a tube of circular cross section and compared under the active (red) and passive (black) regulatory scenarios.

    The axial velocity, ue, in a tube of elliptic cross section can be written as the sum of a steady part, u0,e, and an oscillatory part, uϕ,e,

    ue(x,y,t)=u0,e(x,y)+uϕ,e(x,y,t) (2.20)

    The equation governing the oscillatory part of the velocity is given by [13]

    uϕ,et+1ρpz=μρ(2uϕ,ex2+2uϕ,ey2) (2.21)

    The solution of this equation is facilitated by changing to the confocal elliptic coordinates shown in Figure 2. These were introduced by Lamé, who called them "thermometric coordinates"[23].

    x=dcoshξcosη, y=dsinhξsinη (2.22)

    where the foci are at (±d,0) and ξ, η are the elliptic coordinates. Using these confocal elliptic coordinates, an oscillatory pressure gradient of the form

    pz=k0eiωt, (2.23)

    and separation of variables

    uϕ,e(ξ,η,t)=w(ξ,η)eiωt, (2.24)

    Equation (2.21) can be formulated as an inhomogeneous Helmholtz equation

    2d2(cosh2ξcos2η)(2wξ2+2wη2)iρωμw=k0μ. (2.25)

    Using the translation

    w(ξ,η)=v(ξ,η)k0iρω, (2.26)

    the inhomogeneous term of Eq (2.25) is eliminated and the equation becomes

    (2vξ2+2vη2)i2Λe(cosh2ξcos2η)v=0, (2.27)

    where

    Λe=ρωd2μ (2.28)

    is a nondimensional frequency parameter.

    The boundary conditions are given by

    v(ξ0,η)=k0iρω (no slip at tube wall) (2.29)
    vξ|ξ=0=0 (symmetry) (2.30)
    v(ξ,0)=v(ξ,π) (π periodic in η) (2.31)

    While direct numerical solution of the governing equation (2.27) is also possible, in this paper we pursue a solution based on the use of separation of variables, leading to the use of Mathieu functions. We do this in order to maintain the analytical connection with the classical solution of pulsatile flow in tubes of circular cross sections based on Bessel functions [44]. The use of Mathieu functions is not as straightforward as the use of Bessel functions, however, in part because of numerical difficulties in the evaluation of Mathieu functions of imaginary arguments. Balancing that, this method is spectrally accurate, and does not require many eigenfunctions for the range of q that we consider here. Typically, we need only terms up to about N=6 or N=8. We postpone discussion of how to construct and evaluate the solution until section 3.

    In detail, the method proceeds as follows. The treatment is standard, and we include it mostly for notation and readability. Applying separation of variables to Eq (2.27) then yields two separate equations:

    d2gdη2+(s2qcos2η)g=0 (2.32)
    d2fdξ2(s2qcosh2ξ)f=0 (2.33)

    where s is a separating constant and

    q=iΛe4. (2.34)

    There is some risk of notational confusion because flow rate is often denoted by the variable q; here we will use flow variables with subscripts only, and the undecorated symbol q will refer to the parameter in Eq (2.34). This notation is standard for Mathieu functions, and we believe less confusion results when we use symbols in this fashion.

    Equation (2.32) is the Mathieu equation, and Eq (2.32) is the modified Mathieu equation. These equations are equivalent, with the change of variable ξ=iη. The character of the solutions of the two equations are quite different, however.

    Since in the present problem η varies from 0 to 2π, then η must have periodicity π or 2π, which only occurs for discrete values of s, the eigenvalues of the Mathieu equation, designated by sm in [28]. These eigenvalues are more commonly denoted nowadays with the letters am and bm; see the DLMF https://dlmf.nist.gov/28.2.ii. This results in the set of Mathieu eigenfunctions cem and sem for Eq (2.32). Both cem and sem are periodic, and cem is even, whereas sem is odd. By convention, if m is odd, then the Mathieu functions have period 2π, while if m is even, the Mathieu functions have period π. In our problem, we are only interested in even m values, so as to have π–periodicity, and even functions to satisfy symmetry along both axes*.

    *Symmetry breaking might very well be possible in a physical situation, and we believe it will be worthwhile to investigate this in future work

    The even π-periodic solutions of Eq (2.32) are the ordinary Mathieu functions denoted by ce2m(η,q). We will also need the modified Mathieu functions for the same m and the same value of q, which are solutions of Eq (2.32). The solution that we will compute will then be of the form

    v(ξ,η)=k0iρωm0b2mCe2m(ξ,q)ce2m(η,q) (2.35)

    where the coefficients b2m will be determined by the no-slip boundary conditions, and we have taken the opportunity for a convenient scaling by k0/(iρω).

    For certain values of q, however, such as the Mulholland–Goldstein value q1.4688i (see [4]), the Mathieu equation has double eigenvalues and at those points special care must be taken, because the ordinary Mathieu functions no longer form a complete set of orthogonal functions for expansion. As q tends to the Mulholland–Goldstein point, ce0(η) and ce2(η) coalesce and become the same function, and to ensure that expansion in these eigenfunctions is possible (also known as "completeness"), a generalized eigenfunction must be added to the set of Mathieu functions. In practice, as we will see, these isolated points make little difference to the solution because the overall problem is continuous (indeed analytic) in q, and so it is only the solution process which must be altered at these points. Again, this is discussed in [4].

    Using the no slip boundary condition from Eq (2.29), we find

    v(ξ0,η)=k0iρω=k0iρωm0b2mCe2m(ξ0,q)ce2m(η,q). (2.36)

    Since the Mathieu functions are orthogonal under the bilinear form

    f,g:=2πη=0f(η)g(η)dη (2.37)

    (and if they have period π, the upper limit on the integral can be reduced to π), then multiplying Eq (2.36) by ce2p(η,q) and integrating with respect to η gives

    2π0ce2p(η,q) dη=b2pCe2p(ξ0,q)2π0ce22p(η,q) dη. (2.38)

    We note that the bilinear form does not involve the complex conjugate. Eigenvalues need not be real, and as parameters vary, eigenfunctions can coalesce. Expansion in Mathieu functions is similar to harmonic expansion, but more complicated. In the usual case, when eigenvalues are simple, b2m is given by

    b2m=2π0ce2m(η,q) dηCe2m(ξ0,q)I2m. (2.39)

    Here

    I2m=2π0ce22m(η,q) dη. (2.40)

    We compute these integrals by doing exact integration of the polynomial "blends" interpolating the solution, as described in section 3. The computational cost for this is trivial.

    Remark 2.1. The integral I2m can be zero. In particular, if q=1.4688i (the Mulholland–Goldstein point mentioned earlier) then this integral is zero. In this case, the expansion must be computed by a different method. We ignore this possibility for the moment.

    Remark 2.2. The value of Ce2m(ξ0,q) might a priori be zero. In this case, we would have found a natural frequency of oscillation, and the solution would exhibit resonance. We did not encounter resonance in any of the configurations we tried. It seems that symmetric, even Mathieu functions with purely imaginary values of q have no zeros on the imaginary axis, although we have not proved this.

    With this nonzero integral, Equation (2.26) for w becomes

    w(ξ,η)=k0iρω(1m0b2mCe2m(ξ,q)ce2m(η,q)). (2.41)

    The oscillatory flow velocity in a tube of elliptic cross section is then [13]

    uϕ,e(ξ,η,t)=4ˆu0,eiλe(1m0b2mCe2m(ξ,q)ce2m(η,q))eiωt (2.42)

    where

    λe=12sinh2ξ0tanh2ξ0Λe=2(1ε2)ε2(2ε2)Λe (2.43)

    is a second nondimensional frequency parameter.

    The flow rate is obtained by integrating the oscillatory velocity over the elliptic cross section

    qϕ,e(t)=Duϕ,e(ξ,η,t) dA (2.44)
    =(Dv(ξ,η) dA+4ˆu0,eiλeD dA)eiωt (2.45)

    where D is the region enclosed by the bounding ellipse. The second integral on the right side of Eq (2.45) can be evaluated analytically:

    4ˆu0,eiλeD dA=8q0,eiλe (2.46)

    where q0,e is the steady flow rate in a tube of elliptic cross section (Eq (2.15)). The first integral on the right hand side of Eq (2.45) is then evaluated using v in Eq (2.27)

    Dv(ξ,η) dA=μiρωD2d2(cosh2ξcos2η)(2vξ2+2vη2) dA (2.47)

    If n is an outward pointing normal and ds is an elemental surface, then by Green's theorem it follows that

    μiρωD2d2(cosh2ξcos2η)(2vξ2+2vη2) dA=μiρωDvn ds (2.48)

    where D is the positively oriented bounding curve of D.

    It is shown in McLachlan [24] that ds=δdη and dn=δdξ, where

    δ=d(cosh2ξcos2η)1/2 (2.49)

    Thus, Equation (2.48) becomes

    μiρωDvn ds=μiρω2π0(vξ)ξ=ξ0dη (2.50)

    and

    (vξ)ξ=ξ0=k0iρωm0b2mCe2m(ξ0,q)ce2m(η,q) (2.51)

    where here denotes differentiation with respect to ξ. Because blends are polynomials, differentiation with them is simple, and the code we use provides for this automatically. We thus get (apart from rounding errors) exact derivatives of the interpolants being used to represent the solutions. Because the solutions are so high-order, the derivatives are themselves accurate: while they typically lose an order of accuracy for each derivative taken, if one starts with order 16 then taking one derivative does not do much harm. We remark that with high enough frequency, however, which does occur with large eigenvalues for Mathieu functions, one would need to work to higher precision to maintain this accuracy. For the computations of this paper, we only used higher precision to check the numerics, and found double precision to be perfectly satisfactory.

    It is important to remember that because the code implements Ce2m(ξ,q) by ce2m(iξ,q) one has to use the chain rule and multiply by i: Ce2m(ξ,q)=ice2m(iξ,q).

    Integration of this formula with respect to η is straightforward, using the exact quadrature formula for blendstrings. But in fact we have already integrated each of these functions, in computing the b2m. NB: if the integrals were only to π and not to 2π, one must multiply the following formula by 2.

    Using Eq (2.40) we find that

    μiρω2π0(vξ)ξ=ξ0dη=μk0ρ2ω2m0b22mI2mCe2m(ξ0,q)Ce2m(ξ0,q). (2.52)

    If we further use the relation

    λeμρωSe=1πtanh2ξ0, (2.53)

    this implies that the oscillatory flow rate in a tube of elliptic cross section is given by the following (cf. [13]):

    qϕ,e(t)=8q0,eiλe(11iπλetanh2ξ0m0b22mI2mCe2m(ξ0,q)Ce2m(ξ0,q))eiωt. (2.54)

    By its definition, the wall shear stress is given by

    τϕ,e(η,t)=μ(uϕ,en)D=μ(vn)Deiωt (2.55)

    where uϕ,e is the oscillatory velocity in a tube of elliptic cross section (Eq (2.42)). Using the elemental arc length analysis of McLachlan [24], it can be shown that

    (vn)D=1δ0(vξ)ξ=ξ0 (2.56)

    where δ0 is

    δ0=d(cosh2ξ0cos2η)1/2 (2.57)

    Substituting from Eq (2.51) for the derivative on the right hand side, this becomes

    (vn)D=1δ0(k0iρωm0b2mCe2m(ξ0,q)ce2m(η,q)). (2.58)

    Using Eq (2.17) we can replace k0 by ˆτ0,e(α2+β2)/(α2β), or more conveniently by the limiting case of the circle: k0=2ˆτ0,c/a. Remember that a is the radius of the original circle, and β=fea. After some algebra we obtain the following expression for oscillatory wall shear stress in a tube of elliptic cross section:

    τϕ,e(η,t)=4βfeˆτ0,ciδ0λe(2ε2)(m0b2mCe2m(ξ0,q)ce2m(η,q))eiωt (2.59)

    For reference and comparison, the (constant) oscillatory wall shear stress in a tube of circular cross section is given by the following [42]:

    τϕ,c=2τ0,cΛcJ1(Λc)J0(Λc) (2.60)

    where Jk(z) for k=0, 1 are Bessel functions of the first kind, τ0,c=k0a/2, and

    Λc=(i12)ρωμa. (2.61)

    We will not review all existing numerical methods for computing with Mathieu functions here, but instead refer to [4], which is available as an open-access article. We will, however, summarize the method that we actually used, and give a few more details about the method in a subsection that may be skipped by a reader more concerned with the results, as opposed to how we got them.

    For notational convenience we refer to values of q with positive imaginary part, but because the eigenvalues are the same for q and q in the even and symmetric case (see eg. the DLMF https://dlmf.nist.gov/28.2), this is sufficient for our application (which has negative imaginary part) and saves writing many minus signs.

    The previous work of Haslam and Zamir in [12] used truncations of an infinite tri-diagonal eigenvalue-eigenvector problem to obtain approximations to the eigenvalues a2m. This method goes back at least to the work of Ince, and is widely used [4]. The matrix in question, for the even and symmetric eigenfunctions, is

    [02q0002q4q000q16q000q36q000q64][2A0A2A4A6A8]=λ[2A0A2A4A6A8]. (3.1)

    Truncation at "large enough" dimension gives good estimates of the eigenvalues, but there is a question of exactly how large should we take the matrix, and once the eigenvalues have been computed, how accurate they are. Notice that this is a complex symmetric matrix, not a Hermitian matrix.

    In our computations, we start with the matrix method, but only to get initial estimates of the eigenvalues a2m. We then apply the continued fraction method of Blanch as described in [4] and use Newton's method to refine the eigenvalues to the desired accuracy. This tells us precisely how accurate each eigenvalue is, and is more efficient than computing larger and larger matrices until the eigenvalues converge. Our procedure works well enough for all simple eigenvalues, although sometimes we have to increase precision. For the double eigenvalues, we proceed differently.

    Double eigenvalues occur for purely imaginary q, but (as elsewhere in the complex q-plane) only at isolated points: At the Mulholland–Goldstein point q1.4688i, and (next smallest) q16.47i, and so on. We have pre-computed several of these by the method of Hunter and Guerrieri [15]. They are tabulated in [4] and are also available on-line in the code repository for that paper.

    Given numerical values for the semimajor axis α and semiminor axis β, and given a numerical value for the (purely imaginary) parameter q, we computed up to certain index N (frequently N was 6 and sometimes 8; because this is a spectral method, convergence is very rapid) of the Mathieu eigenvalues a2m(q), for m=0, 1, , N. If the eigenvalues were distinct (which was usually, but not always, the case) then we computed the Mathieu functions ce2m(η) on the interval 0ηπ and the corresponding modified Mathieu functions Ce2m(ξ) on the interval 0ξξ0=invcosh(α/d)=invsech(d/α)=invsech(ε). This is because α=dcoshξ0 or, more simply, ε=sechξ0 gives the value ξ0 of the parameter ξ at the tube wall.

    Here, we are using David Jeffrey's notation for functional inverses: y=invcosh(x) means x=cosh(y), etc. This notation is superior for branched inverses, and superior pedagogically even for simple functions, to the more common overloading of superscripts or use of the inappropriate word "arc", and we hope that it catches on.

    To compute the Mathieu functions and modified Mathieu functions, we used the Hermite-Obreshkov integrator sketched in [4]. We worked in double precision (except where noted explicitly here) and typically used an order 30 or 40 method, with grade§ 15 or 20 Taylor series computed on each marching step, and "blendstrings" as piecewise polynomial interpolants giving the value of the solution (and whatever derivatives were required).

    §The word "grade" means "degree at most". This is convenient because the final Taylor coefficients computed might be zero, but this is still useful information.

    We treat the Mathieu equation (and the modified Mathieu equation) as an initial-value problem (IVP) for an ordinary differential equation (ODE), once both q and the eigenvalue s are fixed. To compute the Mathieu function, we could use almost any standard method to solve the IVP. But the modified Mathieu equation is related to the Mathieu equation by the change of variables ξ=iη. That is, if the standard method chosen for the Mathieu equation could work in the complex plane, then it could also be used for the modified Mathieu equation. This idea restricts us to implementations that work over the complex plane, but because we have a complex parameter q (in fact, purely imaginary in our application), this is necessary anyway.

    We reassure the reader that we do know and highly value the standard general methods, as described for instance in the classic [10,11]. We are also aware of the truly remarkable advances made since then, such as are described in [32]. We have even contributed to the literature and the software ecosystem in the past [37]. But while writing a special-purpose solver for the Mathieu equation—when so many good solvers already exist—might seem quixotic, bear with us for a bit: it turns out to be useful and we believe interesting, and in particular it is reassuring to have the ability to retrospectively measure how accurate the solutions are.

    Also, there is an opportunity for greater efficiency and control. Since the Mathieu equation is linear, special-purpose methods appropriate for linear problems might be used. More, since the Mathieu equation can be written in a "D-finite" or "holonomic" form, Taylor series coefficients can be computed rapidly given the initial values y(ηn) and y(ηn). In fact, we do not use the D-finite form even though it does offer the potential of significant speed-up [26]; this might be pursued in future. Straightforward generation of Taylor coefficients by Cauchy convolution with those of cos2η was fast enough for our purposes.

    This fact was already known to Mathieu, although the names D-finite or holonomic had not been invented yet in 1868. But writing the differential equation in this form allows for faster human computation, too.

    We now explain the interpolants that we use. "Blends", or two-point Hermite interpolants, are described in [6]. In brief, if one knows Taylor coefficients pj for 0jm at one end zk of an interval, and Taylor coefficients qj for 0jn at the other end zk+1 of an interval, and z=zk+sh where h=zk+1zk is the width of the interval so that 0s1, then the following polynomial "blends" the two sets of Taylor coefficients together to form an excellent approximation of the function over the interval: (Hermite, Cours d'Analyse 1873)

    Hm,n(s)=mj=0mjk=0(n+kk)sk+j(1s)n+1pj+nj=0(1)jnjk=0(m+kk)sm+1(1s)k+jqj (3.2)

    has H(j)(0)/j!=pj for 0jm and Hj(1)/j!=qj for 0jn. In this formula, differentiation is with respect to s, and care must be taken to include the correct factors of h from the chain rule when using the formula for the interval [zk,zk+1].

    The error in Hermite interpolation is known; the results on the real line are given in [18] (and the complex results were known to Hermite). Here, the general real results simplify to

    f(s)Hm,n(s)=f(m+n+2)(θ)(m+n+2)!sm+1(s1)n+1 (3.3)

    for some θ=θ(s) between 0 and 1.

    If we have a sequence of nodes, say zk for 0kM, where Taylor coefficients for an analytic function f(z) are known up to grade (say) mk at each node, then it is natural to approximate f(z) on each segment from z=zk to z=zk+1 by the blend determined by those two sets of Taylor coefficients. This gives a piecewise polynomial interpolant, which we call a "blendstring" for short.

    For instance, if Taylor series of only grade 1 are used at each node, then the blendstring is just the more familiar pure piecewise cubic Hermite interpolant on each subinterval, and the result is similar to a cubic spline. Taylor series of only grade 1 do not give us the needed accuracy, though, and we always use much higher order.

    As described in [6], these interpolants are remarkably stable numerically, even for ludicrously high order such as m=500, when implemented in a doubly-recursive Horner form. This turns out to be quite convenient for this application, where we typically use grades of 15 or so but sometimes as high as 40.

    Blends can be integrated exactly, as follows, and this is useful [6]:

    1s=0Hm,n(s)ds=(m+1)!(m+n+2)!mj=0(n+mj+1)!(j+1)(mj)!pj+(n+1)!(m+n+2)!nj=0(n+mj+1)!(j+1)(nj)!(1)jqj. (3.4)

    The numbers showing up in this formula turn out to be smaller for the higher-order Taylor coefficients, as one would expect. Note that the above formula gives (in exact arithmetic) the exact integral of the blend over the whole interval. If the blend is approximating a function f(s), then integrating equation (3.3) gives us

    1s=0f(s)dsF(1)=(1)n+1(m+1)!(n+1)!(m+n+3)!f(m+n+2)(c)(m+n+2)! (3.5)

    where, using the Mean Value Theorem for integrals and the fact that sm+1(1s)n+1 is of one sign on the interval, we replace the evaluation of the derivative at one unknown point θ with another unknown point c on the interval.

    Indeed, as described in [6], one can construct a new blendstring H(z) for the antiderivative F(z) from a blendstring h(z) for f(z), so that H(z)=h(z) exactly (up to roundoff error), and well approximates the antiderivative of f(z). This is useful for the problem at hand.

    The code is available at https://github.com/rcorless/Puiseux-series-Mathieu-double-points in the files ActiveLoopc1p0.maple for the simple eigenvalue case and ActiveDoubles.maple for the double eigenvalue case.

    We chose an implicit marching method based on Taylor series generation**, quite standard in outline, as follows. Taylor series coefficients that have been generated at the current node, say ηn, are supposed to be "known". Specifically, suppose to start with that we have generated a Taylor polynomial of grade m for our desired solution at this point.

    **Taylor series methods for solving IVP for ODE have historically been considered impractical by many people, but in fact this is not so, especially if the series coefficients can be generated easily, as in this case. The quality of the free interpolants that one gets turns out to be a significant benefit. Taylor series methods have other benefits as well: see [29] and its references.

    Suppose also that we have chosen a tentative next node, ηn+1=ηn+h. If our variable η were time, this would be a time step. The stepsize h is tentative at this point. We now generate Taylor coefficients for two independent solutions, satisfying (for one solution)

    y(ηn+1)=1 and y(ηn+1)=0 (3.6)

    and (for the complementary solution)

    y(ηn+1)=0 and y(ηn+1)=1. (3.7)

    Next, we blend the known coefficients at ηn with these independent solutions in the following way. Form a blend of the known coefficients at ηn with the zero Taylor series at ηn+1. Call the result L(η). Form a blend of the first series above at ηn+1 with the zero Taylor series at ηn. Call the result C(η). Form a blend of the second series above with the zero Taylor series at ηn and call the result S(η). Our desired solution will then be a linear combination of these three: say y=AC(η)+BS(η)+L(η). This uses the linearity of the equation, and the linear dependence of blends on their constituent Taylor coefficients.

    We then use collocation at the two points ηn+h/4 and ηn+3h/4 (which are Chebyshev–Lobatto points, not that it matters much at this low order) to give us two equations in the two unknowns A and B. That is, we compute the residuals

    rL(η):=L"+(s2qcos2η)LrC(η):=C"+(s2qcos2η)CrS(η):=S"+(s2qcos2η)S (3.8)

    at those two points, and set the residual for y to zero at those two points:

    0=ArC(ηn+h/4)+BrS(ηn+h/4)+rL(ηn+h/4)0=ArC(ηn+3h/4)+BrS(ηn+3h/4)+rL(ηn+3h/4). (3.9)

    We solve this two-by-two linear system by the exact formula for the inverse (this is as good a method as any, for such a small system) to acquire the coefficients A and B. This system is nonsingular because the solutions are linearly independent at the right endpoint and are well-scaled and well-conditioned in practice, as we observed experimentally.

    Collocation is a well-understood technique for boundary-value problems for ODE [1,2], but it has historically been used successfully for stiff initial-value problems as well [41].

    After having computed A and B and used them to form our tentative solution y, we then sample the residual of y, namely y"+(s2qcos2η)y, at the midpoint ηn+1/2=ηn+h/2. This is (asymptotically as h0) the location of the maximum residual over the step. If this is smaller than our tolerance, we accept the step and continue. Note that if the step is accepted, the Taylor coefficients then become known at ηn+1, being simply the known linear combination of the first and second sets of computed series coefficients. We also use the measured residual (by known step-size control techniques [9]) to predict the next step size hn+1 and thus ηn+2.

    If the step is rejected instead because it does not satisfy the accuracy tolerance, we reduce the stepsize by an amount indicated by the size of the measured residual (taking the order 2m into account), and try again.

    Various known heuristics and safety factors are included in order to be cautious about various contingencies (for instance, the measured residual might be accidentally small, which throws the predicted stepsize off; similarly, the stepsize predictions are determined by assuming that the derivatives involved in the error coefficients "don't change much" from step to step, but this is sometimes violated in practice). Error messages can be generated if too many stepsize reductions are encountered, or if the solver can't find a good starting stepsize††, or if the maximum number of steps is reached, as is usual with IVP solvers.

    ††We start with a pure Taylor series to estimate the initial step size h. This has some potential to go wrong, and sometimes does, because it does not benefit from implicitness, but we have found it satisfactory.

    The reasons we do this, instead of using a more standard method that has already been implemented and tested, include the following.

    1. We work from the beginning over the complex plane (most standard implementations put integration over the real line first).

    2. We can handle the double-eigenvalue case in a straightforward way. To be fair, other methods can also handle this case in a straightforward way, as well, but at least we are not at a disadvantage.

    3. The functions are entire, and therefore Taylor series are defined everywhere for them. Since blendstrings are very smooth (with grade m Taylor coefficients at each knot, they are m times continuously differentiable) they may be expected to be accurate and convenient.

    4. The problem is linear, so the implicitness of the method is simple to deal with (and there are no convergence issues in solving nonlinear equations at each step).

    5. Putting η=ηn+sh, the residual has the error expression‡‡

    ‡‡To show this, notice that the residual is O((zαk)m1) at the left endpoint (not O((zαk)m+1) because we have differentiated twice), is O((zαk+1)m1) at the right endpoint, and vanishes at the Chebyshev–Lobatto points in between.

    r(η)=Kh2msm1(s14)(s34)(1s)m1+O(h2m+1) (3.10)

    for some "constant" K depending on high-order derivatives of the solution, evaluated at some point in the interval. In comparison to an explicit Taylor series method, this gains a factor of 22m+2 in accuracy because the maximum value of the polynomial in s is 22m2. Since we typically take m=15 or higher, this accuracy gain is noticeable.

    6. The effect of the residual on the solution can be analyzed by using the Green's function for the Mathieu equation, which can easily be computed by the same methods:

    G(η,τ)=wI(η)wII(τ)wI(τ)wII(η). (3.11)

    Here we use the notation for the basic solutions as described in the DLMF https://dlmf.nist.gov/28.2.ii. The change in solution produced by a residual r(η) is

    ητ=0G(η,τ)r(τ)dτ. (3.12)

    7. We can re-use standard stepsize heuristics, which are well-known to produce "good" meshes which reflect dynamic changes in the solution.

    8. The Mathieu equation is not "stiff" with the stepsizes and tolerances we are using [38], but is rather oscillatory, and as such benefits somewhat from the implicitness of this method. There is still a stability restriction, but it is not very important compared to the stepsize restriction needed for accuracy, and this implicit method does perform better than a pure explicit Taylor series method.

    9. Using a residual (defect) control is useful even for unstable differential equations. The modified Mathieu equation can be very unstable, exhibiting doubly exponential growth.

    10. These Taylor coefficients are very easy to generate, and the code is quite simple. The fact that the order of accuracy can be chosen more or less arbitrarily is an advantage for very high-precision computation: the cost for accurate solution is polynomial in the number of bits of accuracy [16].

    11. We do want high-precision computation, because we want to be able to state unequivocally that numerical artifacts are not present, and to verify that any given solution is as accurate as the code claims. This is not a given, without an external check, because of the heuristics and safety factors needed in practice for the solver.

    12. Taking derivatives and integrals of blends is very simple, and both of these are needed for subsequent computations with the solution. We are not just interested in the solution, but also in integrals and derivatives of the solution.

    13. It might be true that this method is useful for the numerical solution of some other, similar, equations. In particular this might be of interest for D-finite (holonomic) systems. This application provides a useful test case.

    Because each computed Mathieu function and modified Mathieu function is a smooth piecewise polynomial, it can be differentiated and substituted back into the differential equation. What is left over is sometimes called the "defect" but the more usual name in numerical analysis is the residual. The solutions always had a residual comparable to the tolerance with which the solver was called; typically about 1011 if we were working in double precision, and about 1028 if we were working in 30 decimal digits. This is, of course, not enough to say that the forward error is small: one needs also to compute the Green's function, or otherwise verify that the condition number* is small.

    *By this we do not mean the condition number of a matrix (there are no matrices here) but rather the condition number of the Mathieu differential equation, which since the equation is linear, is equivalent to the maximum value of the Green's function.

    It turns out that for fe near 1, i.e. when the ellipse is nearly circular and the eccentricity ε is small, then ξ0 gets modestly large—and because the modified Mathieu functions grow doubly exponentially, the Green's function does indeed amplify errors in this case. Indeed, the condition number of Ce0(ξ,q) for evaluation, namely C=ξCe0(ξ,q)/Ce0(ξ,q) always grows exponentially with ξ. For q=0.3823i, a typical value of the parameter, the value of the condition number C(4.0) is approximately 120. This is tolerable. For all values of the parameters that we used, except for the stress test when fe=0.9999 (more about this, below), the condition number was similarly modest.

    Another way to see this is to vary the parameters (such as fe or ω) and verify that the solution does not change much. A third way is to do the computations again in higher precision. We found in the end that our computation of the Mathieu functions and modified Mathieu functions was very reliable.

    It is also possible to verify that the underlying PDE is satisfied: one computes (for instance) 2w/ξ2, 2w/η2, and w at one or several or a great many points, and substitues these values back into the partial differential equation (2.25). When we do this we see that what is left over is about the size of the integration tolerance (typically quite near to the unit roundoff level in our runs), over the whole ellipse. Figure 8a shows the result of doing this in one case, for v(ξ,η) using equation (2.27). This kind of a posteriori solution validation is a powerful check against numerical errors. What we have proved by this a posteriori computation is that we have computed the exact solution of a perturbed PDE, where the perturbation is smaller than 1010. Compared to modelling errors (for instance that the true deformed shape is not exactly elliptical, or the even greater modelling error of neglecting the third dimension) this shows definitively that the numerical method has performed satisfactorily. This is a useful guard against blunders, as well: We were reminded to use the chain rule, and also found a typographical error in one equation, when we did this.

    Figure 8.  Guarding against blunders and errors: On the left, the computed residual δ=2v/ξ2+2v/η2iρωd2v/(2μ) with N=6 terms, f=0.6, c=1.0cm, active scenario, ω=3.83722019332829 (a double eigenvalue case), in double precision. This shows that the differential equation is satisfied to better than 1010. On the right, we plot our computed ϵ=k0(v(ξ0,η)1)/(iρω), which ideally should be zero at the boundary, for the same parameter values, real part in black and imaginary part in red, on 0η2π. We see the effect of truncating our series expansion at N=6. Because this is a double eigenvalue case and the Green's function was used, the piecewise polynomial approximation for the solution is not quite periodic: there is a tiny jump between the values at η=0 and at η=2π. Indeed the error is not quite periodic with period π, which it would be ideally.

    Finally one needs to check the boundary conditions. In Figure 8b we see one such check. The oscillatory nature of the error indicates that it is truncation error we are seeing—the effect of taking only N=6 terms in the expansion. With a high enough N, one sees only rounding errors at this stage.

    If q is small, then the continued fraction approach of Blanch becomes somewhat fragile. Blanch had performed a good numerical analysis of the method, and using her methods it can be made to work well in this situation by various adjustments. In contrast, however, the performance of the matrix method improves as q0, so it is simpler just to drop the use of the continued fraction approach when q is "too small". We chose, somewhat arbitrarily, to use just the matrix method if |q|<0.2.

    The solver is meant for use by people willing to adjust parameters experimentally and not, in Blanch's words, simply to be run "in a robot-like fashion." The code has, in particular, an aggressive initial step-size heuristic based on explicit Taylor series. This got into trouble for some of the runs in the passive scenario with circumference 2.0cm and fractions fe closer to 1. We could have adjusted the parameters (tolerance and grade m of Taylor approximation) but it was simpler to use higher precision for those runs, which we did in 30 decimal digits. The time penalty was slight, even though we have not optimized the code for speed. The first step for that, of course, would be to use a production language instead of a prototyping language such as Maple (which saves our time and not the computer's). Even so, the solver is gratifyingly rapid, even at very high precision.

    The only real difficulty that occurs with expansion in Mathieu functions is when the fraction fe is nearly 1. That is, the nearly-circular case is the difficult one for expansion in terms of Mathieu functions. This is because the coordinate transformation used, namely x=dcoshξcosη and y=dsinhξsinη, becomes singular as the focal distance d0, which it must as the ellipse becomes a circle. Another way to think about this is to "zoom out" on confocal ellipses; the larger the diameter, the more nearly circular the confocal ellipses are.

    This singular behaviour shows up in several ways, numerically. For instance, taking fe=0.999 in the active scenario, and choosing ω so q is the Mulholland–Goldstein point, seems that it should not cause problems. But it does, because ξ0=invsech(ε) is then about 4.6. Although that does not seem like much, the modified Mathieu functions grow doubly exponentially; in this case, Ce0(ξ0,q) has magnitude 6.09721045. The modified Mathieu equation becomes very difficult to integrate accurately for large ξ because of this doubly exponential growth.

    One is used to exponential growth, but doubly exponential growth is remarkably difficult to deal with. To see the asymptotics for Ce(ξ,q) look at the DLMF. See https://dlmf.nist.gov/28.25 in particular.

    As a stress-test for the code, we solved the problem at very high precision with fe=0.9999 in both the active and passive scenarios (which wind up being very similar, of course). Using 200 decimal digits of precision, and Taylor series of degree 100 (so the numerical method was of order 200) we were able to solve the problems accurately in only a few seconds. It is ironic that the "difficult case" resembles so strongly the simple case of a tube of circular cross section, which has a direct and natural solution in terms of Bessel functions.

    When we compare this code with that of [37], we see that that standard code performs very well, in fact. Even just with the default solver (a version of RKF45) all scenarios are rapidly solved. The doubly-exponential growth of the modified Mathieu equation is simply taken in stride by the code. We remark that that code, while more than 20 years old, has undergone steady development since then at the hands of Allan Wittkopf of Maplesoft, who has (without publishing papers on the subject) incorporated many speed and reliability enhancements.

    However, access to the internal interpolants used by the code is quite awkward, and it is not easily possible to differentiate the interpolant to compute a residual to validate the solution it produces. With the present code, this is simple (indeed automatic). Secondly, if very high precision is wanted, the higher order of the present code lowers the cost. Indeed, using this method, the cost of solution is polynomial in the number of bits of accuracy requested [16], while for fixed-order Runge–Kutta methods the cost is exponential in the number of bits of accuracy requested. At modest accuracy, or even at double precision accuracy, this is not a problem, of course.

    The third advantage of the present method is the decent numerical properties of the underlying interpolant. In comparison, the monomial basis used by the internal Maple code can suffer more from rounding errors (at high precision), although we have no doubt that the developers have taken steps to minimize the difficulty.

    Something that might have been a fourth advantage, the ease of combining and integrating blends to compute, for instance, the Green's function, is not much of an advantage after all: Maple's dsolve/numeric interface has several flexible features that let one combine solutions, and integrating the solution of a differential equation is merely a matter of integrating the differential equation for the integral in question.

    Still, this present code offers some potential advantages for other applications, and using this problem as a test case for it has proved to be interesting.

    In the case of a double eigenvalue, the previous formulae need to be amended. For the Mathieu equation, double eigenvalues are isolated, and there are no higher-order eigenvalues, so the treatment is relatively straightforward. The theory has been known for a long time [25], but in practice it seems to have been ignored. We therefore give a detailed treatment below.

    We will use a Puiseux series expansion near the double eigenvalue to deduce the analytical form needed for expansion exactly at the double point. We emphasize that the computations in this section are exact computation of series coefficients, and analytical cancellation of large terms will give us the result that we want.

    Now suppose that the coalescing eigenfunctions are ce0(η,q) and ce2(η,q). For our computations this is the one that mattered the most, when q1.4668i is the Mulholland–Goldstein point; but other purely imaginary eigenvalues also occur for larger frequencies, or larger circumference blood vessels; so we give details of the process.

    If q=1.4687686137851i is the Mulholland–Goldstein point, then we may expand the eigenvalues a0 and a2 in Puiseux series to get a0=aα1qq+O(qq) and a2=a+α1qq+O(qq) in a region close to that point. Here a=2.08869890274969540 is real, and known to many decimal places [4]. Similarly,

    α1=1.659487804320+1.659487804320i (4.1)

    is also known to many decimal places, although as we will see it does not appear in the final formulae for the spectral expansion coefficients at the double eigenvalue.

    If we write a Mathieu series expansion for some function, say v(ξ,η), at a point near to the double point, we find that the coefficients of the terms Ce0(ξ,q)ce0(η,q) and Ce2(ξ,q)ce2(η,q) are large and of opposite sign; indeed they have leading behaviour that is O((qq)1/2). Also, all of Ce0, ce0, Ce2, and ce2 can be written as functions of the fundamental solution wI(z,a,q) as follows:

    ce0(η,q)=wI(η,a0,q)Ce0(ξ,q)=wI(iξ,a0,q)ce2(η,q)=wI(η,a2,q)Ce2(ξ,q)=wI(iξ,a2,q). (4.2)

    We will also need the following two new functions:

    u(η):=D2(wI)(η,a,q) (4.3)
    U(ξ):=D2(wI)(iξ,a,q). (4.4)

    Here D2(f)(x,y,z) means the partial derivative with respect to the second variable, and then evaluated at the point (x,y,z). We will show in a following subsection how these can be computed.

    Now suppose that the solution at the point q has the expansion

    v(ξ,η)=b0Ce0(ξ,q)ce0(η,q)+b2Ce2(ξ,q)ce2(η,q)+ (4.5)

    where the terms not included have eigenvalues that will not coalesce, and therefore the previous treatment using orthogonality will suffice to identify their coefficients b2m for m>1. Putting for brevity x=qq and expanding everything in series in x and neglecting terms of size O(x2) or smaller, we have the following:

    a0=aα1x+O(x2)a2=a+α1x+O(x2)b0=A0x+B0+O(x)b2=A2x+B2+O(x)Ce0(ξ,q)=Ce0(ξ,q)α1U(ξ)x+O(x2)ce0(η,q)=ce0(η,q)α1u(η)x+O(x2)Ce2(ξ,q)=Ce0(ξ,q)+α1U(ξ)x+O(x2)ce2(η,q)=ce0(η,q)+α1u(η)x+O(x2). (4.6)

    Now in our case, the coefficients of v(ξ,η) are determined by integration against the constant function 1 at the wall ξ=ξ0, so that for q near to q we have

    2πη=0ce0(η,q)dηα2πη=0u(η)dηx+O(x2) (4.7)

    on the left-hand side, and, expanding everything out and using the fact that

    2πη=0ce20(η,q)dη=0, (4.8)

    we find

    2α1A0Ce0(ξ0,q)2πη=0ce0(η,q)u(η)dη+L1x+O(x2) (4.9)

    on the right-hand side, with

    L1=α21A0Ce0(ξ0,q)2πη=0u2(η)dη2α1(B0Ce0(ξ0,q)α1A0U(ξ0))2πη=0ce0(η,q)u(η)dη. (4.10)

    Since the squared integral of the generalized eigenfunction u(η)=D2(wI)(η,a,q) is not zero, and since the integral of the product of u(η) with ce0(η,q) is not zero, and since Ce0(ξ0,q) is not zero, we may equate the constant terms and the terms linear in x and solve for A0 and for B0. We get

    A0=2πη=0ce0(η,q)dη2α1Ce0(ξ0,q)2πη=0ce0(η,q)u(η)dη (4.11)

    and

    B0=2πη=0u(η)dη+α1A0(2U(ξ0)2πη=0ce0(η,q)u(η)dη+Ce0(ξ0,q)2πη=0u2(η)dη)2Ce0(ξ0,q)2πη=0ce0(η,q)u(η)dη (4.12)

    Similarly, we get

    A2=A0=2πη=0ce0(η,q)dη2α1Ce0(ξ0,q)2πη=0ce0(η,q)u(η)dη (4.13)

    and

    B2=2πη=0u(η)dηα1A2(2U(ξ0)2πη=0ce0(η,q)u(η)dη+Ce0(ξ0,q)2πη=0u2(η)dη)2Ce0(ξ0,q)2πη=0ce0(η,q)u(η)dη (4.14)

    which resembles the formula for B0.

    Putting these formulae into the expansion for v(ξ,η) we get v(ξ,η)=(B0+B2)Ce0(ξ,q)ce0(η,q)2α1A0(U(ξ)ce0(η,q)+Ce0(ξ,q)u(η))+O(x) and thus as qq the expansion of v(ξ,η) becomes

    v(ξ,η)=b0Ce0(ξ,q)ce0(η,q)+ˆb2(U(ξ)ce0(η,q)+Ce0(ξ,q)u(η))+ (4.15)

    where

    b0=2πη=0u(η)dηCe0(ξ0,q)2πη=0ce0(η,q)u(η)dη+(2πη=0ce0(η,q)dη)(Ce0(ξ0,q)2πη=0u2(η)dη+U(ξ0)2πη=0ce0(η,q)u(η)dη)(Ce0(ξ0,q)2πη=0ce0(η,q)u(η)dη)2 (4.16)
    ˆb2=2πη=0ce0(η,q)dηCe0(ξ0,q)2πη=0ce0(η,q)u(η)dη, (4.17)

    and the omitted terms can all be calculated by orthonormality as before.

    Remark 4.1. This result can be derived a different way, by differentiating the original formula with respect to a. With that method, the appearance of U(ξ)ce0(η)+Ce0(ξ)u(η), being the derivative of Ce0(ξ)ce0(η), seems natural. Then one can use the orthogonality of ce0(η) with all ce2m(η) (including itself) to compute ˆb2, and then integrate against u(η) and solve the resulting equation for b0 using the known ˆb2. This leads to the same result, but we feel that the detailed derivation above is more convincing, and explains what happens to the expansion coefficients as qq.

    All that remains is the computation of u(η) and U(ξ). To do this, we compute the Fréchet derivatives of the Mathieu equation and the modified Mathieu equation:

    d2udη2+(a2qcos2η)u+y=0 (4.18)
    d2Udξ2(a2qcosh2ξ)Uy=0. (4.19)

    In the first equation, replace y by ce0(η,q) and solve (we use the Green's function for the Mathieu equation to do so, because algebraic operations and integration are accurate and efficient with blendstrings) and similarly in the second equation replace y by Ce0(ξ,q) and solve. As for initial or boundary conditions, we need to take u(0)=u(2π)=0 to ensure periodicity, and we need to take U(0)=U(0)=0 to ensure symmetry at the line ξ=0.

    This analysis is implicit in the treatment in [25], but does not seem to be widely pursued, and so we have written it down in some detail here.

    Finally, we must amend the formulas for flow rate and oscillatory wall shear stress. Equation (2.54) becomes

    qϕ,e(t)=8q0,eiλe(11iπλetanh2ξ0(b0Ce0(ξ0,q)2πη=0ce0(η,q)+ˆb2(U(ξ0)2πη=0ce0(η,q)dη+Ce0(ξ0,q)2πη=0u(η)dη)+m2b22mI2mCe2m(ξ0,q)Ce2m(ξ0,q)))eiωt. (4.20)

    The integrals appearing in the formula above have already been calculated, in order to find the b2m, but the relationships used to simplify to get b22mI2m as in the rest of the sum no longer obtain because I0=0.

    Equation (2.59) becomes

    τϕ,e(η,t)=2ˆτ0,eiδ0λe(ˆb2(U(ξ0)ce0(η,q)+Ce0(ξ0,q)u(η))+m1b2mCe2m(ξ0,q)ce2m(η,q))eiωt. (4.21)

    Because the solution to the original model equations is continuous (indeed analytic) in the parameters involved in q=iΛe/4, the underlying solution changes continuously as q passes through a value where a double eigenvalue of the Mathieu equation occurs. Therefore, it is only a discontinuity in the representation of the solution, not the solution itself. This means that sampling q "near enough" to the double point would give solutions that are "near enough" to the solution at that point.

    The only difficulty, and this is rather mild, is that the expansion coefficients in Mathieu functions become large and of opposite sign, which might incur some visible rounding error owing to cancellation. Because the size of the coefficients is only O(qq)1/2 this is not typically very severe.

    Nonetheless we feel that it is worthwhile to be able to give the precise solution exactly at a double point for comparison to simple solutions nearby, to be assured that the solutions shown are representative of the model.

    In the results to follow (Figures 913) we consider a change in the cross section of a tube from circular to elliptic, under both passive and active scenarios, and examine the effects of this on the properties of oscillatory flow in a tube of elliptic cross section under the same oscillatory pressure gradient as that in a tube of circular cross section. The effects of physiological interest are those on flow rate and on the distribution of shear stress around the circumference of the tube. The main focus of our study is therefore on these two properties as well on the form of the passive scenario of velocity profiles in a tube of elliptic cross section.

    Figure 9.  Left: A blood vessel of circular cross section with circumference 1.0 cm and radius a=1.0/(2π) cm (blue dashed line) is deformed by imposed vertical forces to ellipses with semiminor axis β=fea with fe=0.6, in two different ways: under the active scenario, neuromuscular control relaxes the vessel wall so it stretches in order to maintain the initial cross-sectional area (red curves), and under the passive scenario, the vessel wall stays at 1.0cm in circumference (black curves), thereby losing some cross-sectional area. Right: The resulting peak shear stress over one cycle of oscillation from equation (2.59), scaled by the steady wall stress ˆτ0,c from equation (2.17) in the case α=β=a, the radius of the circle. Both scenarios are graphed at their peaks in time, together with the wall stress from equation (2.60) (blue dashed line) from the original circle. The wider elliptical vessel has the smallest minimum oscillatory wall shear stress, although the maximum is greater than that of the original circle.
    Figure 10.  One quarter of a cycle of an oscillation starting at a maximum, with each curve showing (half) the velocity profile from equation (2.42) at successive instants. The top curve is at time t=0. As the time t is sampled over the quarter cycle, the curves go down. In the next quarter cycle (not shown) they would go further down to the minimum; then they would come up over the next half-cycle (not shown, because the overlapping lines would be harder to interpret) to the beginning again. The particulars of this picture are that the original circumference was c=2.0 cm, the semiminor axis of the ellipse is a fraction fe=0.95 of the original circular radius, the scenario is active so that the semimajor axis is 1/fe=1.053 of the original circular radius, and the frequency is ω=11.2847265810330Hz (which is slightly higher than for most of our simulations). This is, as it happens, at a double eigenvalue, which is why so many decimal places of the frequency are printed. Fractions fe very close to 1 are harder for the Mathieu expansion to model, as indicated in the text, but this fraction, fe=0.95, is routine, and slight deformations from circular cross section are more likely to occur, so this is a reasonably realistic case.
    Figure 11.  One quarter of a cycle of an oscillation, starting at time t=0 half-way between the minimum and the maximum. Each curve shows (half) the velocity profile from Eq (2.42) at successive instants. Successive curves show the profile going down as the sampled times increase. In the next quarter cycle (not shown) they would go up; then they would go up further to the maximum for the cycle, then come back down over a final quarter cycle to the beginning again (not shown because the overlapping curves would be harder to interpret). The particulars are that the original circular circumference was c=2.0 cm, the semiminor axis is a fraction fe=0.6 of the original radius, the scenario is an active one and so the semimajor axis is 1/fe=1.667 of the original radius, and the frequency ω=0.959305048332072Hz which is well within the frequency range of most of our simulations. This is again, as it happens, at a double eigenvalue, which is why so many decimal places of the frequency are printed. The double eigenvalue solution is very similar to nearby simple eigenvalue solutions, as discussed in the text. The fraction fe=0.6 seems reasonable for the active scenario, because according to Figure ??? the circumference would need to stretch by perhaps 25%.
    Figure 12.  The real part of the oscillatory wall shear stress τe/ˆτ0,e throughout a cycle from equation (2.59) compared in the passive scenario (left) to the active (right). The initial circumference was c=2.0cm. The fraction of the original circular radius is fe=0.3 in both cases, and in both cases ω=1 Hz. The main difference between the two figures, which is slight, occurs at the poles of the major axis, where the shear stress is least. There the wall shear stress is smaller in the active scenario than it is in the passive scenario.
    Figure 13.  Minimum oscillatory wall shear stress from equation (2.59) compared in the passive scenario (left) to the active (right). The initial circumference was c=1.0 cm. The same range of fractions fe are graphed for each scenario, namely 0.3, 0.4, 0.6, 0.8, 0.9, and 0.95 over the same range of frequencies ω. In this figure we see that the different scenarios do not produce greatly different minimum wall shear stresses, with the greatest differences appearing for the smallest fractions fe; and moreover that the different scenarios show a similar dependence on the imposed frequency ω, with a general decrease in shear stress for higher frequency. This dependence is weaker with the more compressed cross sections (smaller fractions fe). We used the simple symbol f instead of fe in the legend to improve readability.

    As noted, the properties of oscillatory flow in a tube of elliptic cross section depend on the nondimensional parameter Λe (Eq 2.28) which involves the frequency of oscillation, ω, as well as the focal distance, d. Thus the effects of tube dimension on oscillatory flow in the tube of elliptic cross section are different at different frequencies and, similarly, the effects of frequency on oscillatory flow in the tube of elliptic cross section are different at different tube dimensions. As a consequence, the effects of tube dimension and of frequency cannot be scaled out and, in the results to follow we examine three specific values of the circumference and thus radii, deformed by forcing them to different fractions fe of their original radii, and several specific values of frequency, as shown in the figures.

    All the results to follow are based on the real part of the oscillatory pressure gradient (Eq (2.23)). We note that the absolute value of the various complex quantities must appear at some point in the flow, possibly with a different phase lag for different locations in the vessel. In our animations (not given here) the differences in phase lags were never very significant.

    The fluid density and viscosity in all the results were taken as 1.0 g/cm3 and 0.04 g/(cms), respectively. We note that we only examine tubes equivalent to those of radius of 0.5 cm, 1.0 cm, and 2.0 cm in both active and passive scenarios for maximum flow rate and maximum wall shear stress.

    The primary factor in the transition of pulsatile blood flow in a vessel of circular cross section to one in a vessel of an elliptic cross section is the loss of radial symmetry of the circular cross section. While from a geometrical perspective this loss of symmetry appears to occur fairly smoothly, from both a mathematical and a hemodynamical perspective it represents a significant change. Geometrically, the change from a circular to an elliptic cross section, however small, introduces "poles" in the cross section, places where the curvature is maximum—at the ends of the major axis—and where the curvature is minimum, at the ends of the minor axis.

    Further, the most convenient coordinate system, namely confocal elliptical coordinates, has singular behaviour in the limit as the focus distance d0. This in turn induces a change in the governing equations of pulsatile flow from Bessel equations to Mathieu equations. Hemodynamically, the change causes a redistribution of shear stress on the vessel boundary, from a uniform distribution in the case of circular cross section to a polarized distribution in the case of elliptic cross section, with maximum shear occurring at the two ends of the minor axis of the ellipse and minimum shear at the two ends of the major axis.

    From the perspective of blood flow regulation, which our study was aimed at, the transition from flow in a vessel of circular cross section to one in an elliptic cross section represents a departure from well known physiological rules of blood flow regulation to a somewhat uncharted territory. A simple change in the diameter of a vessel is well known as the physiological (neurovascular) mechanism used to change the cross sectional area of a blood vessel in order to affect a required change in flow rate.

    Our study was aimed at the question of how this well established rule of blood flow regulation might be altered in the case of an elliptic cross section. Our results indicate that if the regulatory system does not respond to the change from circular to elliptic cross section, which we have dubbed as a "passive scenario", the change from circular to elliptic cross section will occur with no change in the length of circumference of the changing cross section. As a consequence, the cross sectional area available to the flow will then be reduced under this scenario.

    If, on the other hand, the regulatory system intervenes in an attempt to maintain the cross sectional area available to the flow, as it does in the case of a circular cross section, the transition from circular to elliptic cross section will occur while keeping the cross sectional area constant and hence, necessarily, by increasing the length of circumference of the changing cross section. We have dubbed this as an "active scenario".

    This makes a difference even in steady flow in tubes of elliptic cross section, as already seen in Figures 57. The flow quantities tend to be higher under the active scenario. This persists for pulsatile flow, as we will see; but we measure the pulsatile flow quantities relative to their steady counterparts, and so the increase may be hidden. This must be kept in mind.

    While what has been said so far applies to steady flow, the effects of ellipticity in oscillatory flow are further complicated by the acceleration and deceleration of the fluid within the oscillatory cycle. The effects of acceleration and deceleration depend on the volume of fluid being accelerated and decelerated, which in turn depends on the cross sectional area of the tube in which the flow is taking place. We recall that as the circular cross section of a tube becomes elliptic in the passive scenario, the cross sectional area of the elliptic cross section becomes smaller than that of the circular cross section, and therefore a smaller volume of fluid will be accelerated and decelerated in the tube of elliptic cross section than in the circular one. It follows that the maximum flow rate reached at the peak of each oscillatory cycle might actually be higher in the tube of elliptic cross section. This was shown, however, not to be the case for steady flow in Figure 7. The reason for this is that the acceleration and deceleration peaks depend not only on the volume of fluid being accelerated but also on the opposition to that acceleration by the level and distribution of shear stress on the boundary. However, even for steady flow, shear stress at the tube wall can be higher in the tube of elliptic cross section, and is in both scenarios for fe0.6, as shown in Figure 6. For pulsatile flow, Figure 14 shows that this remains true for pulsatile flow.

    Figure 14.  Maximum oscillatory wall shear stress from equation (2.59) compared in the passive scenario (left) to the active (right). The initial circumference was c=1.0 cm. The same range of fractions fe are graphed for each scenario, namely 0.3, 0.4, 0.6, 0.8, 0.9, and 0.95 over the same range of frequencies ω.In this figure we see that the different scenarios produce very similar maximum wall shear stresses, and moreover that the different scenarios show a similar dependence on the imposed frequency ω, with a general decrease in shear stress for higher frequency. This dependence is weaker with the more compressed cross sections (smaller fractions fe). We used the simple symbol f instead of fe in the legend to improve readability.

    The pulsatile flow rates shown in Figures 1517 do not become higher for tubes of elliptic cross section than in circular cross section, even in the active scenario. They show moderate dependence on imposed frequency ω with a general downward trend with increasing frequency. They also show a decrease in maximum flow for tubes with smaller fractions fe of the original radius. They do not show a great effect of the two different scenarios, active versus passive. Indeed, those figures show that there is very little effect on maximum flow rate between the scenarios. The maximum velocity profiles show similar behaviour, and are not plotted here as being redundant.

    Figure 15.  Oscillatory flow rate from equation (2.54) compared in the passive scenario (left) to the active (right). The initial circumference was c=0.5cm. The same range of fractions fe are graphed for each scenario, namely 0.4, 0.6, 0.8, 0.9, and 0.95 over the same range of frequencies ω.
    Figure 16.  Oscillatory flow rate from equation (2.54) compared in the passive scenario (left) to the active (right). The initial circumference was c=1.0 cm. The same range of fractions fe are graphed for each scenario, namely 0.4, 0.6, 0.8, 0.9, and 0.95 over the same range of frequencies ω. These graphs do not at first look greatly different from those at c=0.5cm in Figure 15, but note the vertical scale: over this range of frequencies ω the flow rate drops by about 30% whereas with c=0.5 the flow rate dropped only by about 3%.
    Figure 17.  Oscillatory flow rate from Eq (2.54) compared in the passive scenario (left) to the active (right). The initial circumference was c=2.0 cm. The same range of fractions fe are graphed for each scenario, namely 0.3, 0.4, 0.6, 0.8, 0.9, and 0.95 over the same range of frequencies ω. These graphs show a much greater dependence on ω than those of either c=1.0cm or c=0.5cm in Figures 15 and 16.

    Tubes of elliptic cross section with smaller fe show greater dependence of maximum pulsatile flow rate (and similarly velocity, not shown) on the imposed frequency. Conversely, tubes of smaller fe show weaker dependence of shear stress on the imposed frequency ω.

    As a final consideration, the dependence of oscillatory properties on frequency ω seen in the figures can be interpreted as the way the first few harmonics of a composite pressure wave would be individually affected under the two scenarios being considered. We do not otherwise pursue here the idea of a composite pressure wave.

    The problem of pulsatile flow in tubes of elliptic cross sections is important from a physiological as well as mathematical perspective, and the aim of our study was to examine this problem from both of these perspectives, using a tube of elliptic cross section as a model of a deformed blood vessel. While this is clearly a simplified model of the many ways in which a blood vessel may be deformed, it allowed us to explore a full range of distortions of a tube of circular cross section, from being fully open to almost closed.

    More important than the final form into which a vessel is deformed are the constraints and scenarios under which the transformation from circular to elliptic cross section takes place. The two scenarios which we have considered highlight the mathematical and physiological aspects of the problem and provide useful information on the way the neurovascular control system may respond to the deformation of a blood vessel in the physiological setting. In particular, the ability of the control system to maintain a constant cross sectional area under the active scenario is clearly limited to only small or moderate departures from the circular cross section. When the departure from circular cross section is large (high ellipticity, low fraction fe), a prohibitively large increase in the circumference of the wall would be required to maintain the cross sectional area available for the flow, as illustrated in Figure 4.

    We have extended both the scope and the range of results currently available for this problem by using new methodology to overcome difficulties encountered in the solution of the governing Mathieu equations and in the numerical evaluation of Mathieu functions in the past. Specifically, we used a careful spectral method, including explicit solution in the case of double eigenvalues, for the solution of the governing equations. We used extended precision where necessary to overcome issues of ill-conditioning, for fe very close to 1 which is paradoxically the difficult case. We believe that this novel approach offers a useful new tool in further study of pulsatile blood flow under various pathological conditions.

    Limitations of this study fall into two separate categories.

    In the physiological/clinical category the most obvious limitation is that under pathological conditions a blood vessel may be compressed in a variety of different ways which do not necessarily lead to elliptical cross sections. Under these circumstances, depending on the resulting geometry, a study must clearly be tailored specifically to that geometry. In cases where the geometry is not sufficiently simple to permit solution, or indeed application, of the same governing equations, the present study may be seen as a mathematically tractable approximation.

    In the mathematical category the present study does not deal with computational difficulties arising in the special cases of nearly circular or highly elliptic cross sections. The computational issues and software requirements involved in these cases are fully discussed in [4]. The main focus in the present paper is on hemodynamic consequences.

    The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

    This work was supported by NSERC under grants numbered RGPIN-2020-06438 (RMC) and RGPIN-2019-04749 (MZ); and, while RMC was visiting the Isaac Newton Institute during the programme Complex Analysis: Tools, techniques, and applications, by EPSRC Grant # EP/R014604/1.

    All authors declare no conflicts of interest in this paper.



    [1] W. Chen, T. Zhou, Existence of solutions for p-Laplacian parabolic Kirchhoff equation, Appl. Math. Lett., 122 (2021), 107527. http://dx.doi.org/10.1016/j.aml.2021.107527 doi: 10.1016/j.aml.2021.107527
    [2] I. Lasiecka, J. H. Rodrigues, Weak and strong semigroups in structural acoustic Kirchhoff-Boussinesq interactions with boundary feedback, J. Differ. Equations, 298 (2021), 387–429. http://dx.doi.org/10.1016/j.jde.2021.07.009 doi: 10.1016/j.jde.2021.07.009
    [3] C. Vetro, Variable exponent p(x)-Kirchhoff type problem with convection, J. Math. Anal. Appl., 506 (2022), 125721. http://dx.doi.org/10.1016/j.jmaa.2021.125721 doi: 10.1016/j.jmaa.2021.125721
    [4] N. D. Phuong, N. H. Tuan, Z. Hammouch, R. Sakthivel, On a pseudo-parabolic equations with a non-local term of the Kirchhoff type with random Gaussian white noise, Chaos Soliton. Fract., 145 (2021), 110771. http://dx.doi.org/10.1016/j.chaos.2021.110771 doi: 10.1016/j.chaos.2021.110771
    [5] M. Xiang, D. Hu, D. Yang, Least energy solutions for fractional Kirchhoff problems with logarithmic nonlinearity, Nonlinear Anal., 198 (2020), 111899. http://dx.doi.org/10.1016/j.na.2020.111899 doi: 10.1016/j.na.2020.111899
    [6] M. Xiang, D. Yang, Nonlocal Kirchhoff problems: Extinction and non-extinction of solutions, J. Math. Anal. Appl., 477 (2019), 133–152. http://dx.doi.org/10.1016/j.jmaa.2019.04.020 doi: 10.1016/j.jmaa.2019.04.020
    [7] Y. Han, Q. Li, Threshold results for the existence of global and blow-up solutions to Kirchhoff equations with arbitrary initial energy, Comput. Math. Appl., 75 (2018), 3283–3297. http://dx.doi.org/10.48550/arXiv.1703.09094 doi: 10.48550/arXiv.1703.09094
    [8] B. Guo, H. Zhou, Output feedback stabilization for multi-dimensional Kirchhoff plate with general corrupted boundary observation, Eur. J. Control, 28 (2016), 38–48. http://dx.doi.org/10.1016/j.ejcon.2015.12.004 doi: 10.1016/j.ejcon.2015.12.004
    [9] I. Lasiecka, M. Pokojovy, X. Wand, Global existence and exponential stability for a nonlinear thermoelastic Kirchhoff-Love plate, Nonlinear Anal.-Real, 38 (2017), 184–221. http://dx.doi.org/10.1016/J.NONRWA.2017.04.001 doi: 10.1016/J.NONRWA.2017.04.001
    [10] M. Ghisi, M. Gobbino, Optimal decay-error estimates for the hyperbolic-parabolic singular perturbation of a degenerate nonlinear equation, J. Differ. Equations, 254 (2013), 911–932. http://dx.doi.org/10.1016/j.jde.2012.10.005 doi: 10.1016/j.jde.2012.10.005
    [11] Z. Yang, Longtime behavior of the Kirchhoff type equation with strong damping on RN, J. Differ. Equations, 242 (2007), 269–286. http://dx.doi.org/10.1016/j.jde.2007.08.004 doi: 10.1016/j.jde.2007.08.004
    [12] I. Lasiecka, M. Pokojovy, X. Wan, Long-time behavior of quasilinear thermoelastic Kirchhoff-Love plates with second sound, Nonlinear Anal., 186 (2019), 219–258. http://dx.doi.org/10.48550/arXiv.1811.01138 doi: 10.48550/arXiv.1811.01138
    [13] T. Boudjeriou, M. K. Hamdani, M. Bayrami-Aminloue, Global existence, blow-up and asymptotic behavior of solutions for a class of p(x)-Choquard diffusion equations in RN, J. Math. Anal. Appl., 506 (2022), 125720. http://dx.doi.org/10.1016/j.jmaa.2021.125720 doi: 10.1016/j.jmaa.2021.125720
    [14] Y. Sun, H. Wang, Study of weak solutions for a class of degenerate parabolic variational inequalities with variable exponent, Symmetry, 14 (2022), 1255. http://dx.doi.org/10.3390/sym14061255 doi: 10.3390/sym14061255
    [15] D. Adak, G. Manzini, S. Natarajan, Virtual element approximation of two-dimensional parabolic variational inequalities, Comput. Math. Appl., 116 (2022), 48–70. http://dx.doi.org/10.1016/j.camwa.2021.09.007 doi: 10.1016/j.camwa.2021.09.007
    [16] J. Dabaghi, V. Martin, M. Vohralík, A posteriori estimates distinguishing the error components and adaptive stopping criteria for numerical approximations of parabolic variational inequalities, Comput. Method. Appl. M., 367 (2020), 113105. http://dx.doi.org/10.1016/j.cma.2020.113105 doi: 10.1016/j.cma.2020.113105
    [17] J. Li, C. Bi, Study of weak solutions of variational inequality systems with degenerate parabolic operators and quasilinear terms arising Americian option pricing problems, AIMS Math., 7 (2022), 19758–19769. http://dx.doi.org/10.3934/math.20221083 doi: 10.3934/math.20221083
    [18] X. Chen, F. Yi, Parabolic variational inequality with parameter and gradient constraints, J. Math. Anal. Appl., 385 (2012), 928–946. http://dx.doi.org/10.1016/j.jmaa.2011.07.025 doi: 10.1016/j.jmaa.2011.07.025
    [19] X. Chen, F. Yi, L. Wang, American lookback option with fixed strike price 2-D parabolic variational inequality, J. Differ. Equations, 251 (2011), 3063–3089. http://dx.doi.org/10.1016/j.jde.2011.07.027 doi: 10.1016/j.jde.2011.07.027
  • This article has been cited by:

    1. Robert M. Corless, An Hermite–Obreshkov method for 2nd-order linear initial-value problems for ODE, 2024, 96, 1017-1398, 1109, 10.1007/s11075-023-01738-z
  • Reader Comments
  • © 2023 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(1366) PDF downloads(55) Cited by(6)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog