# A New Approach for Linear Eigenvalue Problems and Nonlinear Euler Buckling Problem

Authors: Meltem Evrenosoglu Adiyaman; Sennur Somali
Journal: Abstract and Applied Analysis
 (2012)

---

## Abstract

We propose a numerical Taylor’s Decomposition method to compute approximate eigenvalues and eigenfunctions for regular Sturm-Liouville eigenvalue problem and nonlinear Euler buckling problem very accurately for relatively large step sizes. For regular Sturm-Liouville problem, the technique is illustrated with three examples and the numerical results show that the approximate eigenvalues are obtained with high-order accuracy without using any correction, and they are compared with the results of other methods. The numerical results of Euler Buckling problem are compared with theoretical aspects, and it is seen that they agree with each other.


---

## Body

## 1. Introduction

We investigate the computation of eigenvalues of regular Sturm-Liouville eigenvalue problems:(1.1)-y′′(x)+r(x)y(x)=λy(x),0≤x0<x<xn,y(x0)=y(xn)=0,
where r(x)∈Cp+q[x0,xn] and p,q∈ℕ and Euler Buckling problem:
(1.2)y′′+λsin⁡y=0,y′(0)=0,y′(1)=0.Regular Sturm-Liouville problems arise in many applications, and many methods are available for their numerical solution Pryce [1].We also examine an elementary, classical problem buckling of an end-loaded rod which possesses a completely soluble continuous model in the form of a nonlinear, second-order boundary value problem as described in elsewhere [2–5]. An essential complete analysis of this problem was provided by Euler [6]. For the nonlinear eigenvalue problem (1.2), one may find that for small λ the only solution is zero solution as in the linear case. But as the eigenvalue λ increases, it reaches a critical value λ1 at which a nonzero solution appears, corresponding to buckling of the rod. For λ>λ1, the nonlinear problem behaves quite differently from the linear problem: for a range of values λ1<λ<λ2, there is exactly one nonzero solution of (1.2) for each λ, and when λ exceeds λ2, a second nonzero solution appears; similarly, there is a value λ3 beyond which there are three nonzero solutions, and so on. Namely, one may give inductively
(1.3)0≤λ≤π2,onlythetrivialsolution,π2<λ≤4π2,onenontrivialsolution,n2π2<λ≤(n+1)2π2,nnontrivialsolutions,
as given by Stakgold [2]. This behavior is a simple example of the phenomenon of bifurcation or branching; it occurs in many different areas of applied mathematics.The method considered here is a Taylor decomposition which was used by Adiyaman and Somali [7] for the solution of certain nonlinear problems. Like classical finite-difference and finite-element methods, this high order method is best suited to the fundamental eigenvalue and small eigenvalues.In Section2, the behavior of eigenvalues and corresponding eigenfunctions for regular Sturm-Liouville problem is obtained by Taylor’s decomposition method, and convergence of the method for regular Sturm-Liouville problem with constant function r(x) is given. We establish a Lemma and a Theorem, and then we give an application of Taylor’s decomposition method to the Euler Buckling problem in Section 3. The technique is illustrated with three examples, and the numerical results of regular Sturm-Liouville problem are given by comparing the results of other methods in Section 4. The numerical results of Euler Buckling problem accompanying the theoretical results and the behavior of solution are also discussed in Section 4. In the conclusion, we summarize the study and present our suggestions regarding future work.


## 2. Application and Error Analysis of Taylor’s Decomposition Method for Regular Sturm-Liouville Eigenvalue Problems

### 2.1. Application of Taylor’s Decomposition on Two Points for Regular Sturm-Liouville Eigenvalue Problems

We consider the regular Sturm-Liouville eigenvalue problem (1.1) by introducing the new depending variable y′(x)=z(x), (1.1) can be written as
(2.1)Y′(x)=A(x)Y(x),C0Y(x0)+C1Y(xn)=0,
where
(2.2)Y(x)=[y(x)z(x)],A(x)=[01r(x)-λ0],C0=[1000],C1=[0010].
From Ashyralyev and Sobolevskii [8], we will consider the application of Taylor’s decomposition of function Y(x) on two points. We need to find Y(j)(x) for any 1≤j≤pandq. Using the equation Y′(x)=A(x)Y(x), we get
(2.3)Y(j)(x)=Aj(x)Y(x),
with
(2.4)A0(x)=I,A1(x)=A(x),Aj(x)=Aj-1′(x)+Aj-1(x)A(x),2≤j≤p,
where I is the 2×2 identity matrix. By using the structure of the matrix A(x), we obtain the entries of the matrix of
(2.5)Aj(x)=[aj(1,1)(λ;x)aj(1,2)(λ;x)aj(2,1)(λ;x)aj(2,2)(λ;x)]
as in the following formulas:
(2.6)aj(1,1)(λ;x)=∂aj-1(1,1)(λ;x)∂x+(r(x)-λ)aj-2(2,2)(λ;x)=aj-1(2,1)(λ;x),aj(2,2)(λ;x)=∂aj-1(2,2)(λ;x)∂x+aj(1,1)(λ;x),aj(1,2)(λ;x)=aj-1(2,2)(λ;x),aj(2,1)(λ;x)=-∂aj(2,2)(λ;x)∂x+aj+1(2,2)(λ;x),
for 2≤j≤p, where
(2.7)a0(1,1)(λ;x)=1,a1(1,1)(λ;x)=0,a0(1,2)(λ;x)=0,a1(1,2)(λ;x)=1,a0(2,1)(λ;x)=0,a1(2,1)(λ;x)=r(x)-λ,a0(2,2)(λ;x)=1,a1(2,2)(λ;x)=0.From the theorem given in Ashyralyev and Sobolevskii [8], we have the following relation:
(2.8)Y(xk)-Y(xk-1)+∑j=1pαjY(j)(xk)hj-∑j=1qβjY(j)(xk-1)hj=(-1)p(p+q)!∫xk-1xk(xk-s)q(s-xk-1)pY(p+q+1)(s)ds
on the uniform grid
(2.9)[x0,xn]h={xk=x0+kh,k=0,1,…,n,nh=xn-x0,n∈N},
where
(2.10)αj=(p+q-j)!p!(-1)j(p+q)!j!(p-j)!,1≤j≤p,βj=(p+q-j)!q!(p+q)!j!(q-j)!,1≤j≤q.
Rewriting (2.8) by neglecting the last term, we obtain the single-step difference scheme of (p+q)-order of accuracy for the approximate solution of problem (2.1):
(2.11)Yk-Yk-1+∑j=1pαjAj(xk)Ykhj-∑j=1qβjAj(xk-1)Yk-1hj=0,
where
(2.12)Yk=[ykzk]
is the approximate value of Y(xk). For the simple computation, let p=q, then we have
(2.13)(I+∑j=1pαjAj(xk)hj)Yk=(I+∑j=1p(-1)jαjAj(xk-1)hj)Yk-1,
where βj=(-1)jαj. Letting M(xk)=(I+∑j=1pαjAj(xk)hj) and N(xk-1)=(I+∑j=1p(-1)jαjAj(xk-1)hj), we write
(2.14)Yk=M-1(xk)N(xk-1)Yk-1.
Since the accuracy and convergence of the method not only depend on h, they also depend on p, we can increase the order of accuracy by increasing p for fixed h. So h is chosen as length of the whole interval as follows. Now, taking h=xn-x0 gives
(2.15)Y1=M-1(xn)N(x0)Y0,
and substituting into the boundary condition of (2.1), we get
(2.16)(C1M-1(xn)N(x0)+C0)Y0=0.
To obtain a nontrivial solution Y0, we must have the following equation:
(2.17)det⁡(C1M-1(xn)N(x0)+C0)=0.
Defining
(2.18)M(xn)=[m11m12m21m22],N(x0)=[n11n12n21n22],
we have the following statement
(2.19)m22n12-m12n22=0.
Since
(2.20)M(xn)=[1+∑j=1pαjaj(1,1)(λ;xn)hj∑j=1pαjaj(1,2)(λ;xn)hj∑j=1pαjaj(2,1)(λ;xn)hj1+∑j=1pαjaj(2,2)(λ;xn)hj],N(x0)=[1+∑j=1p(-1)jαjaj(1,1)(λ;x0)hj∑j=1p(-1)jαjaj(1,2)(λ;x0)hj∑j=1p(-1)jαjaj(2,1)(λ;x0)hj1+∑j=1p(-1)jαjaj(2,2)(λ;x0)hj],
using the entries m12,m22,n12, and n22 of the above matrices and the properties of the entries of Aj(x), we obtain (2.19) in terms of λ:
(2.21)F(λ)=m22n12-m12n22=(1+∑j=1pαjaj(2,2)(λ;xn)hj)(∑j=1p(-1)jαjaj-1(2,2)(λ;x0)hj)-(∑j=1pαjaj-1(2,2)(λ;xn)hj)(1+∑j=1p(-1)jαjaj(2,2)(λ;x0)hj).
Solving nonlinear equation F(λ)=0 by Newton’s method, we find the approximate eigenvalues. This method appears to require a separate calculation for the eigenfunctions.To find the corresponding eigenfunctions of the regular Sturm-Liouville eigenvalue problem (2.1), we substitute the eigenvalue to (2.1) and we solve the obtained boundary value problem by Taylor’s decomposition method on two points xk-1 and xk with the uniform grid [0,1]h for p=q. Then, we get a homogeneous linear equation system of 2n equations with 2n unknown z0,y1,z1,y2,z2…,yn-1,zn-1,zn which are the approximated values of y′(x0),y(x1),y′(x1),y(x2),y′(x2),…,y(xn-1),y′(xn-1),y′(xn), respectively. Solving the 2n×2n homogeneous system, we obtain approximate values of the eigenfunction and its derivative of (1.1) at the point x=xk.


