The Mathieu Differential Equation and Generalizations to Infinite Fractafolds

TL;DR

This paper analyzes the Mathieu differential equation's solutions, explores their stability, and extends the framework to fractal structures like the Sierpinski gasket, combining theoretical and numerical approaches.

Contribution

It introduces a method to generalize the Mathieu equation to fractal domains and studies the stability and properties of solutions in this new setting.

Findings

Fourier series approach effectively approximates solutions

Stability results for classical Mathieu equation confirmed

Extension to fractal domain demonstrates similar solution behavior

Abstract

One of the well-studied equations in the theory of ODEs is the Mathieu differential equation. A common approach for obtaining solutions is to seek solutions via Fourier series by converting the equation into an infinite system of linear equations for the Fourier coefficients. We study the asymptotic behavior of these Fourier coefficients and discuss the ways in which to numerically approximate solutions. We present both theoretical and numerical results pertaining to the stability of the Mathieu differential equation and the properties of solutions. Further, based on the idea of using Fourier series, we provide a method in which the Mathieu differential equation can be generalized to be defined on the infinite Sierpinski gasket. We discuss the stability of solutions to this fractal differential equation and describe further results concerning properties and behavior of these solutions.

Figures

Tables

Table 1. Table 1. The probability P i subscript 𝑃 𝑖 P_{i} ’s.

$i$	$P_i$
1	0.625056436387445
2	0.784428425813594
3	0.845139663995868
4	0.878143787704672
5	0.899154589232086
6	0.913800056779179
7	0.924632580597566
8	0.932988858616457
9	0.939640860147421
10	0.945067300763657
11	0.949581603650156
12	0.953397960289362
13	0.956667925443240
14	0.959502140424685
15	0.961982968671620
16	0.964174006346465
17	0.966133547589692
18	0.967239011961280

Equations

$$u(t)=\sum _{j=0}^{\infty }{a}_j\cos \left(\frac{j}{2}t\right)+\sum _{j=1}^{\infty }{b}_j\sin \left(\frac{j}{2}t\right).$$

$$\text{cosine coefficients}\{\begin{array}{c}\delta a_0+\frac{\epsilon }{2}{a}_2=0,\\ \left(\delta -{\left(\frac{1}{2}\right)}^2\right)a_1+\frac{\epsilon }{2}\left(a_3+a_1\right)=0,\\ (\delta -1)a_2+\epsilon a_0+\frac{\epsilon }{2}{a}_4=0,\\ \left(\delta -{\left(\frac{j}{2}\right)}^2\right)a_j+\frac{\epsilon }{2}\left(a_{j-2}+a_{j+2}\right)=0,(j\ge 3),\\ \end{array}$$

$$\text{sine coefficients}\{\begin{array}{c}\left(\delta -{\left(\frac{1}{2}\right)}^2\right)b_1+\frac{\epsilon }{2}(b_3-b_1)=0,\\ (\delta -1)b_2+\frac{\epsilon }{2}{b}_4=0,\\ \left(\delta -{\left(\frac{j}{2}\right)}^2\right)b_j+\frac{\epsilon }{2}(b_{j-2}+b_{j+2})=0,(j\ge 3).\end{array}$$

$$\underset{A}{\underbrace{\left(\begin{array}{cccccc}\delta & \frac{\epsilon }{2}& & & & \\ \epsilon & \delta -1^2& \frac{\epsilon }{2}& & & \\ & \frac{\epsilon }{2}& \delta -2^2& \frac{\epsilon }{2}& & \\ & & \frac{\epsilon }{2}& \delta -3^2& \frac{\epsilon }{2}& \\ & & & \ddots & \ddots & \ddots \end{array}\right)}}\underset{{\text{a}}_{even}}{\underbrace{\left(\begin{array}{c}{a}_0\\ a_2\\ a_4\\ a_6\\ ⋮\end{array}\right)}}=\underset{\text{0}}{\underbrace{\left(\begin{array}{c}0\\ 0\\ 0\\ 0\\ ⋮\end{array}\right)}}.$$

$$\underset{B}{\underbrace{\left(\begin{array}{ccccc}\delta -1^2& \frac{\epsilon }{2}& & & \\ \frac{\epsilon }{2}& \delta -2^2& \frac{\epsilon }{2}& & \\ & \frac{\epsilon }{2}& \delta -3^2& \frac{\epsilon }{2}& \\ & & \ddots & \ddots & \ddots \end{array}\right)}}\underset{{\text{b}}_{even}}{\underbrace{\left(\begin{array}{c}{b}_2\\ b_4\\ b_6\\ ⋮\end{array}\right)}}=\left(\begin{array}{c}0\\ 0\\ 0\\ ⋮\end{array}\right).$$

$$\underset{C}{\underbrace{\left(\begin{array}{cccccc}\delta -\frac{1}{4}+\frac{\epsilon }{2}& \frac{\epsilon }{2}& & & & \\ \epsilon & \delta -\frac{9}{4}& \frac{\epsilon }{2}& & & \\ & \frac{\epsilon }{2}& \delta -\frac{25}{4}& \frac{\epsilon }{2}& & \\ & & \frac{\epsilon }{2}& \delta -\frac{49}{4}& \frac{\epsilon }{2}& \\ & & & \ddots & \ddots & \ddots \end{array}\right)}}\underset{{\text{a}}_{odd}}{\underbrace{\left(\begin{array}{c}{a}_1\\ a_3\\ a_5\\ a_7\\ ⋮\end{array}\right)}}=\left(\begin{array}{c}0\\ 0\\ 0\\ 0\\ ⋮\end{array}\right).$$

$$\underset{D}{\underbrace{\left(\begin{array}{cccccc}\delta -\frac{1}{4}-\frac{\epsilon }{2}& \frac{\epsilon }{2}& & & & \\ \epsilon & \delta -\frac{9}{4}& \frac{\epsilon }{2}& & & \\ & \frac{\epsilon }{2}& \delta -\frac{25}{4}& \frac{\epsilon }{2}& & \\ & & \frac{\epsilon }{2}& \delta -\frac{49}{4}& \frac{\epsilon }{2}& \\ & & & \ddots & \ddots & \ddots \end{array}\right)}}\underset{{\text{b}}_{odd}}{\underbrace{\left(\begin{array}{c}{b}_1\\ b_3\\ b_5\\ b_7\\ ⋮\end{array}\right)}}=\left(\begin{array}{c}0\\ 0\\ 0\\ 0\\ ⋮\end{array}\right).$$

$$u(t):=\sum _{j=0,2,4,6,...}{a}_j\cos \left(\frac{j}{2}t\right)+\sum _{j=2,4,6,...}{b}_j\sin \left(\frac{j}{2}t\right)$$

$$u(t):=\sum _{j=1,3,5,...}{a}_j\cos \left(\frac{j}{2}t\right)+\sum _{j=1,3,5,...}{b}_j\sin \left(\frac{j}{2}t\right)$$

$$A_m=\left(\begin{array}{cccccc}\delta & \frac{\epsilon }{2}& & & & \\ \epsilon & \delta -1^2& \frac{\epsilon }{2}& & & \\ & \frac{\epsilon }{2}& \delta -2^2& \frac{\epsilon }{2}& & \\ & & \frac{\epsilon }{2}& \delta -3^2& \frac{\epsilon }{2}& \\ & & & \ddots & \ddots & \ddots \\ & & & & \frac{\epsilon }{2}& \delta -(m-1{)}^2\end{array}\right)$$

$$B_m=\left(\begin{array}{cccccc}\delta -1^2& \frac{\epsilon }{2}& & & & \\ \frac{\epsilon }{2}& \delta -2^2& \frac{\epsilon }{2}& & & \\ & \frac{\epsilon }{2}& \delta -3^2& \frac{\epsilon }{2}& & \\ & & \frac{\epsilon }{2}& \delta -4^2& \frac{\epsilon }{2}& \\ & & & \ddots & \ddots & \ddots \\ & & & & \frac{\epsilon }{2}& \delta -m^2\end{array}\right)$$

$$\left(\begin{array}{ccccc}\delta -{\lambda }_1-\gamma \epsilon & {\alpha }_1\epsilon & & & \\ {\beta }_1\epsilon & \delta -{\lambda }_2& {\alpha }_2\epsilon & & \\ & {\beta }_2\epsilon & \delta -{\lambda }_3& {\alpha }_3\epsilon & \\ & & {\beta }_3\epsilon & \delta -{\lambda }_4& \ddots \\ & & & \ddots & \ddots \end{array}\right),$$

$$\left(\begin{array}{ccccc}\delta -{\lambda }_1-\gamma \epsilon & {\alpha }_1\epsilon & & & \\ {\beta }_1\epsilon & \delta -{\lambda }_2& {\alpha }_2\epsilon & & \\ & {\beta }_2\epsilon & \delta -{\lambda }_3& {\alpha }_3\epsilon & \\ & & {\beta }_3\epsilon & \delta -{\lambda }_4& \ddots \\ & & & \ddots & \ddots \end{array}\right)\left(\begin{array}{c}{c}_1\\ c_2\\ c_3\\ c_4\\ ⋮\end{array}\right)=\left(\begin{array}{c}0\\ 0\\ 0\\ 0\\ ⋮\end{array}\right),$$

$$\{\begin{array}{l}(\delta -{\lambda }_1-\gamma \epsilon )c_1+{\alpha }_1\epsilon c_2=0\text{,}\\ {\beta }_{j-1}\epsilon c_{j-1}+(\delta -{\lambda }_j)c_j+{\alpha }_j\epsilon c_{j+1}=0\text{, for }j\ge 2,\end{array}$$

$$\underset{j\to \infty }{\lim }\frac{x_j}{x_{j-1}}=t_1,\underset{j\to \infty }{\lim }\frac{y_j}{y_{j-1}}=t_2.$$

$$\underset{j\to \infty }{\mathrm{limsup}}|c_j{|}^{\frac{1}{j}}=|t_1|$$

$$f_j{c}_{j+1}+d_j{c}_j+g_{j-1}{c}_{j-1}=0(j\ge 1),$$

$$\frac{x_j}{x_{j-1}}\sim -\frac{g_{j-1}}{d_j},\frac{y_j}{y_{j-1}}\sim -\frac{d_{j-1}}{f_{j-1}}.$$

$${\hat{c}}_{j+1}+{\hat{c}}_j+\frac{f_{j-1}{g}_{j-1}}{d_{j-1}{d}_j}{\hat{c}}_{j-1}=0.$$

$$\underset{j\to \infty }{\lim }\frac{{\hat{x}}_j}{{\hat{x}}_{j-1}}=0,\underset{j\to \infty }{\lim }\frac{{\hat{y}}_j}{{\hat{y}}_{j-1}}=-1.$$

$$H_{j+1}^{(i)}{H}_j^{(i)}+H_j^{(i)}+\frac{f_{j-1}{g}_{j-1}}{d_{j-1}{d}_j}=0,i=1,2.$$

$$\{\begin{array}{l}{H}_{j+1}^{(1)}{H}_j^{(1)}{H}_j^{(2)}+H_j^{(1)}{H}_j^{(2)}+\frac{f_{j-1}{g}_{j-1}}{d_{j-1}{d}_j}{H}_j^{(2)}=0,\\ H_{j+1}^{(2)}{H}_j^{(2)}{H}_j^{(1)}+H_j^{(1)}{H}_j^{(2)}+\frac{f_{j-1}{g}_{j-1}}{d_{j-1}{d}_j}{H}_j^{(1)}=0.\end{array}$$

$$H_j^{(1)}{H}_j^{(2)}\left(H_{j+1}^{(1)}-H_{j+1}^{(2)}\right)=\frac{f_{j-1}{g}_{j-1}}{d_{j-1}{d}_j}\left(H_j^{(1)}-H_j^{(2)}\right).$$

$$H_j^{(1)}\sim 0\text{ and }{H}_j^{(2)}\sim -1,$$

$$H_j^{(1)}=\frac{f_{j-1}{g}_{j-1}}{d_{j-1}{d}_j}\frac{1}{H_j^{(2)}}\frac{H_j^{(1)}-H_j^{(2)}}{H_{j+1}^{(1)}-H_{j+1}^{(2)}}\sim -\frac{f_{j-1}{g}_{j-1}}{d_{j-1}{d}_j}.$$

$$\frac{x_j}{x_{j-1}}=\frac{\left({\prod }_{i=1}^{j-1}\frac{d_i}{f_i}\right){\hat{x}}_j}{\left({\prod }_{i=1}^{j-2}\frac{d_i}{f_i}\right){\hat{x}}_{j-1}}=H_j^{(1)}\frac{d_{j-1}}{f_{j-1}}\sim -\frac{g_{j-1}}{d_j},$$

$$\frac{y_j}{y_{j-1}}=\frac{\left({\prod }_{i=1}^{j-1}\frac{d_i}{f_i}\right){\hat{y}}_j}{\left({\prod }_{i=1}^{j-2}\frac{d_i}{f_i}\right){\hat{y}}_{j-1}}=H_j^{(2)}\frac{d_{j-1}}{f_{j-1}}\sim -\frac{d_{j-1}}{f_{j-1}}.\text{hfill}\square$$

$$\frac{c_k}{c_{k-1}}{\left(\frac{x_k}{x_{k-1}}\right)}^{-1}-1=\frac{r_N}{r_{k-1}}\cdot \frac{1-r_k/r_{k-1}}{1-r_N/r_{k-1}},$$

$$\left(\begin{array}{ccccc}\delta -{\lambda }_1-\gamma \epsilon & {\alpha }_1\epsilon & & & \\ {\beta }_1\epsilon & \delta -{\lambda }_2& {\alpha }_2\epsilon & & \\ & {\beta }_2\epsilon & \delta -{\lambda }_3& {\alpha }_3\epsilon & \\ & & {\beta }_3\epsilon & \ddots & \ddots \\ & & & \ddots & \end{array}\right)$$

$$=\delta I-T(\epsilon ),$$

$$T=\left(\begin{array}{ccccc}{d}_1& f_1& & & \\ f_1& d_2& f_2& & \\ & f_2& d_3& f_3& \\ & & f_3& \ddots & \ddots \\ & & & \ddots & \end{array}\right),$$

Full text

1. 1 Introduction

2. 2 Definitions and Methods

3. 2.1 Fourier Expansion and Matrix Form

4. 2.2 Truncation Method

5. 3 Three-Term Recurrence Relations and Truncations of the Infinite Matrices

6. 3.1 Asymptotic Behavior of $\{c_{j}\}_{j=1}^{\infty}$

7. 3.2 The Truncation Method

8. 4 Observations and Analysis for the Mathieu Differential Equation on the Line

9. 4.1 The $\delta$-$\varepsilon$ Plot

10. 4.2 Solutions of the Mathieu Differential Equation

11. 4.3 Explanation on the Behavior of the Solutions

12. 5 The Sierpinski Gasket and the Fractal Laplacian

13. 5.1 Sierpinski Gasket

14. 5.2 Fractal Laplacian

15. 5.3 Spectral Decimation

16. 5.4 Infinite Sierpinski Gasket

17. 6 Extending the Mathieu Differential Equation to Infinite Fractafolds

18. 6.1 Defining the Fractal MDE

19. 6.2 Variants of the Mathieu Differential Equation on the Line

20. 7 Ovservations and Analysis for the Fractal Mathieu Differential Equation

21. 7.1 The $\delta$-$\varepsilon$ Plot

22. 7.2 Solutions of the Mathieu Differential Equation

23. 8 Further Research

24. 9 Appendix

1. Introduction

The Mathieu differential equation, as will be defined in Section 2, takes the form $\frac{d^{2}u}{dt^{2}}+\left(\delta+\varepsilon\cos t\right)u=0$, where $\delta,\varepsilon\in\mathbb{R}$ are fixed, and $u:\mathbb{R}\to\mathbb{R}$. Named after French mathematician Émile Léonard Mathieu (1835-1890), the origin of the Mathieu differential equation stems from real-world phenomena. For example, it describes the motion of a pendulum subject to a periodic driving force. See [21] for more details.

One topic pertaining to the Mathieu differential equation which has been much researched is the stability of solutions. See [2, 10, 16, 17, 28]. Readers may also read [21, 18] for a brief introduction to this area. Most of the existing literature concerns the parameter space, i.e. the space of $\delta$-$\varepsilon$ pairs, whose choices can drastically alter the behavior of solutions. This paper presents research results and new phenomena both on the parameter space and on solutions.

On the other hand, another area of mathematics which has been actively researched in recent years is analysis on fractals, based on J. Kigami’s construction of Laplacians on post-critically finite self-similar sets. See [13, 15, 14, 22]. It is of interest to define and study the Mathieu differential equation on fractal domains, and one suitable domain is the infinite Sierpinski gasket developed by R.S. Strichartz [24]. Existing research on the infinite Sierpinski gasket and on other ‘fractafolds’ be found in [12, 20, 23, 25, 26, 27]. In particular, we will use the important result concerning the spectrum of the Laplacian on the infinite Sierpinski gasket by A. Teplyaev [27]. We will explore a way to generalize the Mathieu differential equation to be defined on the infinite Sierpinski gasket.

This paper is organized as follows. In Section 2, we provide relevant background on the Mathieu differential equation, give definitions which will be used throughout the remainder of the paper, and describe some of the methods used to study the relationship that the values of $\delta$ and $\varepsilon$ have to the solutions of the equation; in doing so we will discuss ‘transition curves’ and the ‘truncation method’ used to study solutions to the Mathieu differential equation. In Section 3 we present a number of proofs for theorems pertinent to the Mathieu differential equation and provide justification for the methods presented in Section 2. Some of our methods and proofs are modified versions from [1, 9]. In Section 4 we provide a discussion of computational results in studying the Mathieu differential equation, including results related to the asymptotic behavior of transition curves and the convergence of solutions. In Section 5 we give an overview of the Sierpinski gasket ($SG$) and provide definitions of various terms in fractal analysis, such as the fractal Laplacian and the infinite Sierpinski gasket ($SG_{\infty}$), to be used in the remaining sections. In Section 6 we extend the content of Sections 2, 3, and 4 to the fractal setting by explaining how the ‘Mathieu differential equation defined on the real line’ can be generalized to a ‘Mathieu differential equation defined on the infinite Sierpinski gasket.’ We describe how solutions to this generalization can be studied by considering solutions expanded as a linear combination of eigenfunctions of the fractal Laplacian. We also describe how the ‘truncation method’ on the line can be used to study the Mathieu differential equation on $SG_{\infty}$. In Section 7 we provide computational results and observations about the shape and asymptotic behavior of the transition curves for the Mathieu differential equation on $SG_{\infty}$ and also about the behavior of solutions. In Section 8 we describe further research that can be done on the Mathieu differential equation and its fractal generalizations. Section 9, the Appendix, describes an alternate approach to the ‘truncation method’ which involves partitioning Fourier coefficients into various equivalence classes.

A website for this research has been created at http://pi.math.cornell.edu/~aac254/. We invite the reader to visit this website, as it contains plots, graphs, data, and other information gathered from the research.

2. Definitions and Methods

In this section we give formal definitions pertaining to the Mathieu differential equation on the real line and then describe the main methods used to study it.

2.1. Fourier Expansion and Matrix Form

Let us suppose we have a periodic solution $u$ with period $2\pi$ or $4\pi$ to the Mathieu differential equation. In this case, we can write the Fourier series expansion of $u$ as

[TABLE]

If we plug this Fourier series for $u$ into the MDE, rearrange terms, and use some trigonometric identities, we find that the function $u$ given above solves the MDE if and only if the two infinite systems of linear homogeneous equations shown below for the cosine and sine coefficients are satisfied. See [18] and [21] for more details.

[TABLE]

and

[TABLE]

To solve these systems of equations, we consider putting them into matrix form. Readers can check that the above two systems of equations are equivalent to the following four equations in matrix form collectively, which correspond to cosine or sine coefficients with even or odd indexes. All the matrices are tridiagonal. For the cosine coefficients with even indexes, we have

[TABLE]

For the sine coefficients with even indexes, we have

[TABLE]

For cosine coefficients with odd indexes, we have

[TABLE]

For sine coefficients with odd indexes, we have

[TABLE]

We can solve the above four matrix equations separately. If either of the first two equations has a nontrivial solution in $\ell^{2}$, then the MDE has a $2\pi$-periodic solution, since

[TABLE]

solves the MDE and is $2\pi$-periodic, comparing with Remark 2.4. Similarly, if either the third or fourth equations has a nontrivial solution in $\ell^{2}$, then the MDE has a $4\pi$-periodic solution, since

[TABLE]

solves the MDE and is $4\pi$-periodic, also comparing with Remark 2.4. Interested readers can also read the appendix for equations in matrix form for periodic solutions with larger periods, say $2N\pi$ with $N$ an arbitrary positive integer. Now we return our discussion to Figure LABEL:deltaepsilonplanefigure, where we plot the stable and unstable regions. According to Theorem LABEL:transitioncurvetheorem, the transition curves consist of $(\delta,\varepsilon)$ pairs with nontrivial $2\pi$- or $4\pi$-periodic solutions. Equivalently, by the above discussion, $(\delta,\varepsilon)$ lies on a transition curve if and only if at least one of the four equations in matrix form discussed above has a nontrivial solution in $\ell^{2}$. Since we develop four different equations in matrix form, we use different colors (orange and black) and line formats (solid and dashed) in Figure LABEL:deltaepsilonplanefigure to distinguish between the four matrix equations:

• •

If $(\delta,\varepsilon)$ falls on a black solid line, then the equation $Ax=0$ has a nontrivial solution.

• •

If $(\delta,\varepsilon)$ falls on a black dashed line, then the equation $Bx=0$ has a nontrivial solution.

• •

If $(\delta,\varepsilon)$ falls on an orange solid line, then the equation $Cx=0$ has a nontrivial solution.

• •

If $(\delta,\varepsilon)$ falls on an orange dashed line, then the equation $Dx=0$ has a nontrivial solution.

2.2. Truncation Method

In finite-dimensional linear algebra, a homogeneous matrix system of equations $Mx=0$, where $M$ is an $m\times m$ matrix, has a nontrivial solution if and only if the determinant of $M$ is zero. Although our four systems of matrix equations above are infinite, we hope to use techniques from finite-dimensional linear algebra to study the properties of these infinite systems. For each infinite matrix, fix $m\in\mathbb{N}$ and consider its $m\times m$ leading principal submatrix. Below we give the truncated matrices $A_{m}$ and $B_{m}$ of $A$ and $B,$ respectively, as examples. $C_{m}$ and $D_{m}$ are defined similarly.

[TABLE]

[TABLE]

Taking the determinant of each truncated matrix yields an algebraic expression involving $\delta$ and $\varepsilon$. Setting each expression equal to zero, we can then plot each equation in the $\delta$-$\varepsilon$ plane and obtain a set of algebraic curves in the variables $\delta$ and $\varepsilon$. If we choose $m$ sufficiently large, the $\delta$-$\varepsilon$ curves derived from the truncated matrices will be very close to the true transition curves corresponding to the infinite matrices. This statement is made precise by Ikebe et al. in [1], and we present their results in Section 3.2.

3. Three-Term Recurrence Relations and Truncations of the Infinite Matrices

In this section, we will provide theoretical foundation for the methods we use. In doing so we will discuss the asymptotic convergence of Fourier coefficients and provide an error estimate for the truncation method we use for approximating the transition curves. (We will extend these techniques to the fractal MDE in later sections.) In this section, we work in generality by considering matrices of the form

[TABLE]

where $\alpha_{j},\beta_{j},\lambda_{j}\in\mathbb{R}$ for $j\geq 1$, and $\gamma\in\mathbb{R}$; further, we assume that $\{\alpha_{j}\}_{j=1}^{\infty}$ and $\{\beta_{j}\}_{j=1}^{\infty}$ are bounded sequences and that $\lim_{j\to\infty}\lambda_{j}=\infty$. Note that matrices $A,B,C,$ and $D$ as defined in Section 2.2 are special cases of (3.1). In addition, we always assume $\alpha_{j},\beta_{j}\neq 0$, for all $j\geq 0$, in this section. We will consider equations of the form

[TABLE]

where $\bm{c}=(c_{1},c_{2},\cdots)^{T}\in\ell^{2}:=\left\{(x_{1},x_{2},...)^{T}:\sum_{j=1}^{\infty}\left|x_{j}\right|^{2}<\infty\right\}$ is a real sequence.

3.1. Asymptotic Behavior of $\{c_{j}\}_{j=1}^{\infty}$

Three-term recurrence relations (TTRRs) and difference equations have been well-studied throughout the last century. Since (3.2) naturally gives a TTRR,

[TABLE]

we can use existing theorems to study the asymptotic behavior of the sequence $\{c_{j}\}_{j=1}^{\infty}$ as $j\to\infty$. Below we first present the well-known Poincaré Theorem (see [9]), also called the Poincaré-Perron Theorem, which describes the asymptotic behavior of solutions to general TTRRs, and then we will see how the theorem applies to the equations above.

Theorem 3.1** (Poincaré-Perron Theorem, [9]).**

For $j\geq 1$, let $c_{j+1}+b_{j}c_{j}+a_{j}c_{j-1}=0$ with $\lim_{j\to\infty}b_{j}=b$ and $\lim_{j\to\infty}a_{j}=a$. Let $t_{1}$ and $t_{2}$ denote the zeros of the characteristic equation $t^{2}+bt+a=0$. Then, if $|t_{1}|\neq|t_{2}|$, the difference equation has two linearly independent solutions $\{x_{n}\}$ and $\{y_{n}\}$ satisfying

[TABLE]

If $|t_{1}|=|t_{2}|$, then

[TABLE]

for any nontrivial solution $\{c_{j}\}$ to the recurrence relation. ∎

With this theorem in mind, we consider the TTRR

[TABLE]

where $\{f_{j}\}$ and $\{g_{j}\}$ neither of which contains [math] as one of its terms, are uniformly bounded, and where $|d_{j}|\to\infty$ as $j\to\infty$. Notice that, with our assumptions, the solution of Equation 3.4 is uniquely determined by specifying the first two terms $c_{1}$ and $c_{2}$. We can derive the following proposition from the Poincaré-Perron Theorem, adapting its proof directly from the proof of the Poincaré-Perron Theorem in [9].

Proposition 3.2**.**

Assume $\{f_{j}\}$ and $\{g_{j}\}$, neither of which is equal to the zero sequence, are uniformly bounded, and assume that $|d_{j}|\to\infty$ as $j\to\infty$. Then the TTRR (3.4) has two linearly independent solutions $\{x_{j}\}$ and $\{y_{j}\}$ satisfying

[TABLE]

Proof. The proof is modified directly from the proof for the Poincaré-Perron Theorem in [9]. For convenience, we assume $d_{j}\neq 0$ for any $j$, since otherwise we can consider the TTRR defined for $j\geq M$ for some large $M$. Let $\hat{c}_{j}=\left(\prod_{i=1}^{j-1}\frac{f_{i}}{d_{i}}\right)c_{j}$. Then the TTRR becomes

[TABLE]

By the Poincaré-Perron theorem, there exists a pair of linearly independent solutions of this new recurrence, $\{\hat{x}_{j}\}$ and $\{\hat{y}_{j}\}$, satisfying

[TABLE]

For $j$ sufficiently large, $\hat{x}_{j}$ does not vanish. Indeed, if $\hat{x}_{j_{0}}=0$ for some $j_{0}\in\mathbb{N}$ then the TTRR implies $\frac{\hat{x}_{j_{0}+2}}{\hat{x}_{j_{0}+1}}=-1$, which cannot happen if $j_{0}$ were sufficiently large (since $\frac{\hat{x}_{j}}{\hat{x}_{j-1}}\to 0$). Similarly, for $j$ sufficiently large $\hat{y}_{j}$ does not vanish, for if $\hat{y}_{j_{0}}=0$ for some $j_{0}\in\mathbb{N}$ then $\frac{\hat{y}_{j_{0}}}{\hat{y}_{j_{0}-1}}=0$, which cannot happen for $j_{0}$ sufficiently large (since $\frac{\hat{y}_{j}}{\hat{y}_{j-1}}\to-1$). Let $H^{(1)}_{j}=\hat{x}_{j}/\hat{x}_{j-1}$ and $H^{(2)}_{j}=\hat{y}_{j}/\hat{y}_{j-1}$. Then

[TABLE]

Multiplying each side by $H_{j}^{(l)},l\neq i$, we get

[TABLE]

Subtracting one equation from the other, we obtain

[TABLE]

So, since

[TABLE]

we have

[TABLE]

Therefore,

[TABLE]

and

[TABLE]

Let $\{z_{j}\}$ be an arbitrary solution to the TTRR 3.4. We say $\{z_{j}\}$ is a minimal solution of the TTRR if it obeys $\frac{z_{j}}{z_{j-1}}\sim-\frac{g_{j-1}}{d_{j}}$, and we say $\{z_{j}\}$ is a dominant solution of the TTRR if it satisfies $\frac{z_{j}}{z_{j-1}}\sim-\frac{d_{j-1}}{f_{j-1}}$. Note that $\lim_{j\to\infty}z_{j}=0$ if $\{z_{j}\}$ is a minimal solution and $\lim_{j\to\infty}z_{j}=\infty$ if $\{z_{j}\}$ is a dominant solution. We can apply Proposition 3.2 to study Equation 3.3, since Equation 3.3 consists of an equation for $c_{1},c_{2}$ and a TTRR of the form 3.4.

Corollary 3.3**.**

Assume $\varepsilon\neq 0$. (a). Equation 3.2 has a unique solution up to multiplication by a constant. (b). Let $\bm{c}=(c_{1},c_{2},c_{3},...)^{T}$ be a nontrivial solution to Equation (3.2). If $\bm{c}\in\ell^{2}$, then $c_{j}$ will decay with the asymptotic behavior $\frac{c_{j}}{c_{j-1}}\sim\frac{\beta_{j-1}\varepsilon}{\lambda_{j}-\delta}$. If $\bm{c}\notin\ell^{2}$, then $c_{j}$ will diverge with the asymptotic behavior $\frac{c_{j}}{c_{j-1}}\sim\frac{\lambda_{j-1}-\delta}{\alpha_{j-1}\varepsilon}$.

Proof. It is equivalent to studying Equation 3.3. (a). Given $c_{1}\in\mathbb{R}$, we can solve $c_{2},c_{3},\cdots$ inductively, which gives us a solution to Equation 3.3. In this way, any solution to Equation 3.3 is uniquely determined by $c_{1}$. (b) Any solution $\textbf{c}=\{c_{j}\}_{j=1}^{\infty}$ can be written as a linear combination of $\{x_{j}\}$ and $\{y_{j}\}$ as given in in Proposition 3.2, and hence there exist $p,q\in\mathbb{R}$ such that $c_{j}=px_{j}+qy_{j}$ for all $j$. Note that if $q\neq 0$ then $\{z_{j}\}$ is a dominant solution, and if $q=0$ then $\{z_{j}\}$ is a minimal solution. If $\{c_{j}\}\in\ell^{2}$ then $\lim_{j\to\infty}c_{j}=0$ and thus $\{c_{j}\}$ cannot be a dominant solution; thus, we must have $p=0$ and hence $\frac{c_{j}}{c_{j-1}}\sim\frac{\beta_{j-1}\varepsilon}{\lambda_{j}-\delta}$. On the other hand, if $\{c_{j}\}\notin\ell^{2}$ then $p\neq 0$ and hence $\frac{c_{j}}{c_{j-1}}\sim\frac{\lambda_{j-1}-\delta}{\alpha_{j-1}\varepsilon}$. $\square$

Corollary 3.3 shows that there is a sharp contrast between the two possible types of behavior that a solution to Equation (3.2) can have. Since $\{\alpha_{j}\}$ and $\{\beta_{j}\}$ are bounded sequences, when $(\delta,\varepsilon)$ are fixed and properly chosen any solution $\{c_{j}\}$ to (3.2) will converge to [math] very rapidly; otherwise, all nontrivial solutions will tend to infinity very fast. If $c_{1}\in\mathbb{R}$ is an initial value which corresponds to a minimal solution, Proposition 3.2 shows that numerically computing the solution with the forward recursion method is unstable, since a small error in computation will lead to a dominant solution, in which $\{c_{j}\}$ explodes. On the other hand, the backward recursion method can give a good approximation for a minimal solution. In [9], detailed discussions are given on this method, as well as a corresponding error estimate. We briefly state the result here.

Proposition 3.4** ([9]).**

Fix a large $N$, and set $c_{N+1}=0,c_{N}=1$. Compute $c_{N-1},c_{N-2},\cdots,c_{1}$ one-by-one backward with TTRR (3.4). If $\{x_{k}\}$ is an arbitrary minimal solution and $\{y_{k}\}$ is an arbitrary dominant solution, then $\{c_{j}\}$ satisfy

[TABLE]

where $r_{k}=\frac{x_{k}}{y_{k}}$ and $1\leq k\leq N$.

3.2. The Truncation Method

In this part, we introduce the truncation method for determining which $(\delta,\varepsilon)$ pairs yield nontrivial solutions of (3.2). We start from the observation

[TABLE]

where $I$ is the identity matrix and the shorthand notation $T(\varepsilon)$ is used in the last equality. Clearly, for a fixed $\varepsilon$, equation (3.2) has nontrivial solutions if and only if $\delta$ is an eigenvalue of $T(\varepsilon)$. The eigenvalue problem for infinite tridiagonal matrices has been widely studied. In [1], there is a result on the error between eigenvalues of the truncated matrices and those of the corresponding infinite matrix. We state the result in Theorem 3.5 below. Here we consider the infinite symmetric tridiagonal matrix $T$ of the form

[TABLE]

where $d_{n}\to\infty$ as $n\to\infty$ and $\{f_{n}\}$ bounded. Denote its $n\times n$ truncation by $T_{n}$, i.e.,

[TABLE]

Theorem 3.5**.**

*([1]) Let $T$ and $T_{n}$ $(n\geq 1)$ be given as above.

(a). $T$ has pure point spectrum.

(b). If $\delta$ is a given simple eigenvalue of $T$, then there exists, for each $n\in\mathbb{N}$, an eigenvalue $l_{n}$ of $T_{n}$ such that the sequence $\{l_{n}\}_{n=1}^{\infty}$ satisfies $l_{n}\to\delta$ as $n\to\infty$. For any such sequence the error is given by*

[TABLE]

where $\bm{c}=(c_{1},c_{2},\cdots)^{T}\in\ell^{2}$ is an eigenvector corresponding to $\delta$. ∎

We can extend the above result to nonsymmetric tridiagonal matrices.

Theorem 3.6**.**

Let $T$ be an infinite tridiagonal matrix of the form

[TABLE]

where $d_{n}\to\infty$ as $n\to\infty$, and $\{f_{n}\}_{n=1}^{\infty}$ and $\{g_{n}\}_{n=1}^{\infty}$ are bounded, positive, and nonzero. Let $T_{n}$ be the $n$ $\times$ $n$ truncation of $T$. If $\delta$ is a given simple eigenvalue of $T$, then there exists, for each $n\in\mathbb{N}$, an eigenvalue $l_{n}$ of $T_{n}$ such that the sequence $\{l_{n}\}_{n=1}^{\infty}$ satisfies $l_{n}\to\delta$ as $n\to\infty$. For any such sequence the error is given by

[TABLE]

where $\kappa_{j}:=\sqrt{\prod_{i=1}^{j-1}\frac{g_{i}}{f_{i}}}$ and where $\bm{c}=(c_{1},c_{2},\cdots)^{T}\in\ell^{2}$ is an eigenvector corresponding to $\delta$. ∎

Proof. (a) Let

[TABLE]

We observe that

[TABLE]

It is clear that $T^{\prime}$ is a self-adjoint operator from $\ell^{2}$ to $\ell^{2}$ with pure point spectrum by Theorem 3.5. Let

[TABLE]

where $\bm{c}=(c_{1},c_{2},\cdots)^{T}$. Clearly, the eigenvalue equation $(\delta I-T)\bm{c}=0$ holds if and only if $(\delta I-T^{\prime})\bm{c}^{\prime}=0$ in a pointwise sense, where $I$ is the identity transformation. In addition, from Corollary 3.3 we deduce that there is a one-to-one correspondence between eigenvectors of $T$ and $T^{\prime}$ in $\ell^{2}$. Indeed, if $\bm{c}$ is an eigenvector of $T$ in $\ell^{2}$, we have $\frac{c^{\prime}_{j}}{c^{\prime}_{j-1}}=\frac{\kappa^{-1}_{j}c_{j}}{\kappa^{-1}_{j-1}c_{j-1}}\sim-\frac{\sqrt{f_{j-1}g_{j-1}}}{d_{j}}$ because $\frac{c_{j}}{c_{j-1}}\sim-\frac{g_{j-1}}{d_{j}}$ by Corollary 3.3, which implies $\bm{c}^{\prime}$ is an eigenvector of $T^{\prime}$ in $\ell^{2}$, and for the same reason if $\bm{c}^{\prime}\in\ell^{2}$ is an eigenvector of $T^{\prime}$ then $\bm{c}=diag(\kappa_{1},\kappa_{2},\cdots)\bm{c^{\prime}}$ is an eigenvector of $T$ in $\ell^{2}$. Now, using Theorem 3.5 we deduce that there is a sequence of eigenvalues of the truncated $T^{\prime}_{n}$ that converges to the eigenvalue $\delta$ of $T^{\prime}$ and hence of $T$, since $T$ and $T^{\prime}$ have the same eigenvalues. In addition, noticing that, for each $n\in\mathbb{N}$, $T_{n}$ has the same eigenvalue as $T^{\prime}_{n}$, we see that there is a sequence of eigenvalues of $T_{n}$ which converges to $\delta$. Lastly, the error estimate comes by applying Theorem 3.5 to $T^{\prime}$ with the eigenvector $\bm{c}^{\prime}=\text{diag}(\kappa^{-1}_{1},\kappa^{-1}_{2},\cdots)\bm{c}$. $\square$

4. Observations and Analysis for the Mathieu Differential Equation on the Line

The Mathieu differential equation has been an important topic in differential equations due to its numerous real-world applications. However, most of the existing work on the MDE focuses on theoretical analysis of the stable and unstable regions of the $\delta$-$\varepsilon$ plane, such as the asymptotic behavior of the transition curves, and there has not been much computational work done on the stability curves as well as on the solutions to the Mathieu differential equation. In this section, in addition to describing some of our theoretical results, we will explain our computational results on the intricate shape of the $\delta$-$\varepsilon$ curves and on the converging behavior of the solutions.

4.1. The $\delta$-$\varepsilon$ Plot

In this part, we use the truncation method introduced in Section 2 and Section 3 to study the $\delta$-$\varepsilon$ curves in more detail. As discussed above, the transition curves are found via $2\pi$-periodic and $4\pi$-periodic Fourier expansions of solutions. However, one might wonder what the curves would be if the same process is undertaken for Fourier expansions of larger periods, say periods $8\pi,16\pi,32\pi,...$. The answer is that curves corresponding to expansions of larger periods in fact ”fill in” the stable regions whose boundary is obtained using the $2\pi$- and $4\pi$-periodic expansions. See Figure 4.1 for an illustration.

In fact, this “fill in” property remains valid for even larger periods, and a proof is given below.

Proposition 4.1**.**

Let $D$ be the set of $\delta$-$\varepsilon$ pairs such that the corresponding MDE has a solution of period $2^{k}\pi$ for some $k\geq 3$. Then $D$ is dense in the stable region.

Proof. First, if $(\delta,\varepsilon)$ is a stable pair, then by the discussion given in Chapter VI of [21] we can find two linearly independent solutions $u_{1}$ and $u_{2}$ of the MDE, such that

[TABLE]

for some $\theta\in[0,\frac{1}{2}]$. $\theta$ is uniquely determined by the pair $(\delta,\varepsilon)$, so we may view $\theta=\theta(\delta,\varepsilon)$ as a function from the stable region to $[0,1/2]$. In fact, by the Floquet Theorem, there exists a $2\times 2$ matrix $A$ depending only on $(\delta,\varepsilon)$ such that, for any solution $u$ of the MDE,

[TABLE]

So $(\delta,\varepsilon)$ is a stable pair only if the two eigenvalues of $A$ both have norm $1$, which can be written as $e^{i2\pi\theta}$ and $e^{-i2\pi\theta}$. Readers can find details in [21]. Clearly, $A$ depends continuously on the parameter $(\delta,\varepsilon)$, which shows that $\theta=\theta(\delta,\varepsilon)$ is a continuous function on the stable region. Next, we fix $\theta_{0}\in[0,\frac{1}{2}]$ and $\varepsilon\in\mathbb{R}$, and show that only countably many $\delta$’s can be found such that $\theta_{0}=\theta(\delta,\varepsilon)$. Recall that for each such $\delta$ we have a solution $u_{1}$ to the corresponding MDE so that $e^{-i\theta_{0}t}u_{1}(t)$ is a $2\pi$-periodic function. We have the Fourier series expansion $e^{-i\theta_{0}t}u_{1}(t)=\sum_{j=-\infty}^{\infty}c_{j}e^{ijt}$, which implies

[TABLE]

Plugging this into the Mathieu differential equation, we get the following equation for the coefficients,

[TABLE]

From this we can see that there are only countably many values of $\delta$ such that the corresponding equation has a nontrivial solution $\{c_{j}\}_{j=-\infty}^{\infty}\in\{\bm{c}\in\mathbb{C}^{\mathbb{Z}}|\sum_{j=-\infty}^{\infty}|c_{j}|^{2}<\infty\}$. Now, we prove the proposition by contradiction. Suppose $D$ is not dense in the stable region. Then, there is a ball in the $\delta$-$\varepsilon$ plane where $\theta$ does not take any dyadic rational value, noticing that when $\theta$ is a dyadic rational the MDE has $2^{k}\pi$-periodic solutions $u_{1},u_{2}$ for some $k$. This means $\theta$ equals a constant $\theta_{0}$ on that ball, since $\theta$ is a continuous function of $(\delta,\varepsilon)$. However, for this fixed $\theta=\theta_{0}$, and a fixed $\varepsilon$ in the ball, there are only countably many possible values for $\delta$, which can not fill in the ball. This gives a contradiction. $\square$

Next we will discuss the asymptotic behavior of the stable and unstable regions. Existing works including [10], [16], and [28] cover two important observations:

1. (i)

First, as illustrated in Figure 4.2, the width of each stable band gets thinner and thinner as $|\varepsilon|$ tends to $\infty$. In [28], it is shown that the width of the $k$th stable band will decrease exponentially on $|\varepsilon|\geq k^{2}$ as $|\varepsilon|\to\infty$.

2. (ii)

Secondly, for each fixed $\varepsilon$, the width of $k$th unstable band becomes smaller as $k$ increases. An estimate is given in [2] and [10] with different methods, both yielding

[TABLE]

where $d_{k}$ is the width of the $k$th unstable band with $\varepsilon$ fixed.

Since the study of the type of convergence in the second observation is comparatively complete, we focus on the type in the first observation, which concerns the width of the stable band as $|\varepsilon|\to\infty$. In particular, we study the first ten stable bands by computing the width of each band from $\varepsilon=0$ to $\varepsilon=50$ in increments of $0.1$. A graph for each stable band, plotting width vs epsilon, are shown below. Then we use the curve fitting toolbox in MATLAB to estimate the fitting curves between $\varepsilon$ and width in each stable band. Please see Figure 4.3 and Figure 4.4 for the result and fitting curves.

The last topic we will discuss before moving on to talk about solutions to the MDE is to approximate and simplify the irregular boundary between stable and unstable regions for practical use. We can see that most $(\delta,\varepsilon)$ pairs in the first and the fourth quadrants which lie below by the line $\varepsilon=\delta$ and above by the line $\varepsilon=-\delta$ are stable, while most $(\delta,\varepsilon)$ values in those quadrants outside that region form unstable pairs. For any $w>0$, let $R_{w}:=\{(\delta,\varepsilon)\in\mathbb{R}^{2}:0<\varepsilon<w\text{ and }-\delta<\varepsilon<\delta\}$. We are interested in determining the probability, for various values of $w$, that a $\delta$-$\varepsilon$ pair is stable, given that it lies in the triangle $R_{w}.$ So, we numerically compute the probabilities with different choices of $w$. Notice that the line $\varepsilon=\delta$ (as well as the line $\varepsilon=-\delta$) alternatively passes through stable and unstable regions. So, we define $w_{i}$ to be the $\delta$-coordinate of the point of intersection between the line $\varepsilon=\delta$ and the right boundary of the $i$-th unstable region; also, for each $i,$ we define $P_{i}$ to be the probability that a $\delta$-$\varepsilon$ pair is stable, given that it is in $R_{w_{i}}$. Readers can see figure 4.5 for an illustration of how we divide the regions into triangles.

$P_{i}$ can also be considered as the probability of getting a stable $(\delta,\varepsilon)$ pair within the triangular regions, which can characterize how well the irregular boundaries can be approximated by lines $\varepsilon=\pm\delta$. We list the $P_{i},1\leq i\leq 18$ in Table 1.

In addition, we also obtain a surprisingly good fitting curve of the form $P_{i}=\frac{ai+b}{i+c}$, with $a=0.9938,b=-0.1424,c=0.3608$. See figure 4.6 for details.

4.2. Solutions of the Mathieu Differential Equation

Now we discuss the periodic solutions for the Mathieu differential equation. In Corollary 3.3, we showed that the sequence $\{c_{j}\}$ of Fourier coefficients to periodic solutions corresponding to $\delta$-$\varepsilon$ pairs on the transition curves converge rapidly. Hence, given a $\delta$-$\varepsilon$ pair on a transition curve, we can calculate explicitly the coefficients of the corresponding solution with a finite truncation of the Fourier series and use those coefficients to plot, with very high accuracy, solutions to the MDE. One question one might investigate is as follows. Suppose one fixes a transition curve and considers periodic solutions corresponding to various $(\delta,\varepsilon)$ points along that transition curve. How do properties of solutions change as $\varepsilon$ varies? To be more precise, we need some new notation. Recall from Section 2.2 that, when solving $2\pi$- or $4\pi$- periodic solutions, we can rewrite the MDE as a system of linear equations for Fourier coefficents, which collectively can be rewritten as four independent homogeneous matrix equations, with corresponding matrices $A,B,C,D$. The transition curves consist of the $\delta$-$\varepsilon$ pairs such that one of these matrices degenerates. In this way, a transition curve can be labeled with a matrix and an integer $k\geq 1$, so that for each $\delta$-$\varepsilon$ pair on this transition curve, $\delta$ is the $k$th smallest real number such that the matrix degenerates with the fixed $\varepsilon$. We want a more convenient way to label these transition curves, and to achieve this, we proceed as follows. Fix a transition curve, let $\varepsilon=0$ on the transition curve, and solve the corresponding equation in matrix form. According to the discussion in Section 2.2, we can then construct a periodic solution to the MDE, and the solution is clearly an eigenfunction of $\frac{d^{2}}{dt^{2}}$. For example, if we consider the transition curve labeled by $A$ and $k=2$ we will get $\cos t$ as a solution with the above process. In this way, a transition curve can be labeled with an eigenfunction of $\frac{d^{2}}{dt^{2}}$. In light of this, we define $p(\varphi,\varepsilon)$ to be the point in the $\delta$-$\varepsilon$ plane with prescribed $\varepsilon$-coordinate which lies on the transition curve corresponding to the eigenfunction $\varphi$. As shown in Section 2, the transition curves can be organized into four different classes, made up of points $\{p(\sin kt,\varepsilon)\},\{p(\cos kt,\varepsilon)\},\{p(\sin\frac{2k+1}{2}t,\varepsilon)\},\{p(\cos\frac{2k+1}{2}t,\varepsilon)\}$ with our new notation. It is natural to study them separately. First we consider the class $\{p(\sin kt,\varepsilon)\}$. Because of the symmetry of the transition curves across the $\delta$-axis, it suffices to only consider positive values of $\varepsilon$. In fact, if $u$ is a solution corresponding to the pair $(\varepsilon,\delta)$, then $u(t+\pi)$ is a solution corresponding to the pair $(-\varepsilon,\delta)$ since $-\cos t=\cos(t+\pi)$. The normalized solutions are plotted in Figures 4.7-4.9, where by ‘normalized solution’ we mean a solution $u$ for which $\max_{t\in\mathbb{R}}u(t)=1$. In plotting the graphs we have used the particular values $k=1,2,3$ and $\varepsilon=5,10,20,40,80,160$. $15\times 15$ truncated matrices are used to compute the solutions. The horizontal axis is the $t$-axis, and the vertical axis is the $u$-axis.

An interesting pattern can be observed from Figure 4.7. The sequence of $t$-coordinates of these maximal points appears to approach [math] monotonically as $\varepsilon$ increases. Figures 4.8 and 4.9 show the same behavior as well. We now use curves to fit a curve of the $t$-coordinate of these relative maxima in $[0,\pi]$ as a function of $\varepsilon$. We still pick the same three transition curves, and take $\varepsilon$ from $1$ to $200$, in increments by $1$. The relationship of $t$ coordinate of maximal points and $\varepsilon$ on the three curves are shown in figure 4.10, and the fitting curve is of the form $t=\frac{a\varepsilon+b}{\varepsilon^{2}+c\varepsilon+d}$.

For the case of $p(\sin 2t,\varepsilon)$ and $p(\sin 3t,\varepsilon)$, does the same convergence behavior occur if one considers the sequence of the first minima when $t>0$? We also do the computation, and show the results in figure 4.11, where the same kind of fitting curve also works.

In case of $p(\sin 3t,\varepsilon)$, the same observation applies to the second maxima when $t>0$. See figure 4.12 for details.

One can also consider the behavior of the $u$-coordinate of these various sequences of local extrema. See figure 4.13 for the behavior of $u$-coordinate of the first minimal points on the curve $p(\sin 2t,\varepsilon)$, where we can fit the points well with a rational function of the form $\frac{a\varepsilon^{2}+b\varepsilon+c}{\varepsilon^{2}+b\varepsilon+c}$.

For the curve $p(\sin 3t,\varepsilon)$, there are three local extremes on each solution, and we also do the computation for the $u$-coordinate for the first minimal and the second maximal value. The fitting curve is a little more complicated, of the form $t=\frac{a\varepsilon^{3}+b\varepsilon^{2}+c\varepsilon+d}{\varepsilon^{3}+d\varepsilon^{2}+e\varepsilon+g}$. See figure 4.14.

Next, we turn to the solutions on curves $p(\cos kt,\varepsilon)$, and we do experiments for $k=0,1,2$. Please look at figure 4.15,4.16,4.17 for the results.

We can see that the shapes of these curves are different from those in figure 4.7 4.8,4.9. The first curve $p(\cos 0t,\varepsilon)$ is of course very special, as it corresponds to the constant functions. It turns out the solution is always positive on this curve, as a consequence of Sturm’s Theorem [17]. Also, we have the minimal values of the solutions at $\pi$, see figure 4.18 for the minimal values for different choices of $\varepsilon$. We use a fitting curve of the form $u=\frac{a\varepsilon+b}{\varepsilon^{2}+c\varepsilon+d}$.

The solutions on the second curve is a little more complicated. As shown in figure 4.16, we can see that minimal value is achieved at $t=\pi$ in $[0,2\pi]$ for $\varepsilon=0$ and $\varepsilon=1$, but as $\varepsilon$ becomes larger, the minimum point splits into two different minimum points. This is an interesting phenomenon, and we have an explanation for this in Theorem LABEL:thm44. We do computations on the $t$-coordinate of the minimum points as $\varepsilon$ gets larger. Please see figure 4.19 for the data, and we use a fitting curve of the form $t=\frac{a\varepsilon+b}{\varepsilon^{2}+c\varepsilon+d}$.

Another thing we need to highlight here is that the maximal absolute value does not occur at [math] or $\pi$, but at some other points. So, as usual, we normalize the solutions so that their minimum value are $-1$, and we compute values of solutions at the two other local peaks at [math] and $\pi$, which is shown in figure 4.20. It is also interesting to see in figure 4.21 that the local maximal value at [math] drop rapidly near $\varepsilon=0$, and then increase quickly.

The solutions on the curve $p(\cos 3t,\varepsilon)$ is even more complicated. Still, we observe the split of one peak. When $\varepsilon=0,1,2,3,4,5$, we only see one maximum point, but as $\varepsilon$ becomes larger, we see two maximum points. We also do computation on the $t$ coordinate of one of these splited maximum point. See figure 4.22 for the computation result, where we use fitting curves of the form $t=\frac{a\varepsilon^{2}+b\varepsilon+c}{\varepsilon^{2}+d\varepsilon+e}$.

We also do experiment on the t coordinate of the first minimum points and the $u$-coordinate of the peaks on the curve $p(\cos 2t,\varepsilon)$.

4.3. Explanation on the Behavior of the Solutions

Before the end of this section, we give a proof on the convergence of the t-coordinates of the peaks. This interesting point is that the behaviors are closely related to the asymptotic behavior of the transition curves, i.e. how $\delta$ behaves as $\varepsilon$ increases to $\infty$. There are several works on this problem, and here we refer [17] for one version. Readers can find other dicussions in [28].

Theorem 4.2** (W.S. Loud).**

For fixed $\varepsilon>0$, let $\delta$ be the $k$th smallest value on the transition curves. Then we have

[TABLE]

as $\varepsilon\to+\infty$. ∎

Notice that in the above theorem, we fix a transition curve and let $\varepsilon\to\infty$, which is exactly the way we study the behavior of the solutions. Now, we apply the theorem to derive some interesting facts that we have observed in our experiments.

Lemma 4.3**.**

For fixed $\varepsilon>0$, let $\delta$ be the $k$th smallest value on the transition curves, and $u$ be a nontrivial periodic solution for the corresponding MDE. Then, for large $\varepsilon>0$, there is no local extremum of $u$ in $[c_{\varepsilon},\pi)$ where $c_{\varepsilon}\sim\varepsilon^{-1/4}$.

Proof. Let $c_{\varepsilon}$ be the zero of $\delta+\varepsilon\cos t$ in $(0,\pi)$, where theorem 4.2 guarantee the existence of $c_{\varepsilon}$ when $\varepsilon$ is large. In addition, we have the estimate by using Taylor expansion of $\cos t$ locally at [math],

[TABLE]

(b) We only need to understand the behavior of $u$ at $\pi$. The idea is essentially the same as the proof of Lemma 4.3, where we look at the sign of $\delta+\varepsilon\cos t$. It is also clear that $u(\pi)\neq 0$ and $u^{\prime}(\pi)=0$ in this case. Notice that $u^{\prime\prime}(\pi)=-(\delta-\varepsilon)u(\pi)$, and we set $u(\pi)>0$ by multipling a constant, so that we can conveniently look at the absolute value of $u$. If $\delta+\varepsilon\cos\pi=\delta-\varepsilon>0$, then $\pi$ is a local maximum of $|u|$; if $\delta-\varepsilon<0$, then $\pi$ is a local minimum of $|u|$. For the case that $\varepsilon=\delta$, we can see that $\delta-\varepsilon\cos t>0$ on $(0,2\pi)\setminus\{\pi\}$ except for the case $\delta=\varepsilon=0$, which only happens on the curve characterized by $\cos 0t$. So $\pi$ is still a local maximum of $|u|$. In addition, $\delta-\varepsilon$ is a strictly decreasing function of $\varepsilon$ on $(0,\infty)$, so we can find a critical point $\alpha$ such that $\delta-\varepsilon>0$ on $(0,\alpha)$, and $\delta-\varepsilon<0$ on $(\alpha,\infty)$. To show that $\delta-\varepsilon$ is strictly decreasing, notice that $\delta$ is the $k$th smallest value such that the MDE has a $2\pi$- or $4\pi$-periodic solution, which is equivalent to say that $\delta$ is the $k$th smallest eigenvalue of the self-adjoint operator $-\Delta-\varepsilon\cos t$ on $L^{2}(\mathbb{R}/4\pi\mathbb{Z})$. Let $\varepsilon>\varepsilon^{\prime}>0$, and $\delta,\delta^{\prime}$ be the corresponding eigenvalues. Then using Raleigh quotient, we get

[TABLE]

where $M_{k,\varepsilon^{\prime}}$ is the space spanned by the smallest $k$ eigenfunctions of $-\Delta-\varepsilon^{\prime}\cos t$ and the infimum in the last line is taken over all the $k$ dimensional subspaces of $dom\Delta$ on $\mathbb{R}/4\pi\mathbb{Z}$. $\square$

Remark. In Theorem 4.2, Lemma 4.3 and Theorem 4.4, we considered the case that $\varepsilon\geq 0$. Since $-\cos t=\cos(t+\pi)$, for the case that $\varepsilon\leq 0$, the MDE can be rewritten as

[TABLE]

with $\hat{u}(t)=u(t-\pi)$. We can still use Theorem 4.2, Lemma 4.3 and Theorem 4.4 to study solutions to the MDE with coefficient $\varepsilon<0$, by applying a shift of $\pi$.

5. The Sierpinski Gasket and the Fractal Laplacian

In this section we give preliminary definitions and results concerning analysis on fractals and provide a basic introduction to the ‘infinite’ Sierpinski gasket will be given. This serves to set up the discussion of the generalization of the MDE to the fractal setting described in in Section 6.

5.1. Sierpinski Gasket

Consider the three contraction mappings $\{F_{i}:\mathbb{R}^{2}\to\mathbb{R}^{2}\}_{i=0,1,2}$ given by

[TABLE]

where $x\in\mathbb{R}^{2}$. Then $\{F_{i}\}_{i=0,1,2}$ form an ‘iterated function system’ (see page 133 of [7]). By Theorem 9.1 in [7], there exists a unique nonempty compact set $K\subset\mathbb{R}^{2}$ such that (see [11])

[TABLE]

Then $K$ is defined to be the Sierpinski gasket, often denoted $SG$. See Figure 5.1.

In studying $SG$ it is useful to use its graph approximations, constructed as follows. Let $q_{0},q_{1},$ and $q_{2}$ be the unique fixed points of $F_{0},F_{1},$ and $F_{2},$ respectively. Define a vertex set $V_{0}:=\{q_{0},q_{1},q_{2}\}\subset SG$. Then $F_{i}(SG)\cap F_{j}(SG)=F_{i}V_{0}\cap F_{j}V_{0}$,for any $i\neq j$. We refer to $V_{0}$ as the boundary of $SG$. We further define vertex sets $V_{n}\subset SG$ for $n\geq 1$ inductively by $V_{n}:=\bigcup_{i=0}^{2}F_{i}(V_{n-1})$ and let $V_{*}:=\bigcup_{n=0}^{\infty}V_{n}$ be the set of all vertices. Note that $V_{*}$ is dense in $SG$. For an $m$-tuple $w=\left(w_{1},w_{2},\dotsb,w_{m}\right)$, where $w_{j}\in\{0,1,2\}$ for each $w_{j}\text{ }(1\leq j\leq m)$, we define $F_{w}$ by

[TABLE]

and say that $w$ is a word of length $\left|w\right|=m$. With this, an edge relation $\sim_{m}$ on $V_{m}$ can be introduced as follows: for $x,y\in V_{m}$ we say $x\sim_{m}y$ if and only if there exists a word $w$ of length $m$ and unequal indices $i,j\in\{0,1,2\}$ such that $x=F_{w}(q_{i})$ and $y=F_{w}(q_{j})$. This relation on $V_{m}$ gives a sequence of graphs $\Gamma_{m}$ approximating SG, with the vertex set $V_{m}$ and the edge set $E_{m}=\{\{x,y\}|x\sim_{m}y\}$. See Figure 5.2 for $\Gamma_{0},\Gamma_{1},\Gamma_{2}$.

5.2. Fractal Laplacian

Now we are ready to define the Laplacian on $SG$. Suppose $u:SG\to\mathbb{R}.$ We define the level-$m$ discrete Laplacian $\Delta_{m}$ by

[TABLE]

Then we define the continuous Laplacian $\Delta$ by

[TABLE]

If the limit above converges uniformly on $V_{*}\setminus V_{0}$ to a continuous function, we say $u\in\text{dom}\Delta$ . In this case, we extend $\Delta u$ to all of $SG$, including points not in $V_{*}\setminus V_{0}$, by continuity (recall that $V_{*}\setminus V_{0}$ is dense in $SG$). The continuous Laplacian on SG is the analog of the usual ‘second-order derivative’ on the line. We shall mention now the the following proposition derived by O. Ben-Bassat, R.S. Strichartz,and A. Teplyaev in [3], which will be of interest in the next section.

Theorem 5.1** (O. Ben-Bassat, R.S. Strichartz & A. Teplyaev).**

Let $u$ be a nonconstant function in dom$\Delta.$ Then $u^{2}$ is not in dom$\Delta$. ∎

A function $u$ satisfying $-\Delta u=\lambda u$ for some number $\lambda$ is called an eigenfunction of $\Delta$ with eigenvalue $\lambda$. If a function $u:SG\to\mathbb{R}$ satisfies $u|_{V_{0}}=0,$ then we say that $u$ satisfies the Dirichlet boundary condition, and the eigenvalue problem

[TABLE]

is called the Dirichlet eigenvalue problem. A function $u$ satisfying both of these equations is called a Dirichlet eigenfunction of $\Delta.$ Similarly, we have a notion of a ‘Neumann condition’ as follows. Define the normal derivative $\partial_{n}(q_{i})$ for $i\in\{0,1,2\}$ by

[TABLE]

where we have identified indices modulo 3. If a function $u:SG\to\mathbb{R}$ satisfies $\partial u(q_{i})=0$ for $i=0,1$ and $2,$ then we say that $u$ satisfies the Neumann boundary condition, and the eigenvalue problem

[TABLE]

is called the Neumann eigenvalue problem. A function $u$ satisfying both of these equations is called a Neumann eigenfunction of $\Delta.$ Dirichlet and Neumann eigenfunctions on $SG$ are the analog of sine and cosine functions on the line.

5.3. Spectral Decimation

A method for explicitly computing all possible eigenvalues and eigenfunctions of $\Delta$ was introduced in [8] using a process called spectral decimation. Below we briefly discuss some results from spectral decimation we will use. Readers can find detailed discussion on spectral decimation in [8] and [22].

Proposition 5.2**.**

Suppose $\lambda_{m}\neq 2,5,6$, and $\lambda_{m-1}$ is given by

[TABLE]

(a) If $u$ is a $\lambda_{m-1}$-eigenfunction of $\Delta_{m-1}$ on $V_{m-1}$, then it can be uniquely extended to be a $\lambda_{m}$-eigenfunction of $\Delta_{m}$ defined on $V_{m}$. (b) Conversely, if $u$ is a $\lambda_{m}$-eigenfunction of $\Delta_{m}$ on $V_{m}$, then $u|_{V_{m-1}}$ is a $\lambda_{m-1}$-eigenfunction of $\Delta_{m-1}$ on $V_{m-1}$.

If we want to extend an eigenfunction of $\Delta_{m-1}$ with eigenvalue $\lambda_{m-1}$ to an eigenfunction of $\Delta_{m}$ using the proposition above, we have two choices, except when $\lambda_{m-1}=6$ (as we will see below), in which to extend the eigenfunction:

[TABLE]

For convenience, define the functions $\psi_{+}(x)$ and $\psi_{-}(x)$ by

[TABLE]

The numbers $2,5,6$ are called forbidden eigenvalues, and it turns out that each Dirichlet eigenfunction of $\Delta$ comes from a 2-, 5-, or 6-eigenfunction of $\Delta_{m_{0}}$ for some $m_{0}\geq 0$, while all the Neumann eigenfunctions come from 5- or 6-eigenfunctions of $\Delta_{m_{0}}$ for some $m_{0}\geq 0$. If $u$ is a Dirichlet or Neumann eigenfunction, we call $m_{0}$ the generation of birth and $u|_{V_{m_{0}}}$ the initial function. Suppose $u$ is an eigenfunction of $\Delta$ with eigenvalue $\lambda$ arising from initial eigenvalue $\lambda_{m_{0}}.$ Then we say that $\lambda$ is a 2-series (resp., 5-series or 6-series) eigenvalue if $\lambda_{m_{0}}=2$ (resp., $\lambda_{m_{0}}=5$ or 6). With a fixed generation of birth $m_{0}$ and initial function $f$, we can extend the function level-by-level. We can first fix a sequence $\{\varepsilon_{m}\}_{m=1}^{\infty}$ of $\pm$, with only finitely many $-$, and then let $\lambda_{m}=\psi_{\varepsilon_{m-m_{0}}}(\lambda_{m-1})$ for $m>m_{0}$ inductively. Then the function is extended to be an eigenfunction of $\Delta$, with the corresponding eigenvalue

[TABLE]

In fact, all the eigenfunctions with a given generation of birth and initial function can be generated by the above recipe. Also, if the initial eigenvalue is $6$, we can only choose $\varepsilon_{1}=+1$, as $\psi_{-}(6)=2$ is a forbidden eigenvalue. Now, suppose we fix a generation of birth $m_{0}$, fix an initial function on $V_{m_{0}}$, and let

[TABLE]

where $\emptyset$ denotes the empty sequence. Define $\psi_{e}=\psi_{e_{1}}\psi_{e_{2}}\cdots\psi_{e_{|e|}}$, where $|e|$ is the length of $e$. In particular, $\psi_{\emptyset}(x)=x$. Let $\Psi(x)=\frac{3}{2}\lim_{l\to\infty}5^{l}\psi^{l}_{-}(x)$. Then, we can deduce the following possibilities:

1. If $\lambda_{m_{0}}=2$ or $5$, then all the possible eigenvalues of $\Delta$ having generation of birth $m_{0}$ are given by

[TABLE]

2. If $\lambda_{m_{0}}=6$, then all the possible eigenvalues of $\Delta$ having generation of birth $m_{0}$ are given by

[TABLE]

Here we remark that when we use the notation $\lambda_{e}$, we always assume that we have a fixed generation of birth $m_{0}$ and initial eigenvalue $\lambda_{m_{0}}$. There is a method, given by [5], in which to arrange in increasing order the set of eigenvalues arising a fixed generation of birth and initial eigenvalue. The idea is to translate each finite sequence in $E$ into a binary number. The process is as follows. Given $e\in E$ of length $|e|=n$, let $d(e)$ be the integer with binary (base-2) expansion $d(e)=\sum_{i=1}^{n}2^{n-i}d_{i}(e)$, where

[TABLE]

In addition, we set $d(\emptyset)=0$. Then $\lambda_{e}$ is the \big{(}d(e)+1\big{)}-th smallest eigenvalue. For example, if $e=(e_{1},e_{2},e_{3},e_{4},e_{5},e_{6},e_{7})=(+,-,+,+,+,-,-)$, then $d(e)=d(+,-,+,+,+,-,-)=1101000_{2}=104,$ where $1101000_{2}$ is written in base-2. Thus, $\lambda_{e}$ corresponding to the sequence $e$ is the 105th smallest eigenvalue. Also, note that, in particular, $\lambda_{e=\emptyset}$ is the smallest eigenvalue, $\lambda_{e=(+)}$ is the 2nd smallest eigenvalue, and $\lambda_{e=(+,+)}$ is the 3rd smallest eigenvalue. We can, of course, reverse this process so that, given an in integer $d$ we can find the $e\in E$ corresponding to the $d$-th smallest eigenvalue. Another fact is that the sequence of eigenvalues $\{\lambda_{n,m_{0}}\}_{n=1}^{\infty}$, where $\lambda_{1,m_{0}}<\lambda_{2,m_{0}}<\lambda_{3,m_{0}}<\lambda_{4,m_{0}}<...$, corresponding to a fixed generation of birth $m_{0}$ and a fixed initial eigenfunction grow according to the power law $n^{\log 5/\log 2}$. In addition, if we define the eigenvalue counting function $\rho_{m_{0}}:\mathbb{R}_{+}\to\mathbb{N}$ to be

[TABLE]

then we have the following proposition concerning the asymptotic behavior of $\rho_{m_{0}}$.

Proposition 5.3**.**

For a fixed generation of birth $m_{0}$ and fixed initial function on $V_{m_{0}}$, there exists a $\log 5$-periodic continuous function $g(t)$ such that

[TABLE]

Proof. To avoid the high multiplicity of eigenfunctions corresponding to a same eigenvalue, we fix an initial eigenfunction instead of just fixing an initial eigenvalue. In our setting, the eigenfunction is unique for each eigenvalue, so we only need to count the number of eigenvalues. For convenience, we prove the proposition for generation of birth $m_{0}$ and initial eigenvalue $5$. For initial eigenvalue $2$ and $6$, the arguments are essentially the same. We will show that $\frac{\rho_{m_{0}}(5^{n}x)}{2^{n}x^{\log 2/\log 5}}$ converges to $g(\log x)$ uniformly on some interval $[e^{c},5e^{c})$ as $n\to\infty$. First, we have the following observations. Observation 1: For eigenvalues of generation of birth $m_{0}$ and initial eigenvalue $5$, if $n$ is fixed then we have

[TABLE]

for all $|e|>n$ where $\psi_{l}=\psi_{l_{1}}\circ\psi_{l_{2}}\cdots\psi_{l_{n}}$ for each word l=(l_{1},l_{2},\cdots,l_{n})\in\{-,+\}^{n}=\big{\{}l=\{l_{j}\}_{j=1}^{n}:l_{j}=+\text{ or }-\big{\}}. Notice that $5^{-{|e|}}\lambda_{e}\in 5^{m_{0}}\Psi\circ\psi_{e}([0,5])$, and the fact $\psi_{-}([0,5])\subset[0,5],\psi_{+}([0,5])\subset[0,5]$, we have $5^{-{|e|}}\lambda_{e}\in 5^{m_{0}}\Psi\circ\psi_{e_{1}e_{2}\cdots e_{n+1}}([0,5])$. The observation follows. For each $l\in\bigcup_{n=0}^{\infty}\{-,+\}^{n}$, we denote the endpoints of $5^{m_{0}}\Psi\psi_{+}\psi_{l}([0,5])$ by $a_{l}$ and $b_{l}$, i.e., $5^{m_{0}}\Psi\psi_{+}\psi_{l}([0,5])=[a_{l},b_{l}]$. By some easy computation, we can get the following observation. Observation 2: Assume $l=(l_{1},l_{2},\cdots,l_{n})$. Then

[TABLE]

To show observation 2, we need to consider two cases. First, for any $e$ with $|e|\leq k-1$, we have $\lambda_{e}<5^{k}b_{l}$, which counts for $2^{k-1}$ eigenvalues. Second, we have $2^{k-n-1}$ eigenvalues $\lambda_{e}$ in each interval of the form $5^{k}5^{m_{0}}\Psi\psi_{+}\psi_{l^{\prime}}([0,5])$, where $l^{\prime}\in\{+,-\}^{n}$. In fact, for $|e|=k$, $\lambda_{e}\in 5^{k}5^{m_{0}}\Psi\psi_{+}\psi_{l^{\prime}}([0,5])$ if and only if $(e_{1},e_{2},\cdots,e_{n+1})=(+,l^{\prime}_{1},l^{\prime}_{2},\cdots,l^{\prime}_{n})$, and we have $2$ free choices for each $e_{j},n+2\leq j\leq k$. There are $d(+l)+1-2^{n}$ intervals in $[0,b_{l}]$, noticing that $b_{l^{\prime}}\leq b_{l}$ if and only if $d(+l^{\prime})\leq d(+l)$. Combining the above facts, we get the second term 2^{k-n-1}\big{(}1+d(+l)-2^{n}\big{)}. Now, fix $n$ and consider $l\in\{-,+\}^{n}$. It is easy to see that $\frac{\rho_{m_{0}}(5^{k}b_{l})}{2^{k}b_{l}^{\log 2/\log 5}}$ converges as $k\to\infty$, since it is a constant for $k\geq n+1$. We denote the limit $g(\log b_{l})$. Next, we look at general $x$. Note that we can find some constant $c$ such that

[TABLE]

since $5^{m_{0}-1}\Psi\psi_{+}(5)<5^{m_{0}}\Psi\psi_{+}(0)$. We want to show that $\frac{\rho_{m_{0}}(5^{n}x)}{2^{n}x^{\log 2/\log 5}}$ converges to some function $g(\log x)$ uniformly on $c\leq\log x\leq c+\log 5$ as $n\to\infty$. In fact, if we fix $n$ and look at $b_{l}\leq x\leq b_{l^{\prime}}$ for some $l,l^{\prime}\in\{+,-\}^{n}$ such that $d(+l^{\prime})=d(+l)+1$, we have

[TABLE]

for any $k\geq n+1$. In addition, if $e^{c}\leq x\leq\min_{l\in\{+,-\}^{n}}a_{l}$, we have

[TABLE]

where $l=(-,-,-,\cdots)\cdots\in\{+,-\}^{n}$. Similarly, we have for $\max_{l\in\{+,-\}^{n}}b_{l}\leq x\leq 5e^{c}$,

[TABLE]

where $l=(+,-,-,\cdots)\in\{+,-\}^{n}$. The above discussions shows that $\frac{\rho_{m_{0}}(5^{n}x)}{2^{n}x^{\log 2/\log 5}}$ converges to some function $g(\log x)$ uniformly on $c\leq\log x\leq c+\log 5$ as $n\to\infty$. Also we can easily see that $g$ is continusous with $g(c)=g(\log 5)$ from the estimates. In fact, for each $x$, we can find a small neighbourhood such that for $n>k$ and any $y$ in the neighborhood, $\left|\frac{\rho(5^{k}y)}{2^{k}y^{\log 2/\log 5}}-\frac{\rho(5^{k}x)}{2^{k}x^{\log 2/\log 5}}\right|<2^{-n}x^{-\log 2/\log 5}$. The estimate obviously holds for the limit function. We can extend $g$ to be periodic on $\mathbb{R}$, and the theorem follows immediately. $\square$

5.4. Infinite Sierpinski Gasket

In the last part of this section, we introduce the infinite Sierpinski gasket ($SG_{\infty}$). It is a particular example of fractal blow-ups introduced in [24] by R. S. Strichartz. Recall that the Sierpinski gasket is defined by the self-similar identity, $SG=\bigcup_{i=0}^{2}F_{i}(SG)$, where each $F_{i}$ is a contraction mapping $\mathbb{R}^{2}\to\mathbb{R}^{2}$ of contraction ratio $\frac{1}{2}$ for $i=0,1,2$, as defined earlier in this section. The infinite Sierpinski gasket is constructed as follows.

Definition 5.4**.**

Suppose a sequence $\mathcal{K}=\{k_{n}\}_{n\geq 1}$, $k_{n}\in\{0,1,2\}$, is fixed. Define $SG_{M}=F^{-1}_{\mathcal{K},M}SG$, where $F_{\mathcal{K},M}=F_{k_{1}}F_{k_{2}}\cdots F_{k_{M}}$. Then the infinite Sierpinski Gasket $SG_{\infty}$ is defined by $SG_{\infty}=\cup_{M\geq 1}SG_{M}$. ∎

The Laplacian $\Delta_{\infty}$ on $SG_{\infty}$ can be defined locally with graph approximation in a same way as on $SG$. In [27], a Sierpinski lattice was introduced to describe the infinite graphs that approximate $SG_{\infty}$. Define

[TABLE]

and say $x\sim_{(m)}y$ if $F_{\mathcal{K},M}(x)\sim_{m+M}F_{\mathcal{K},M}(y)$ for some $M$. Then the resulting infinite graph is called a Sierpinski lattice. We can still define the discrete Laplacian on the lattices by

[TABLE]

Then the continuous Laplacian $\Delta$ is defined by

[TABLE]

One of the most important results on $SG_{\infty}$ was A. Teplyaev’s theorem (see [27]) below showing that the Laplacian $\Delta_{\infty}$ on $SG_{\infty}$ has pure point spectrum, which means the eigenfunctions of the Laplacian form a complete set.

Theorem 5.5** (A. Teplyaev).**

The Laplacian $\Delta_{\infty}$ is self-adjoint in $L^{2}(SG_{\infty},\mu)$, where $\mu$ is the Hausdoff measure on $SG_{\infty}$. The spectrum of $\Delta_{\infty}$ is pure point (i.e., the eigenfunctions of $\Delta_{\infty}$ form a basis of $L^{2}(SG_{\infty},\mu)$) and each eigenvalue has infinite multiplicity. The set of eigenfunctions with compact support is complete in $L^{2}(SG_{\infty},\mu)$. ∎

As a result of the theorem, spectral decimation still works on $SG_{\infty}$. Each of the eigenfunctions of $\Delta_{\infty}$ is an extension of an eigenfunction of $\Delta_{(m_{0})}$ with eigenvalue $5$ or $6$ by spectral decimation. The only difference here is that the generation of birth $m_{0}$ takes values in $\mathbb{Z}$ instead of $\mathbb{N}$. All the results concerning eigenvalues from a same generation of birth and initial function in the previous section, including Proposition 5.3, still hold on $SG_{\infty}$.

6. Extending the Mathieu Differential Equation to Infinite Fractafolds

Now we are ready to discuss how we will define the Mathieu differential equation on an infinite fractafold.

6.1. Defining the Fractal MDE

Recall that the MDE, defined on the real line, is given by

[TABLE]

where $u$ is a function from $\mathbb{R}$ to $\mathbb{R}$. The first questions we wish to address are “What should the fractal space be that replaces the line?” and “what should ‘periodic function’ mean?”. We choose to consider the infinite Sierpinski gasket $SG_{\infty}$ to be our domain. By “periodic function,” we mean a function on $SG_{\infty}$ which is identical on all the copies of $SG$ of the same size. In particular, if we are given a function $u$ on $SG$ with the boundary conditions

[TABLE]

we can get a periodic function on $SG_{\infty}$ by translating the function to other copies. But what about the differential equation? The first step in finding a fractal analog of the MDE defined on the line is to replace $\frac{d^{2}}{dt^{2}}$ with the fractal Laplacian $\Delta$, since $\Delta$ is the analog of the second-derivative operator. Now, what to do with the $\varepsilon\left(\cos x\right)u$ term? Recall from Theorem 5.1 above that the multiplication of two nonconstant functions in dom$\Delta$ may result in a function which is not in dom$\Delta$. Thus, we cannot simply replace $\cos x$ by a function in dom$\Delta$, and so we must figure out a suitable analog of multiplication by cosine. Recall that, in the line case, we sought solutions $u$ in the form of a Fourier expansion in terms of cosines and sines. Note, however, that cosines and sines on the line are Neumann and Dirichlet eigenfunctions, respectively, of $\frac{d^{2}}{dt^{2}}$. Hence, we will adopt a form of Mathieu’s equation which is compatible with functions $u$ that can be written as a linear combination of Neumann eigenfunctions, motivated by Equation (6.1). We choose to only consider functions which have Neumann eigenfunction expansions of the form

[TABLE]

where each $\varphi_{j}$ is a Neumann eigenfunction function defined as follows. Fix a generation of birth, a series (5-series or 6-series), and an initial eigenfunction $\varphi$ such that $\Delta_{m_{0}}\varphi(x)=\lambda_{m_{0}}\varphi(x)$ for all $x\in V_{m_{0}}$ and $\varphi(q_{0})=\varphi(q_{1})=\varphi(q_{2})$. With the spectral decimation algorithm introduced in Section 5.3, we obtain a set of Neumann eigenfunctions $\varphi_{i}$ of $\Delta$ extended from $\varphi$, and we write $\lambda_{i}$ for the eigenvalue corresponding to $\varphi_{i}$. We still take the order $\lambda_{1}<\lambda_{2}<\lambda_{3}<\cdots$ as in Section 5.3. Readers can also find more details on Neumann eigenfunctions on $SG$ in [22]. The reason we fix a common initial Neumann eigenfunction is that, if we do not fix such an initial eigenfunction and instead consider the set of all Neumann eigenfunctions of $\Delta$, then we cannot order their eigenvalues in a discrete way as above. To this end, suppose we have a function $u$ on $SG$ which can be written as a linear combination of Neumann eigenfunctions as in Equation 6.2, where the $\varphi_{j}$ ($j\geq 0$) satisfy the definition in the previous paragraph. Then, for any $\varphi_{j}$ with $j\geq 2$ we define multiplication by cosine as follows:

[TABLE]

The motivation for this definition comes from the fact that, for the line case, the Neumann eigenfunction $\cos(jt)$ obeys the following trigonometric property when multiplied by $\cos t$:

[TABLE]

As for $j\geq 1$ we consider two possibilities, each of which will be described in Section 6.2. In addition, we will consider another two variant versions in the following subsection. Using our developments thus far, the fractal Mathieu differential equation would say

[TABLE]

where $\lambda_{j}$ is the eigenvalue of $\Delta$ corresponding to eigenfunction $\varphi_{j}$. Since the set $\{\varphi_{j}\}$ of eigenfunctions is linearly independent, we must have

[TABLE]

Putting these equations into matrix form we obtain

[TABLE]

Note that this matrix takes the form of Equation 3.1, which reproduce below:

[TABLE]

6.2. Variants of the Mathieu Differential Equation on the Line

We will consider 4 different ‘versions’, $M_{1}$, $M_{2}$, $M_{3}$, and $M_{4}$, of the coefficient matrix for the fractal MDE $M_{i}x=0$ ($i=1,2,3,4$):

• •

Version 1:

[TABLE]

This matrix is reminiscent of the cosine matrices for the line case. Note that the first term in the second row is $\varepsilon$, not $\frac{1}{2}\varepsilon$. The recursion relation for the coefficients $c_{j}$ becomes

[TABLE]

• •

Version 2:

[TABLE]

This matrix is reminiscent of the sine matrices for the line case. Note that the eigenvalues start from $\lambda_{1}$ instead of $\lambda_{0}$. The recursion relation for the coefficients $c_{j}$ becomes

[TABLE]

• •

Versions 3 and 4: For Version 3 and for Version 4 we take the following approach. We now consider a variant of the Mathieu differential equation, given by

[TABLE]

where $A$ is the operator analogous to multiplication by cosine. For $j\geq 2$, define $A\cdot\varphi_{j}:=\alpha_{j-1}\varphi_{j-1}+\beta_{j}\varphi_{j+1},$ with $\alpha_{j}$ and $\beta_{j}$ satisfying

[TABLE]

for $j\geq 1$. One can solve system 6.7 to find the solutions of $\alpha_{j}$ and $\beta_{j}$

[TABLE]

Note that $0<\alpha_{j}<1$ and $0<\beta_{j}<1$ for all $j$, and hence the sequences $\{\alpha_{j}\}$ and $\{\beta_{j}\}$ are each uniformly bounded. The motivation for the setup above is that, if we plug in $\lambda_{j}=j^{2}$, which corresponds to the eigenvalues on the line, Equation 6.8 yields $\alpha_{j}=\beta_{j}=\frac{1}{2}$. For $j=0$ and $j=1$, we still consider two cases, which give us Version 3 and Version 4 as follows. We let Version 3 be

[TABLE]

where $\alpha_{j}$ and $\beta_{j}$ are as given in Equation 6.8. This matrix is reminiscent of the cosine matrices for the line case. Note that the first term in the second row has an extra factor of $2$. The recursion relation for the coefficients $c_{j}$ becomes

[TABLE]

We let Version 4 be

[TABLE]

where $\alpha_{j}$ and $\beta_{j}$ are as given in Equation 6.8. This matrix is reminiscent of the sine matrices for the line case. Again, note that the eigenvalues start from $\lambda_{1}$ instead of $\lambda_{0}$. The recursion relation for the coefficients $c_{j}$ becomes

[TABLE]

7. Ovservations and Analysis for the Fractal Mathieu Differential Equation

In this section we parallel our results presented in Section 4 by giving a discussion of the asymptotic behavior of the transition curves (defined below) for the fractal MDE and of the convergence of solutions. We also describe a phenomenon wherein the $SG_{\infty}$ transition curves form a prominent ‘diamond’ pattern.

7.1. The $\delta$-$\varepsilon$ Plot

In Section 6, we introduced 4 different Versions $M_{i}x=0$, ($i=1,2,3,4$) of the fractal MDE in matrix form on $SG_{\infty}$. If one fixes one of the four Versions of the MDE as above, then the points $(\delta,\varepsilon)$ in the $\delta$-$\varepsilon$ plane for which the $M_{i}x=0$ has a nontrivial solution in $\ell^{2}$ make up the transition curves for that Version of the MDE. In Figures 7.1-7.8, we show the transition curves for each of the four Versions.

As can be viewed from Figures 7.1-7.8, in each plot there is a ‘diamond’ pattern formed between adjacent transition curves. That is, the curves seem to ‘bounce’ off each other before parting in separate directions. This pattern does not appear in the transition curves for the Mathieu differential equation on the line. This is one way in which the transition curves for $SG_{\infty}$ are different than those for the line. We give a short proof here concerning the asymptotic behavior of the transition curves for Versions 1 and 2.

Proposition 7.1**.**

For version 1 and version 2, fix one transition curve. The ratio $\frac{\delta}{\varepsilon}$ converges to $-1$ if $\varepsilon\to\infty$. The ratio $\frac{\delta}{\varepsilon}$ converge to $1$ if $\varepsilon\to-\infty$.

Proof. In fact, $M_{1}$ and $M_{2}$ are of the form $A+\varepsilon\cdot B$, where $A$ is a diagonal matrix, and $B$ has one of the following two forms,

[TABLE]

The first matrix is real symmetric with continuous spectrum $[-1,1]$. To see this, we only need observe that $B=U\cdot\cos t\cdot U^{-1}$, where $U:\ell^{2}\to L^{2}([0,\pi])$ is defined as $U\bm{c}=\sum_{n=1}^{\infty}c_{n}\sin nt$. The second one is not symmetric, but we can study the matrix $\kappa^{-1}(A-\varepsilon\dot{B})\kappa=A-\kappa^{-1}B\kappa$ where $\kappa$ is a diagonal matrix, we did in the proof for theorem 4. After the transformation, we have

[TABLE]

which is also a real symmetric matrix with continuous spectrum $[-1,1]$. Without loss of generality, we just assume $B$ is symmetric in the follows. In the following proof, we only consider $\varepsilon\to+\infty$ case, and the other half is then immediate by symmetry. On the $k$-th curve, we know that

[TABLE]

where the infimum takes over all the $k$ dimensional subspace $M_{k}$ of dom$A$. It is easy to see the following lower bound of $\delta$

[TABLE]

Next, we give an upper bound. For any $l<1$, we can find a $k$ dimensional subspace $N_{k}\subset domA$, such that $-1<B|_{N_{k}}<l$, which means

[TABLE]

In fact, observe that $B$ has continuous spectrum $[-1,1]$, we can always find a $k$ dimensional subspace $\tilde{N}_{k}\subset\ell^{2}$ such that $-1<B|_{\tilde{N}_{k}}<l$. But $\tilde{N}_{k}$ is not necessarily in $domA$. So we pick a large $L$, and define

[TABLE]

Obviously, $N_{k}\subset domA$, and $-1<B|_{N_{k}}<l$ as long as $L$ is large enough. Then we get the following upper bound

[TABLE]

Combining the upper and lower bounds of $\delta$, we get

[TABLE]

for any $l<1$. Thus $\lim_{\varepsilon\to+\infty}\frac{\delta}{\varepsilon}=-1$. By symmetry, we also get $\lim_{\varepsilon\to-\infty}\frac{\delta}{\varepsilon}=1$. $\square$

7.2. Solutions of the Mathieu Differential Equation

Now we investigate solutions to the fractal MDE. From Corollary 3.3 we deduce that, for $\varepsilon\neq 0,$ a nontrivial solution $\textbf{c}=(c_{1},c_{2},...)\in\ell^{2}$ (or $\bm{c}=(c_{0},c_{1},c_{2},...)\in\ell^{2}$) to any of the four systems $M_{i}x=0$ $(i=1,2,3,4)$ in section 6.2 will decay with asymptotic behavior $\frac{c_{j}}{c_{j-1}}\sim\frac{\beta_{j-1}\varepsilon}{\lambda_{j}-\delta}.$ Hence, a solution

[TABLE]

to any of the four versions of the fractal MDE may be well-approximated by the sum of just the first few terms in the series, and such an approximation will yield only a small error to the actual infinite-series solution. Recall that in Section 4 we investigated, for the MDE on the line, how properties of solutions change if one fixes a transition curve and considers periodic solutions corresponding to various $(\delta,\varepsilon)$ points along that transition curve. We will consider the same question for the $SG_{\infty}$ transition curves for the fractal MDE as well. We investigate this question by first studying the solutions plotted in Figures 7.10-7.29. In constructing the solutions we choose eigenfunctions of generation of birth 2 and of initial eigenvalue $5$ or $6$, arising from either of the initial functions shown in Figure 7.9, Figures 7.10-7.15 correspond to solutions of Version 1, Figures 7.16-7.21 correspond to solutions of Version 2, Figures 7.22-7.25 correspond to solutions of Version 3, and Figures 7.26-7.29 correspond to solutions of Version 4. Each Figure includes two images. The one on the left is a plot of three different solutions corresponding to different $(\delta,\varepsilon)$ pairs along the first transition curve; the set of three solutions in each Figure corresponds to either 5-series or 6-series, and each solution been normalized so that the maximum value it attains is 1. The image on the right in each Figure shows an overlook of the three solutions, where the color at each point indicates which solution takes on the biggest value among the three. From the plots, it is evident that these normalized solutions seem to converge as $\varepsilon$ tends to infinity along the first transition curve. The second way in which we investigate the question is by observing the behavior of the location of relative maxima of solutions as $\varepsilon$ traverses along a transition curve. Here we still choose to study $(\delta,\varepsilon)$ pairs and corresponding solutions on the first transition curve for each version. See Figure 7.30-7.37. This is analogous to our study of the relative maxima of solutions to the MDE on the line in Section 4.2. Each of the Figures contains two images indicating the position of relative maxima: the image on the left shows how the positions of the maximas move as $\varepsilon$ becomes larger, and the image on the right is an overlook of the left one where we can clearly see the locations of the peaks. We make several observations:

• •

The first thing we can see is that the movement of the peaks are not large, and we do not observe the peaks converging to the boundary.

• •

The second observation is that the $5$ and $6$ series behave quite differently. We can observe jumps of peaks in all the versions for $5$ series, but the number of peaks is usually 3. A special case is figure 7.36, where many peaks occur when $\varepsilon$ is quite large. However, for $6$ series, we can observe many more peaks when $\varepsilon$ is very small, but most of them do not appear when $\varepsilon$ is large. Until now, we do not have explanation for these behaviors.

8. Further Research

Now we discuss further research questions that can be investigated in this area.

1. (i)

Asymptotic behavior of $SG_{\infty}$ transition curves One could further investigate the asymptotic behavior of the $SG_{\infty}$ transition curves. For the line case, we have Theorem 4.2 due to W.S. Loud, stating that

[TABLE]

as $\varepsilon\to+\infty$ on the $k-$th transition curve. We can investigate whether estimates of a similar form hold on $SG_{\infty}$. This would include extending Proposition 7.1 to Version 3 and Version 4 the MDE on $SG_{\infty}$. 2. (ii)

Other modifications of MDE matrix Further research could investigate different matrix versions of the fractal MDE, aside from Versions 1-4 presented in Section 6. 3. (iii)

Other Fractal Domains One could investigate how the Mathieu differential equation could be extended to other infinite fractafolds. Such infinite fractafolds may or may not be based on $SG$. 4. (iv)

The Hill Equation The Mathieu differential equation is actually a special case of the so-called Hill differential equation. The Hill differential equation is given by

[TABLE]

where $f(t)$ is an arbitrary periodic function. Readers can read [21, 16, 17, 28] for details. Further research could investigate ways in which to extend, for various choices of $f(x)$, the Hill differential equation to be defined on fractal domains.

9. Appendix

In this appendix we derive equations in matrix form for solutions to the MDE on the line which have period $2N\pi$, $N\in\mathbb{N}$. We will use a similar procedure as in Section 2 and will obtain tridiagonal matrices. Let $u$ be a solution with period $2N\pi$. We can write the Fourier expansion of $u$ as

[TABLE]

Plugging in this Fourier series for $u$ into the Mathieu differential equation, we obtain two infinite systems of linear homogeneous equations for the cosine and sine coefficients, respectively, as follows:

[TABLE]

and

[TABLE]

We can immediately write the equations in the matrix form, but the result is unwieldy. So, we make further classifications of the coefficients, just as we did in Section 2.2, where coefficients are separated into two classes—those with even indices and those with and odd indexes. A natural criterion is to have coefficients that appear in the same equation be in a same class. For example, $a_{0}$ and $a_{N}$ should be in a same class. Based on this idea, we say that $a_{j}$ and $a_{j^{\prime}}$ (or $b_{j}$ and $b_{j^{\prime}}$) are in the ‘same class’ if and only if there is a finite sequence

[TABLE]

such that $a_{j_{l}}$ and $a_{j_{l+1}}$ appear in a same equation for all $0\leq l\leq L-1$. For example, $a_{0}$ and $a_{3N}$ are in the same class, since we $a_{0}$ and $a_{N}$ are in the same equation, $a_{N}$ and $a_{2N}$ are in the same equation, and $a_{2N}$ and $a_{3N}$ are in the same equation. It is easy to check that the the property of being in the ‘same class’ is an equivalence relation on the set $\{a_{j}:j\geq 0\}$ and on the set $\{b_{j}:j\geq 1\}$. This partitions $\{a_{j}:j\geq 0\}$ into $[\frac{N}{2}]+1$ different equivalence classes and partitions $\{b_{j}:j\geq 0\}$ into $[\frac{N}{2}]+1$ different equivalence classes. So we can write equations for different equivalence classes separately. Below, we discuss all the possible cases. First we give a discussion for the cosine coefficients $\{a_{j}\}$. As we will see, the case for the sine coefficients $\{b_{j}\}$ is similar. There are three possible forms that the matrix corresponding to any particular equivalence class can take: 1). The first possible equivalence class is $\{a_{0},a_{N},a_{2N},\cdots\}$, with the corresponding matrix equation

[TABLE]

Note that this is just the equation $Ax=0$ in Section 2.2. If there is a nontrivial solution to the above equation, then the MDE has a $2\pi$ periodic solution, since $u(t)=\sum_{l=0}^{\infty}a_{lN}\cos(lt)$ solves the MDE. 2). The second form that an equivalence class can take is

[TABLE]

where $1\leq k<\frac{N}{2}$. The corresponding equations can be written in the form

[TABLE]

If there is a nontrivial solution to the above equation, then $u(t)=\sum_{l=-\infty}^{\infty}a_{k+lN}\cos(\frac{k+lN}{N}t)$ is a $\frac{2N\pi}{\text{gcd}(N,K)}$-periodic solution to the MDE, where $\text{gcd}(N,k)$ denotes the greatest common divisor of the pair $(N,k)$. In particular, we get a $2N\pi$-periodic solution to the MDE from the equation above if and only if $k$ and $N$ are coprime. Thus, the matrices of the second form with $k$ and $N$ coprime yield $2N\pi$-periodic solutions. 3) The third form that an equivalence class can take is $\{a_{\frac{N}{2}},a_{\frac{3}{2}N},a_{\frac{5}{2}N},\cdots\}.$ This form can only occur if $N$ is even. The corresponding equations can be written in the form

[TABLE]

The above matrix is exactly the matrix $C$ in Section 2.2, which yield $4\pi$-periodic solutions.

The matrix forms for sine coefficients are quite similar to those for the cosine coefficients. For solutions of period $2N\pi,N\geq 3$, the equations for sine coefficents have the following form

[TABLE]

where $k$ can by any integer such that $k$ and $N$ are coprime.

Acknowledgements

This research was hosted by the Cornell University Department of Mathematics through its 2018 Summer Program for Undergraduate Research. Anthony Coniglio’s participation in this research was partly supported by Indiana University Bloomington. Xueyan Niu’s participation in this research was partly supported by the Overseas Research Fellowship (ORF) of Faculty of Science, University of Hong Kong.