### 2.2. Error Analysis for Regular Sturm-Liouville Problem Whenr(x)=c

In this section, we will show the convergence of the method for eigenfunctions with the constant functionr(x)=c by obtaining approximate value of eigenfunction at the point x¯∈[x0,xn] of the problem (1.1). Without loss of generality, we may choose r(x)=0, then Aj(x)=Aj, that is, aj(2,2)(λ;xn)=aj(2,2)(λ;0)=aj(2,2)(λ). Using (2.6), we can find explicit values of aj(1,1), aj(2,2) as follows:
(2.22)a2j(1,1)=(-1)jλj,a2j(2,2)=(-1)jλj,a2j+1(1,1)=0,a2j+1(2,2)=0,j≥0.
This yields
(2.23)m22=1+∑j=1⌊p/2⌋α2j(-1)jλjh2j,n22=m22,m12=∑j=0⌊(p-1)/2⌋α2j+1(-1)jλjh2j+1,n12=-m12,m11=m22,m21=-λm12.
Using (2.14) for k=1, we have
(2.24)Y1=M-1(x¯)N(x0)Y0,
where Y0 and Y1 are the approximated values of Y(x0) and Y(x¯), respectively, with the stepsize h=x¯-x0:(2.25)Y1=z0det⁡(M)[-2m22m12-λ(m12)2+(m22)2].
The first component of the above vector (2.25) gives the approximate eigenfunction y1, and the second component of the above vector (2.25) gives the derivative of the approximate eigenfunction z1 of the regular Sturm-Liouville problem (1.1) at x¯. Now, we will show that y1 and z1 converge to exact functions y(x¯) and y′(x¯), respectively, as p→∞.Using the Stirling’s Formulan!≈2πn(n+1)/2e-n for αj in (2.10), we obtain
(2.26)αj=(2p-j)!p!(-1)j(2p)!j!(p-j)!≈(-1)j1j!12j(p-j/2p-j)(p-j+1)/2(p-j/2p)p/2.
This gives
(2.27)lim⁡p→∞αj=(-1)j1j!12j.
Thus,
(2.28)lim⁡p→∞m22=lim⁡p→∞(1+∑j=1⌊p/2⌋α2j(-1)jλjh2j)=1+∑j=1∞1(2j)!122j(-1)jλjh2j=∑j=0∞(-1)j(λh2)2j1(2j)!=cos⁡(λ)h2.
By using the same idea, we obtain
(2.29)lim⁡p→∞m12=lim⁡p→∞(∑j=0⌊(p-1)/2⌋(-1)jα2j+1λjh2j+1)=1λsin⁡(λ)h2.
It follows from (2.28) and (2.29) that
(2.30)lim⁡p→∞det⁡(M)=m222+λm122=cos⁡2(λ)h2+λ(1λsin⁡(λ)h2)2=1.
Hence, for r(x)=0, the approximate eigenfunction of (1.1) to the corresponding eigenvalue λ converges to exact eigenfunction:
(2.31)lim⁡p→∞y1=2z0det⁡(M)1λ(cos⁡(λh2))(sin⁡(λh2))=z0λsin⁡(λ(x¯-x0)).
Since we have z(x)=y′(x), the derivative of approximate eigenfunction of (1.1) to the corresponding eigenvalue λ converges to derivative of the exact solution:
(2.32)lim⁡p→∞z1=z0det⁡(M)((-λ)1λsin⁡2(λh2)+cos⁡2(λh2))=z0cos⁡(λ(x¯-x0)),
where λ=k2π2,k=1,2,….This demonstration shows that approximate eigenfunction and eigenvalue converges to exact one asp→∞ for fixed step-size “h.”


### 2.3. Taylor’s Decomposition Method to the Euler Buckling Problem

For convenience, we introduce the following notations as in (2.1) and Adiyaman and Somali [7]:
(2.33)Y(x)=[y(x)z(x)],F(Y(x))=[f1(0)(y,z)f2(0)(y,z)],C0=[0100],C1=[0001],f1(0)(y,z)=z,f2(0)(y,z)=-λsin⁡y.
Thus, the Euler Buckling Problem (1.2) can be written in the form:
(2.34)Y′(x)=F(Y(x)),0<x<1,C0Y(0)+C1Y(1)=0,
Defining the following recurrence relations for j=1,…,2p:
(2.35)fi(j)(y,z)=z∂fi(j-1)(y,z)∂y-λsin⁡y∂fi(j-1)(y,z)∂z,i=1,2,
we obtain
(2.36)Y(j)(x)=[f1(j-1)(y,z)f2(j-1)(y,z)]=[f2(j-2)(y,z)f2(j-1)(y,z)]forj=2,…,2p+1.
We first give the following lemma which defines f2(j-1)(y,z) explicitly.Lemma 2.1.
Forj=0,…,2p, let f2(j)(y,z) satisfy the recurrence relation (2.35) with f2(0)(y,z)=-λsiny. Then
(2.37)f2(2m)(y,z)=∑i=0mλi+1z2m-2i∑k=0⌊i/2⌋(-1)m+1-ka2m+1,i,k(cos⁡y)i-2k(sin⁡y)2k+1,
where
(2.38)a2m+1,i,k={1,i=0,k=0,(2k+2)a2m,i,k+1+(i+1-2k)a2m,i,k+(2m+1-2i)a2m,i-1,k,1≤i≤m,0≤k≤i2,0,else,
for m=0,…,p, and
(2.39)f2(2m+1)(y,z)=∑i=0mλi+1z2m+1-2i∑k=0⌊(i+1)/2⌋(-1)m+1-ka2m+2,i,k(cos⁡y)i+1-2k(sin⁡y)2k,
where
(2.40)a2m+2,i,k={1,i=0,k=0,(2k+1)a2m+1,i,k+(i-2k)a2m+1,i,k-1+(2m+2-2i)a2m+1,i-1,k-1,1≤i≤m,0≤k≤i20,else,
for m=0,…,p-1.Proof.
The proof follows induction argument based on (2.35).Theorem 2.2.
Iff1(j)(y,z) and  f2(j)(y,z), are sufficiently smooth and satisfy (2.36), (2.37), and (2.39) then the following relations hold: (a) it holds that
(2.41)y1-y0+∑j=1pβj[(-1)jf1(j-1)(y1,0)-f1(j-1)(y0,0)]=-2(y0+∑j=1pβjf1(j-1)(y0,0)),∑j=1pβj[(-1)jf2(j-1)(y1,0)-f2(j-1)(y0,0)]=0,
for y1=-y0, (b) it holds that
(2.42)∑j=1pβj[(-1)jf2(j-1)(y1,0)-f2(j-1)(y0,0)]=-2∑j=1pβjf2(j-1)(y0,0),y1-y0+∑j=1pβj[(-1)jf1(j-1)(y1,0)-f1(j-1)(y0,0)]=0,
for y1=y0.Proof.
Letj=2m+1 for m=0,…,p, then f1(2m+1)(y,z) becomes
(2.43)f1(2m)(y,z)=f2(2m-1)(y,z)=∑i=0(m-1)λi+1z2m-2i-1∑k=0⌊(i+1)/2⌋(-1)m-ka2m,i,k(cos⁡y)i+1-2k(sin⁡y)2k,
by Lemma 2.1. Since all terms of previous sum contain z, f12m(y,0)=f22m-1(y,0)=0 for m=1,…,p, hence, we get the following equations:
(2.44)y1-y0+∑j=1pβj[(-1)jf1(j-1)(y1,0)-f1(j-1)(y0,0)]=y1-y0+∑j=1pβj[f1(j-1)(y1,0)-f1(j-1)(y0,0)],∑j=1pβj[(-1)jf2(j-1)(y1,0)-f2(j-1)(y0,0)]=∑j=1pβj[-f2(j-1)(y1,0)-f2(j-1)(y0,0)].
Letting j=2m+2 for m=0,…,p-1 and using (2.37), we get
(2.45)f1(2m+2-1)(y,z)=f2(2m)(y,z)=∑i=0mλi+1z2m-2i∑k=0⌊i/2⌋(-1)m+1-ka2m+1,i,k(cos⁡y)i-2k(sin⁡y)2k+1.
Substituting the value z=0 into (2.45), we obtain
(2.46)f1(2m+1)(-y0,0)=f22m(-y0,0)=λm+1∑k=0⌊m/2⌋(-1)m+1-ka2m+1,m,k(cos⁡y0)i-2k[-(sin⁡y0)2k+1]=-f2(2m)(y0,0),
which gives the following relations:
(2.47)f1(2m+1)(-y0,0)=-f1(2m+1)(y0,0),f2(2m)(-y0,0)=-f2(2m)(y0,0).
Using (2.47) for y1=-y0, we obtain the following relations:
(2.48)y1-y0+∑j=1pβj[(-1)jf1(j-1)(y1,0)-f1(j-1)(y0,0)]=-y0-y0+∑j=1pβj[f1(j-1)(-y0,0)-f1(j-1)(y0,0)]=-2(y0+∑j=1pβjf1(j-1)(y0,0))∑j=1pβj[-f2(j-1)(y1,0)-f2(j-1)(y0,0)]=0.
Similarly for y1=y0 using (2.47), we observe that
(2.49)y1-y0+∑j=1pβj[(-1)jf1(j-1)(y1,0)-f1(j-1)(y0,0)]=0,∑j=1pβj[(-1)jf2(j-1)(y1,0)f2(j-1)(y0,0)]=-2∑j=1pβjf2(j-1)(y0,0).
So, our assertions (a) and (b) are proved.
Again, we consider the application of Taylor’s decomposition method to (2.34) on two points xk andxk-1:
(2.50)Yk-Yk-1+∑j=1pαjYk(j)hj-∑j=1qβjYk-1(j)hj=0,
where Yk(j) is the approximate value of Yk(j)(xk). For the computation of the eigenvalues of (1.2), putting h=1 and p=q, the approximation (2.50) gives
(2.51)Y1-Y0+∑j=1p(-1)jβjY1(j)-∑j=1pβjY0(j)=0,
where αj=(-1)jβj. Writing (2.51) with respect to the components and imposing the boundary conditions z0=z(0)=y′(0)=0 and z1=z(1)=y′(1)=0, we have the following equations
(2.52)y1-y0+∑j=1pβj[(-1)jf1(j-1)(y1,0)-f1(j-1)(y0,0)]=0,(2.53)∑j=1pβj[(-1)jf2(j-1)(y1,0)-f2(j-1)(y0,0)]=0.
Using Theorem 2.2(a) for y1=-y0, (2.52) becomes
(2.54)G1(y0,λ)=-2(y0+∑j=1pβjf1(j-1)(y0,0))=0,
and (2.53) is satisfied. For y1=y0, (2.52) is satisfied by Theorem 2.2(b) and (2.53) becomes
(2.55)G2(y0,λ)=-2∑j=1pβjf2(j-1)(y0,0)=0.
From Table1, we observe that there is only trivial initial condition for 0≤λ≤π2, there is one nontrivial initial condition from (2.54) for π2<λ≤4π2, there are n nontrivial initial conditions for n2π2<λ≤(n+1)2π2. These results show that the numerical results obtained using Taylor’s decomposition method agree with the theoretical results of Euler buckling problem given in [2].
Now, we find an approximate solution to the problem(2.56)Y′(x)=F(Y(x)),Y(0)=Y0,
Which corresponds to Euler buckling problem (1.2) for an eigenvalue λ and the initial value y0. Using Taylor’s Decomposition on two points xk-1, xk for p=q then y0≃y(0), z0=z(0). Solving the obtained nonlinear system by Newton’s method, we obtain the approximate value yk of the eigenfunction y(x) at x=xk with O(h2p).
It is clear thatf2(0)(y,z)=siny is Lipschitz in y in 2-dimensional box D. Using the results (Adiyaman and Somali [7, Lemma  2 and Theorem  3]), the global error for (2.50) is bounded by
(2.57)‖Y(xk)-Yk‖≤C0‖Y(0)-Y0‖+C1ξh2pMp+1(2p)!,
where C0=ex¯((2LB(h))/(1-LB(h))), C1=const(C0/L)(1/(1+(β2/β1)h+⋯+βp/β1hp-1)), M=max(y,z)∈D{|f1(0)(y,z)|,|f2(0)(y,z)|}, D is 2-dimensional box in R2, ξ=max{∑k=0⌊(i+1)/2⌋aj,i,k}, j=1,…,2p, i=0,…,p, const is a constant independent of h,p,∥·∥ denotes ∥·∥∞, L=max1≤j≤p{l1,j,l2,j} with l1,j=max1≤j≤p{d1,j,s1,j}, l2,j=max1≤j≤p{d2,j,s2,j}, dk,j=max(y,z)∈D|(∂fk(j)(y,z)/∂y)|, sk,j=max(y,z)∈D|(∂fk(j)(y,z)/∂z)|, k=1,2, and B(h)=L∑j=1pβjhj-1 for some x¯>0.Table 1
Corresponding to the initial valuesy1,0, y2,0, y3,0, and y4,0 for various λ obtained from (2.54) and (2.55).
λy1,0y2,0y3,0y4,015>π21.7471———45>4π22.85781.0092——90>9π23.07182.34130.3236—160>16π23.12722.79992.02390.3771


## 2.1. Application of Taylor’s Decomposition on Two Points for Regular Sturm-Liouville Eigenvalue Problems

We consider the regular Sturm-Liouville eigenvalue problem (1.1) by introducing the new depending variable y′(x)=z(x), (1.1) can be written as
(2.1)Y′(x)=A(x)Y(x),C0Y(x0)+C1Y(xn)=0,
where
(2.2)Y(x)=[y(x)z(x)],A(x)=[01r(x)-λ0],C0=[1000],C1=[0010].
From Ashyralyev and Sobolevskii [8], we will consider the application of Taylor’s decomposition of function Y(x) on two points. We need to find Y(j)(x) for any 1≤j≤pandq. Using the equation Y′(x)=A(x)Y(x), we get
(2.3)Y(j)(x)=Aj(x)Y(x),
with
(2.4)A0(x)=I,A1(x)=A(x),Aj(x)=Aj-1′(x)+Aj-1(x)A(x),2≤j≤p,
where I is the 2×2 identity matrix. By using the structure of the matrix A(x), we obtain the entries of the matrix of
(2.5)Aj(x)=[aj(1,1)(λ;x)aj(1,2)(λ;x)aj(2,1)(λ;x)aj(2,2)(λ;x)]
as in the following formulas:
(2.6)aj(1,1)(λ;x)=∂aj-1(1,1)(λ;x)∂x+(r(x)-λ)aj-2(2,2)(λ;x)=aj-1(2,1)(λ;x),aj(2,2)(λ;x)=∂aj-1(2,2)(λ;x)∂x+aj(1,1)(λ;x),aj(1,2)(λ;x)=aj-1(2,2)(λ;x),aj(2,1)(λ;x)=-∂aj(2,2)(λ;x)∂x+aj+1(2,2)(λ;x),
for 2≤j≤p, where
(2.7)a0(1,1)(λ;x)=1,a1(1,1)(λ;x)=0,a0(1,2)(λ;x)=0,a1(1,2)(λ;x)=1,a0(2,1)(λ;x)=0,a1(2,1)(λ;x)=r(x)-λ,a0(2,2)(λ;x)=1,a1(2,2)(λ;x)=0.From the theorem given in Ashyralyev and Sobolevskii [8], we have the following relation:
(2.8)Y(xk)-Y(xk-1)+∑j=1pαjY(j)(xk)hj-∑j=1qβjY(j)(xk-1)hj=(-1)p(p+q)!∫xk-1xk(xk-s)q(s-xk-1)pY(p+q+1)(s)ds
on the uniform grid
(2.9)[x0,xn]h={xk=x0+kh,k=0,1,…,n,nh=xn-x0,n∈N},
where
(2.10)αj=(p+q-j)!p!(-1)j(p+q)!j!(p-j)!,1≤j≤p,βj=(p+q-j)!q!(p+q)!j!(q-j)!,1≤j≤q.
Rewriting (2.8) by neglecting the last term, we obtain the single-step difference scheme of (p+q)-order of accuracy for the approximate solution of problem (2.1):
(2.11)Yk-Yk-1+∑j=1pαjAj(xk)Ykhj-∑j=1qβjAj(xk-1)Yk-1hj=0,
where
(2.12)Yk=[ykzk]
is the approximate value of Y(xk). For the simple computation, let p=q, then we have
(2.13)(I+∑j=1pαjAj(xk)hj)Yk=(I+∑j=1p(-1)jαjAj(xk-1)hj)Yk-1,
where βj=(-1)jαj. Letting M(xk)=(I+∑j=1pαjAj(xk)hj) and N(xk-1)=(I+∑j=1p(-1)jαjAj(xk-1)hj), we write
(2.14)Yk=M-1(xk)N(xk-1)Yk-1.
Since the accuracy and convergence of the method not only depend on h, they also depend on p, we can increase the order of accuracy by increasing p for fixed h. So h is chosen as length of the whole interval as follows. Now, taking h=xn-x0 gives
(2.15)Y1=M-1(xn)N(x0)Y0,
and substituting into the boundary condition of (2.1), we get
(2.16)(C1M-1(xn)N(x0)+C0)Y0=0.
To obtain a nontrivial solution Y0, we must have the following equation:
(2.17)det⁡(C1M-1(xn)N(x0)+C0)=0.
Defining
(2.18)M(xn)=[m11m12m21m22],N(x0)=[n11n12n21n22],
we have the following statement
(2.19)m22n12-m12n22=0.
Since
(2.20)M(xn)=[1+∑j=1pαjaj(1,1)(λ;xn)hj∑j=1pαjaj(1,2)(λ;xn)hj∑j=1pαjaj(2,1)(λ;xn)hj1+∑j=1pαjaj(2,2)(λ;xn)hj],N(x0)=[1+∑j=1p(-1)jαjaj(1,1)(λ;x0)hj∑j=1p(-1)jαjaj(1,2)(λ;x0)hj∑j=1p(-1)jαjaj(2,1)(λ;x0)hj1+∑j=1p(-1)jαjaj(2,2)(λ;x0)hj],
using the entries m12,m22,n12, and n22 of the above matrices and the properties of the entries of Aj(x), we obtain (2.19) in terms of λ:
(2.21)F(λ)=m22n12-m12n22=(1+∑j=1pαjaj(2,2)(λ;xn)hj)(∑j=1p(-1)jαjaj-1(2,2)(λ;x0)hj)-(∑j=1pαjaj-1(2,2)(λ;xn)hj)(1+∑j=1p(-1)jαjaj(2,2)(λ;x0)hj).
Solving nonlinear equation F(λ)=0 by Newton’s method, we find the approximate eigenvalues. This method appears to require a separate calculation for the eigenfunctions.To find the corresponding eigenfunctions of the regular Sturm-Liouville eigenvalue problem (2.1), we substitute the eigenvalue to (2.1) and we solve the obtained boundary value problem by Taylor’s decomposition method on two points xk-1 and xk with the uniform grid [0,1]h for p=q. Then, we get a homogeneous linear equation system of 2n equations with 2n unknown z0,y1,z1,y2,z2…,yn-1,zn-1,zn which are the approximated values of y′(x0),y(x1),y′(x1),y(x2),y′(x2),…,y(xn-1),y′(xn-1),y′(xn), respectively. Solving the 2n×2n homogeneous system, we obtain approximate values of the eigenfunction and its derivative of (1.1) at the point x=xk.


## 2.2. Error Analysis for Regular Sturm-Liouville Problem Whenr(x)=c

In this section, we will show the convergence of the method for eigenfunctions with the constant functionr(x)=c by obtaining approximate value of eigenfunction at the point x¯∈[x0,xn] of the problem (1.1). Without loss of generality, we may choose r(x)=0, then Aj(x)=Aj, that is, aj(2,2)(λ;xn)=aj(2,2)(λ;0)=aj(2,2)(λ). Using (2.6), we can find explicit values of aj(1,1), aj(2,2) as follows:
(2.22)a2j(1,1)=(-1)jλj,a2j(2,2)=(-1)jλj,a2j+1(1,1)=0,a2j+1(2,2)=0,j≥0.
This yields
(2.23)m22=1+∑j=1⌊p/2⌋α2j(-1)jλjh2j,n22=m22,m12=∑j=0⌊(p-1)/2⌋α2j+1(-1)jλjh2j+1,n12=-m12,m11=m22,m21=-λm12.
Using (2.14) for k=1, we have
(2.24)Y1=M-1(x¯)N(x0)Y0,
where Y0 and Y1 are the approximated values of Y(x0) and Y(x¯), respectively, with the stepsize h=x¯-x0:(2.25)Y1=z0det⁡(M)[-2m22m12-λ(m12)2+(m22)2].
The first component of the above vector (2.25) gives the approximate eigenfunction y1, and the second component of the above vector (2.25) gives the derivative of the approximate eigenfunction z1 of the regular Sturm-Liouville problem (1.1) at x¯. Now, we will show that y1 and z1 converge to exact functions y(x¯) and y′(x¯), respectively, as p→∞.Using the Stirling’s Formulan!≈2πn(n+1)/2e-n for αj in (2.10), we obtain
(2.26)αj=(2p-j)!p!(-1)j(2p)!j!(p-j)!≈(-1)j1j!12j(p-j/2p-j)(p-j+1)/2(p-j/2p)p/2.
This gives
(2.27)lim⁡p→∞αj=(-1)j1j!12j.
Thus,
(2.28)lim⁡p→∞m22=lim⁡p→∞(1+∑j=1⌊p/2⌋α2j(-1)jλjh2j)=1+∑j=1∞1(2j)!122j(-1)jλjh2j=∑j=0∞(-1)j(λh2)2j1(2j)!=cos⁡(λ)h2.
By using the same idea, we obtain
(2.29)lim⁡p→∞m12=lim⁡p→∞(∑j=0⌊(p-1)/2⌋(-1)jα2j+1λjh2j+1)=1λsin⁡(λ)h2.
It follows from (2.28) and (2.29) that
(2.30)lim⁡p→∞det⁡(M)=m222+λm122=cos⁡2(λ)h2+λ(1λsin⁡(λ)h2)2=1.
Hence, for r(x)=0, the approximate eigenfunction of (1.1) to the corresponding eigenvalue λ converges to exact eigenfunction:
(2.31)lim⁡p→∞y1=2z0det⁡(M)1λ(cos⁡(λh2))(sin⁡(λh2))=z0λsin⁡(λ(x¯-x0)).
Since we have z(x)=y′(x), the derivative of approximate eigenfunction of (1.1) to the corresponding eigenvalue λ converges to derivative of the exact solution:
(2.32)lim⁡p→∞z1=z0det⁡(M)((-λ)1λsin⁡2(λh2)+cos⁡2(λh2))=z0cos⁡(λ(x¯-x0)),
where λ=k2π2,k=1,2,….This demonstration shows that approximate eigenfunction and eigenvalue converges to exact one asp→∞ for fixed step-size “h.”


## 2.3. Taylor’s Decomposition Method to the Euler Buckling Problem

For convenience, we introduce the following notations as in (2.1) and Adiyaman and Somali [7]:
(2.33)Y(x)=[y(x)z(x)],F(Y(x))=[f1(0)(y,z)f2(0)(y,z)],C0=[0100],C1=[0001],f1(0)(y,z)=z,f2(0)(y,z)=-λsin⁡y.
Thus, the Euler Buckling Problem (1.2) can be written in the form:
(2.34)Y′(x)=F(Y(x)),0<x<1,C0Y(0)+C1Y(1)=0,
Defining the following recurrence relations for j=1,…,2p:
(2.35)fi(j)(y,z)=z∂fi(j-1)(y,z)∂y-λsin⁡y∂fi(j-1)(y,z)∂z,i=1,2,
we obtain
(2.36)Y(j)(x)=[f1(j-1)(y,z)f2(j-1)(y,z)]=[f2(j-2)(y,z)f2(j-1)(y,z)]forj=2,…,2p+1.
We first give the following lemma which defines f2(j-1)(y,z) explicitly.Lemma 2.1.
Forj=0,…,2p, let f2(j)(y,z) satisfy the recurrence relation (2.35) with f2(0)(y,z)=-λsiny. Then
(2.37)f2(2m)(y,z)=∑i=0mλi+1z2m-2i∑k=0⌊i/2⌋(-1)m+1-ka2m+1,i,k(cos⁡y)i-2k(sin⁡y)2k+1,
where
(2.38)a2m+1,i,k={1,i=0,k=0,(2k+2)a2m,i,k+1+(i+1-2k)a2m,i,k+(2m+1-2i)a2m,i-1,k,1≤i≤m,0≤k≤i2,0,else,
for m=0,…,p, and
(2.39)f2(2m+1)(y,z)=∑i=0mλi+1z2m+1-2i∑k=0⌊(i+1)/2⌋(-1)m+1-ka2m+2,i,k(cos⁡y)i+1-2k(sin⁡y)2k,
where
(2.40)a2m+2,i,k={1,i=0,k=0,(2k+1)a2m+1,i,k+(i-2k)a2m+1,i,k-1+(2m+2-2i)a2m+1,i-1,k-1,1≤i≤m,0≤k≤i20,else,
for m=0,…,p-1.Proof.
The proof follows induction argument based on (2.35).Theorem 2.2.
Iff1(j)(y,z) and  f2(j)(y,z), are sufficiently smooth and satisfy (2.36), (2.37), and (2.39) then the following relations hold: (a) it holds that
(2.41)y1-y0+∑j=1pβj[(-1)jf1(j-1)(y1,0)-f1(j-1)(y0,0)]=-2(y0+∑j=1pβjf1(j-1)(y0,0)),∑j=1pβj[(-1)jf2(j-1)(y1,0)-f2(j-1)(y0,0)]=0,
for y1=-y0, (b) it holds that
(2.42)∑j=1pβj[(-1)jf2(j-1)(y1,0)-f2(j-1)(y0,0)]=-2∑j=1pβjf2(j-1)(y0,0),y1-y0+∑j=1pβj[(-1)jf1(j-1)(y1,0)-f1(j-1)(y0,0)]=0,
for y1=y0.Proof.
Letj=2m+1 for m=0,…,p, then f1(2m+1)(y,z) becomes
(2.43)f1(2m)(y,z)=f2(2m-1)(y,z)=∑i=0(m-1)λi+1z2m-2i-1∑k=0⌊(i+1)/2⌋(-1)m-ka2m,i,k(cos⁡y)i+1-2k(sin⁡y)2k,
by Lemma 2.1. Since all terms of previous sum contain z, f12m(y,0)=f22m-1(y,0)=0 for m=1,…,p, hence, we get the following equations:
(2.44)y1-y0+∑j=1pβj[(-1)jf1(j-1)(y1,0)-f1(j-1)(y0,0)]=y1-y0+∑j=1pβj[f1(j-1)(y1,0)-f1(j-1)(y0,0)],∑j=1pβj[(-1)jf2(j-1)(y1,0)-f2(j-1)(y0,0)]=∑j=1pβj[-f2(j-1)(y1,0)-f2(j-1)(y0,0)].
Letting j=2m+2 for m=0,…,p-1 and using (2.37), we get
(2.45)f1(2m+2-1)(y,z)=f2(2m)(y,z)=∑i=0mλi+1z2m-2i∑k=0⌊i/2⌋(-1)m+1-ka2m+1,i,k(cos⁡y)i-2k(sin⁡y)2k+1.
Substituting the value z=0 into (2.45), we obtain
(2.46)f1(2m+1)(-y0,0)=f22m(-y0,0)=λm+1∑k=0⌊m/2⌋(-1)m+1-ka2m+1,m,k(cos⁡y0)i-2k[-(sin⁡y0)2k+1]=-f2(2m)(y0,0),
which gives the following relations:
(2.47)f1(2m+1)(-y0,0)=-f1(2m+1)(y0,0),f2(2m)(-y0,0)=-f2(2m)(y0,0).
Using (2.47) for y1=-y0, we obtain the following relations:
(2.48)y1-y0+∑j=1pβj[(-1)jf1(j-1)(y1,0)-f1(j-1)(y0,0)]=-y0-y0+∑j=1pβj[f1(j-1)(-y0,0)-f1(j-1)(y0,0)]=-2(y0+∑j=1pβjf1(j-1)(y0,0))∑j=1pβj[-f2(j-1)(y1,0)-f2(j-1)(y0,0)]=0.
Similarly for y1=y0 using (2.47), we observe that
(2.49)y1-y0+∑j=1pβj[(-1)jf1(j-1)(y1,0)-f1(j-1)(y0,0)]=0,∑j=1pβj[(-1)jf2(j-1)(y1,0)f2(j-1)(y0,0)]=-2∑j=1pβjf2(j-1)(y0,0).
So, our assertions (a) and (b) are proved.
Again, we consider the application of Taylor’s decomposition method to (2.34) on two points xk andxk-1:
(2.50)Yk-Yk-1+∑j=1pαjYk(j)hj-∑j=1qβjYk-1(j)hj=0,
where Yk(j) is the approximate value of Yk(j)(xk). For the computation of the eigenvalues of (1.2), putting h=1 and p=q, the approximation (2.50) gives
(2.51)Y1-Y0+∑j=1p(-1)jβjY1(j)-∑j=1pβjY0(j)=0,
where αj=(-1)jβj. Writing (2.51) with respect to the components and imposing the boundary conditions z0=z(0)=y′(0)=0 and z1=z(1)=y′(1)=0, we have the following equations
(2.52)y1-y0+∑j=1pβj[(-1)jf1(j-1)(y1,0)-f1(j-1)(y0,0)]=0,(2.53)∑j=1pβj[(-1)jf2(j-1)(y1,0)-f2(j-1)(y0,0)]=0.
Using Theorem 2.2(a) for y1=-y0, (2.52) becomes
(2.54)G1(y0,λ)=-2(y0+∑j=1pβjf1(j-1)(y0,0))=0,
and (2.53) is satisfied. For y1=y0, (2.52) is satisfied by Theorem 2.2(b) and (2.53) becomes
(2.55)G2(y0,λ)=-2∑j=1pβjf2(j-1)(y0,0)=0.
From Table1, we observe that there is only trivial initial condition for 0≤λ≤π2, there is one nontrivial initial condition from (2.54) for π2<λ≤4π2, there are n nontrivial initial conditions for n2π2<λ≤(n+1)2π2. These results show that the numerical results obtained using Taylor’s decomposition method agree with the theoretical results of Euler buckling problem given in [2].
Now, we find an approximate solution to the problem(2.56)Y′(x)=F(Y(x)),Y(0)=Y0,
Which corresponds to Euler buckling problem (1.2) for an eigenvalue λ and the initial value y0. Using Taylor’s Decomposition on two points xk-1, xk for p=q then y0≃y(0), z0=z(0). Solving the obtained nonlinear system by Newton’s method, we obtain the approximate value yk of the eigenfunction y(x) at x=xk with O(h2p).
It is clear thatf2(0)(y,z)=siny is Lipschitz in y in 2-dimensional box D. Using the results (Adiyaman and Somali [7, Lemma  2 and Theorem  3]), the global error for (2.50) is bounded by
(2.57)‖Y(xk)-Yk‖≤C0‖Y(0)-Y0‖+C1ξh2pMp+1(2p)!,
where C0=ex¯((2LB(h))/(1-LB(h))), C1=const(C0/L)(1/(1+(β2/β1)h+⋯+βp/β1hp-1)), M=max(y,z)∈D{|f1(0)(y,z)|,|f2(0)(y,z)|}, D is 2-dimensional box in R2, ξ=max{∑k=0⌊(i+1)/2⌋aj,i,k}, j=1,…,2p, i=0,…,p, const is a constant independent of h,p,∥·∥ denotes ∥·∥∞, L=max1≤j≤p{l1,j,l2,j} with l1,j=max1≤j≤p{d1,j,s1,j}, l2,j=max1≤j≤p{d2,j,s2,j}, dk,j=max(y,z)∈D|(∂fk(j)(y,z)/∂y)|, sk,j=max(y,z)∈D|(∂fk(j)(y,z)/∂z)|, k=1,2, and B(h)=L∑j=1pβjhj-1 for some x¯>0.Table 1
Corresponding to the initial valuesy1,0, y2,0, y3,0, and y4,0 for various λ obtained from (2.54) and (2.55).
λy1,0y2,0y3,0y4,015>π21.7471———45>4π22.85781.0092——90>9π23.07182.34130.3236—160>16π23.12722.79992.02390.3771


## 3. Numerical Results

### 3.1. Numerical Results for Regular Sturm-Liouville Eigenvalue Problems

We consider three regular Sturm-Liouville eigenvalue problems, one of them has polynomial coefficients and the others have periodic coefficients taken from Bujurke et al. [9] and Andrew [10].Example 3.1.
Consider the Titchmarch equation:(3.1)y′′+(λ-x2n)y(x)=0,y(0)=y(1)=0,
where n is a nonnegative integer. We obtain the numerical solutions taking n=0,2. The accuracy of the method is tested by comparing with the exact solution which exists when n=0 and finite-difference method (FDM) solution and Haar wavelet series method (HWSM) solution when n=2.
Tables2 and 3 give computed eigenvalues and solution y(x) of Titchmarch problem using Taylor’s decomposition method (TDM) with different values of p, HWSM and FDM for n=0,2, the integer parameter in Titchmarch problem. In Table 2, it is easily seen that the error between approximate eigenfunction and exact eigenfunction decreases as p increases or the step-size decreases or both happen. So, we can find good approximation to eigenfunctions for relatively large step-sizes by increasing p. In Table 3m is the number of intervals. Table 4 gives the errors between exact and approximate eigenvalues for fixed step-size h=1 for n=0. Notice that, as p increases, the accuracy of approximation almost doubles in digits which demonstrates a fast convergence.Table 2
Comparison of the first eigenvalue and solutions of Example3.1 using Taylor’s decomposition method, exact values, and Table 4 from [9], when n=0.
xHWSMFDMTDMExactErrors of TDMErrors of TDMErrors of  TDM(h=1/16)(h=1/16)(p=3) (h=1/16)(p=3)   (h=1/16)(p=4) (h=1/8)(p=5) (h=1/4)000000000.06250.275210.2785990.2759990.2758991.54E-10——0.1250.541810.5411960.5411960.5411962.91E-101.13E-11—0.18750.785490.7856950.7856950.7856953.93E-10——0.251.004821114.45E-101.74E-116.87E-120.31251.178511.175881.175881.175884.37E-10——0.3751.312851.306561.306561.306563.61E-101.41E-11—0.43751.383761.387041.387041.387042.15E-10——0.51.411031.414211.414211.414212.22E-162.22E-1600.56251.383761.387041.387041.387042.76E-10——0.6251.312851.306561.306561.306566.03E-102.35E-11—0.68751.178511.175881.175881.175889.63E-10——0.751.004821111.33E-95.22E-112.06E-110.81250.785490.7856950.7856950.7856951.70E-9——0.8750.541810.5411960.5411960.5411962.03E-97.96E-11—0.93750.275210.2758990.2759990.2758992.31E-9——10000000λ1=10.9334 (HWSM), 10.8379 (FDM), 10.8696 (TDM), 10.8696 (Exact).Table 3
Comparison of the first eigenvalue and solutions of Example3.1 using TDM and Table 4 from [9], when p=16, n=2, and h=0.0625.
xHWSMFDMTDM00000.06250.275210.277560.2775630.1250.541810.544340.5443370.18750.789490.789960.7899530.251.004851.004881.004870.31251.181531.180751.180740.3751.312861.310821.310760.43751.393721.389961.389940.51.421021.415271.415290.56251.393711.385911.385980.6251.312851.303231.303340.68751.181541.180661.170810.751.00480.993610.9937920.81250.779490.779170.7793570.8750.534810.535770.5359340.93750.277260.272770.2728781000λ1=10.3452 (HWSM), 9.95067 (FDM), 9.98317 (TDM).Table 4
The errors between exact and approximate fundamental eigenvalue for variousp and for h=1 and n=0 in Example 3.1.
p2346791011Errors2.13E01.30E-15.49E-32.60E-63.38E-82.60E-121.59E-141.77E-15Example 3.2.
Consider the Mathieu’s equation:(3.2)y′′+(λ-2θcos⁡(2x))y=0,y(0)=y(π)=0.
We will solve these two problems approximately using Taylor’s decomposition method (TDM), and we will compare our results with the results in Bujurke et al. [9]. Bujurke et al. [9] solved Examples 3.1 and 3.2 approximately using Haar wavelets. They transform the interval [0,π] to [0,1] because of the properties of Haar wavelets. So to compare the results we normalize the interval [0,π] by using x=πt, Mathieu’s equation in Example 3.2 transformed into
(3.3)y′′+(π2λ-2π2θcos⁡(2πt))y=0,y(0)=y(1)=0.
The eigenvalues for a fixed value forθ=5 are obtained in Table 6 which gives the asymptotic behavior of higher eigenvalues of Mathieu’s equation, and these eigenvalues are λn=n2+O(1). This result agrees with the classical theorem on asymptoticity of the eigenvalues limn→∞λn1/2/n=1 from van Brunt [11]. Figure 1 demonstrates that the nth eigenfunction has n-1 zeros in (0,1) which is consistent with the relevant graph in Bujurke et al. [9]. The selected values of parameter θ shifts the symmetry of the solutions and this property is given in Figure 2.Figure 1
Higher eigenfunctions of Mathieu’s equation for a fixed parameterθ=5.Figure 2
Solutions of Mathieu’s equation for different parameters ofθ.Example 3.3.
Consider the equation(3.4)-y′′+(2βcos⁡2x+β2sin⁡22x)y=λy,y′(0)=y′(π)=0.
We give the comparison of approximate eigenvalues obtained using Taylor’s Decomposition method with the approximate eigenvalues obtained using Numerov’s method Andrew [10] for β=10 and β=20 in Tables 7 and 8. The values shown as the “exact” λk and the corrected approximate eigenvalues Λk obtained using Numerov’s method for step-sizes h=40 and h=80 in Tables 7 and 8 are taken from Andrew [10]. The values shown as μk are evaluated using Taylor’s Decomposition method for p=70 and 80 in Table 7 and for p=100 and 110 in Table 8. From the tables, it can be seen that Taylor’s Decomposition method approximates small eigenvalues with high-order accuracy without using any correction.
In comparison to Example3.1, the estimation of eigenvalues for Examples 3.2 and 3.3 is more complicated. But Example 3.1 is important to show the high accuracy of the method while calculating the eigenfunctions for relatively large step-sizes. Other two examples show the accuracy of the method while calculating the eigenvalues for large step-sizes which equal to whole interval.
In Table5, the observed orders ord(h) are computed using the following formula
(3.5)ord(h)=log⁡⁡((y4h-y2h)/(y2h-yh))log⁡2,
where y4h, y2h, and yh are the approximated value of eigenfunctions at xk to the corresponding eigenvalue λ when the problems are solved with step sizes 4h, 2h, and h respectively. The observed orders given in Table 5 well confirm the theoretical results. That is, the order of TDM is order of 2p.
The numerical calculations and all figures in this work are performed using Mathematica.Table 5
Observed orders of Example3.1 for n=2 at x=1/2 using Taylor’s decomposition method.
p=2p=3p=4ord(1/8)4.232477.9060710.6574ord(1/16)3.861535.221418.13719ord(1/32)3.986745.930747.63932Table 6
Comparison of higher eigenvalues for Mathieu’s equation obtained from FDM, HWSM, and TDM corresponding toθ=5.
nn2λn(FDM)λn (HWSM)λn (TDM)11−5.7311−5.4665−5.79008242.09922.61612.09946399.23659.42279.2363341616.64816.370716.648252525.51124.147125.510863636.35936.657736.3589λ1=-5.46653 (HWSM), −5.73115 (FDM), −5.79008 (TDM).Table 7
Comparison of the exact eigenvalues with approximate eigenvalues obtained from Numerov’s method (Λk) with correction and TDM (μk) for Example 3.3 corresponding to β=10.
kλkλk-Λk(40)λk-Λk(80)λk-μk70λk-μk8010.00000002.32E-31.43E-43.63E-37.60E-5237.75962857.52E-34.42E-49.15E-43.81E-6337.80590021.61E-39.93E-55.66E-71.12E-7437.85259956.88E-34.43E-48.97E-43.85E-5570.54750973.22E-21.97E-31.57E-23.74E-4692.65381772.25E-21.37E-32.25E-32.20E-5796.20581591.20E-27.28E-42.18E-31.78E-48102.2543473.45E-22.09E-34.38E-32.44E-4h=π for Taylor’s decomposition method.Table 8
Comparison of the exact eigenvalues with approximate eigenvalues obtained from Numerov’s method (Λk) with correction and TDM (μk) for Example 3.3 corresponding to β=20.
kλkλk-Λk(40)λk-Λk(80)λk-μk100λk-μk11010.00000001.93E-21.17E-31.04E-17.64E-3277.91619431.22E-17.38E-34.37E-26.98E-2377.91619571.63E-29.95E-49.89E-41.67E-4477.91619721.62E-29.89E-41.13E-17.48E-25151.4632243.90E-12.33E-27.04E-11.56E-1h=π for Taylor’s decomposition method.


### 3.2. Numerical Results for Euler Buckling Problem

The approximate solutions of Euler Buckling problem forλ=15, λ=45, λ=90, and λ=160 generated using Taylor’s Decomposition method for step size h=1/20 are illustrated in Figures 3, 4, 5, and 6, respectively.Figure 3
Solution of (1.2) corresponding to the initial values y0 for π2≤λ=15<4π2.Figure 4
Solution of (1.2) corresponding to the initial values y0 for 4π2≤λ=45<9π2.Figure 5
Solution of (1.2) corresponding to the initial values y0 for 9π2≤λ=90<16π2.Figure 6
Solution of (1.2) corresponding to the initial values y0 for 16π2≤λ=160<25π2.


## 3.1. Numerical Results for Regular Sturm-Liouville Eigenvalue Problems

We consider three regular Sturm-Liouville eigenvalue problems, one of them has polynomial coefficients and the others have periodic coefficients taken from Bujurke et al. [9] and Andrew [10].Example 3.1.
Consider the Titchmarch equation:(3.1)y′′+(λ-x2n)y(x)=0,y(0)=y(1)=0,
where n is a nonnegative integer. We obtain the numerical solutions taking n=0,2. The accuracy of the method is tested by comparing with the exact solution which exists when n=0 and finite-difference method (FDM) solution and Haar wavelet series method (HWSM) solution when n=2.
Tables2 and 3 give computed eigenvalues and solution y(x) of Titchmarch problem using Taylor’s decomposition method (TDM) with different values of p, HWSM and FDM for n=0,2, the integer parameter in Titchmarch problem. In Table 2, it is easily seen that the error between approximate eigenfunction and exact eigenfunction decreases as p increases or the step-size decreases or both happen. So, we can find good approximation to eigenfunctions for relatively large step-sizes by increasing p. In Table 3m is the number of intervals. Table 4 gives the errors between exact and approximate eigenvalues for fixed step-size h=1 for n=0. Notice that, as p increases, the accuracy of approximation almost doubles in digits which demonstrates a fast convergence.Table 2
Comparison of the first eigenvalue and solutions of Example3.1 using Taylor’s decomposition method, exact values, and Table 4 from [9], when n=0.
xHWSMFDMTDMExactErrors of TDMErrors of TDMErrors of  TDM(h=1/16)(h=1/16)(p=3) (h=1/16)(p=3)   (h=1/16)(p=4) (h=1/8)(p=5) (h=1/4)000000000.06250.275210.2785990.2759990.2758991.54E-10——0.1250.541810.5411960.5411960.5411962.91E-101.13E-11—0.18750.785490.7856950.7856950.7856953.93E-10——0.251.004821114.45E-101.74E-116.87E-120.31251.178511.175881.175881.175884.37E-10——0.3751.312851.306561.306561.306563.61E-101.41E-11—0.43751.383761.387041.387041.387042.15E-10——0.51.411031.414211.414211.414212.22E-162.22E-1600.56251.383761.387041.387041.387042.76E-10——0.6251.312851.306561.306561.306566.03E-102.35E-11—0.68751.178511.175881.175881.175889.63E-10——0.751.004821111.33E-95.22E-112.06E-110.81250.785490.7856950.7856950.7856951.70E-9——0.8750.541810.5411960.5411960.5411962.03E-97.96E-11—0.93750.275210.2758990.2759990.2758992.31E-9——10000000λ1=10.9334 (HWSM), 10.8379 (FDM), 10.8696 (TDM), 10.8696 (Exact).Table 3
Comparison of the first eigenvalue and solutions of Example3.1 using TDM and Table 4 from [9], when p=16, n=2, and h=0.0625.
xHWSMFDMTDM00000.06250.275210.277560.2775630.1250.541810.544340.5443370.18750.789490.789960.7899530.251.004851.004881.004870.31251.181531.180751.180740.3751.312861.310821.310760.43751.393721.389961.389940.51.421021.415271.415290.56251.393711.385911.385980.6251.312851.303231.303340.68751.181541.180661.170810.751.00480.993610.9937920.81250.779490.779170.7793570.8750.534810.535770.5359340.93750.277260.272770.2728781000λ1=10.3452 (HWSM), 9.95067 (FDM), 9.98317 (TDM).Table 4
The errors between exact and approximate fundamental eigenvalue for variousp and for h=1 and n=0 in Example 3.1.
p2346791011Errors2.13E01.30E-15.49E-32.60E-63.38E-82.60E-121.59E-141.77E-15Example 3.2.
Consider the Mathieu’s equation:(3.2)y′′+(λ-2θcos⁡(2x))y=0,y(0)=y(π)=0.
We will solve these two problems approximately using Taylor’s decomposition method (TDM), and we will compare our results with the results in Bujurke et al. [9]. Bujurke et al. [9] solved Examples 3.1 and 3.2 approximately using Haar wavelets. They transform the interval [0,π] to [0,1] because of the properties of Haar wavelets. So to compare the results we normalize the interval [0,π] by using x=πt, Mathieu’s equation in Example 3.2 transformed into
(3.3)y′′+(π2λ-2π2θcos⁡(2πt))y=0,y(0)=y(1)=0.
The eigenvalues for a fixed value forθ=5 are obtained in Table 6 which gives the asymptotic behavior of higher eigenvalues of Mathieu’s equation, and these eigenvalues are λn=n2+O(1). This result agrees with the classical theorem on asymptoticity of the eigenvalues limn→∞λn1/2/n=1 from van Brunt [11]. Figure 1 demonstrates that the nth eigenfunction has n-1 zeros in (0,1) which is consistent with the relevant graph in Bujurke et al. [9]. The selected values of parameter θ shifts the symmetry of the solutions and this property is given in Figure 2.Figure 1
Higher eigenfunctions of Mathieu’s equation for a fixed parameterθ=5.Figure 2
Solutions of Mathieu’s equation for different parameters ofθ.Example 3.3.
Consider the equation(3.4)-y′′+(2βcos⁡2x+β2sin⁡22x)y=λy,y′(0)=y′(π)=0.
We give the comparison of approximate eigenvalues obtained using Taylor’s Decomposition method with the approximate eigenvalues obtained using Numerov’s method Andrew [10] for β=10 and β=20 in Tables 7 and 8. The values shown as the “exact” λk and the corrected approximate eigenvalues Λk obtained using Numerov’s method for step-sizes h=40 and h=80 in Tables 7 and 8 are taken from Andrew [10]. The values shown as μk are evaluated using Taylor’s Decomposition method for p=70 and 80 in Table 7 and for p=100 and 110 in Table 8. From the tables, it can be seen that Taylor’s Decomposition method approximates small eigenvalues with high-order accuracy without using any correction.
In comparison to Example3.1, the estimation of eigenvalues for Examples 3.2 and 3.3 is more complicated. But Example 3.1 is important to show the high accuracy of the method while calculating the eigenfunctions for relatively large step-sizes. Other two examples show the accuracy of the method while calculating the eigenvalues for large step-sizes which equal to whole interval.
In Table5, the observed orders ord(h) are computed using the following formula
(3.5)ord(h)=log⁡⁡((y4h-y2h)/(y2h-yh))log⁡2,
where y4h, y2h, and yh are the approximated value of eigenfunctions at xk to the corresponding eigenvalue λ when the problems are solved with step sizes 4h, 2h, and h respectively. The observed orders given in Table 5 well confirm the theoretical results. That is, the order of TDM is order of 2p.
The numerical calculations and all figures in this work are performed using Mathematica.Table 5
Observed orders of Example3.1 for n=2 at x=1/2 using Taylor’s decomposition method.
p=2p=3p=4ord(1/8)4.232477.9060710.6574ord(1/16)3.861535.221418.13719ord(1/32)3.986745.930747.63932Table 6
Comparison of higher eigenvalues for Mathieu’s equation obtained from FDM, HWSM, and TDM corresponding toθ=5.
nn2λn(FDM)λn (HWSM)λn (TDM)11−5.7311−5.4665−5.79008242.09922.61612.09946399.23659.42279.2363341616.64816.370716.648252525.51124.147125.510863636.35936.657736.3589λ1=-5.46653 (HWSM), −5.73115 (FDM), −5.79008 (TDM).Table 7
Comparison of the exact eigenvalues with approximate eigenvalues obtained from Numerov’s method (Λk) with correction and TDM (μk) for Example 3.3 corresponding to β=10.
kλkλk-Λk(40)λk-Λk(80)λk-μk70λk-μk8010.00000002.32E-31.43E-43.63E-37.60E-5237.75962857.52E-34.42E-49.15E-43.81E-6337.80590021.61E-39.93E-55.66E-71.12E-7437.85259956.88E-34.43E-48.97E-43.85E-5570.54750973.22E-21.97E-31.57E-23.74E-4692.65381772.25E-21.37E-32.25E-32.20E-5796.20581591.20E-27.28E-42.18E-31.78E-48102.2543473.45E-22.09E-34.38E-32.44E-4h=π for Taylor’s decomposition method.Table 8
Comparison of the exact eigenvalues with approximate eigenvalues obtained from Numerov’s method (Λk) with correction and TDM (μk) for Example 3.3 corresponding to β=20.
kλkλk-Λk(40)λk-Λk(80)λk-μk100λk-μk11010.00000001.93E-21.17E-31.04E-17.64E-3277.91619431.22E-17.38E-34.37E-26.98E-2377.91619571.63E-29.95E-49.89E-41.67E-4477.91619721.62E-29.89E-41.13E-17.48E-25151.4632243.90E-12.33E-27.04E-11.56E-1h=π for Taylor’s decomposition method.


## 3.2. Numerical Results for Euler Buckling Problem

The approximate solutions of Euler Buckling problem forλ=15, λ=45, λ=90, and λ=160 generated using Taylor’s Decomposition method for step size h=1/20 are illustrated in Figures 3, 4, 5, and 6, respectively.Figure 3
Solution of (1.2) corresponding to the initial values y0 for π2≤λ=15<4π2.Figure 4
Solution of (1.2) corresponding to the initial values y0 for 4π2≤λ=45<9π2.Figure 5
Solution of (1.2) corresponding to the initial values y0 for 9π2≤λ=90<16π2.Figure 6
Solution of (1.2) corresponding to the initial values y0 for 16π2≤λ=160<25π2.


## 4. Conclusion

In this paper, we have described Taylor’s Decomposition method for regular Sturm-Liouville eigenvalue problems with Dirichlet and Neumann boundary conditions to obtain approximate eigenvalues and eigenfunctions and for Euler Buckling Problem to obtain approximate initial values and eigenfunctions. The obtained results for Euler Buckling problem give the behavior of eigenvalues and corresponding eigenfunctions with high-order accuracy without using small stepsize. We have seen that these results agree with the theoretical aspects. This method can be extended to solve regular Sturm-Liouville eigenvalue problems with Robin (mixed) boundary conditions and to some nonlinear eigenvalue problems to investigate the behavior of the eigenvalues and eigenfunction. However, this method is best suited to find small eigenvalues for the other nonlinear problems in literature. One possible method of improving its efficiency for higher eigenvalues may be to follow the ideas of [10, 12, 13] and for eigenvalue problems for partial differential equations given in elsewhere [14–18].



---

*Source: 697013-2012-06-05.xml*