Solving an Inverse Sturm-Liouville Problem by a Lie-Group Method

Volume 2008, Article ID 749865,18pages doi:10.1155/2008/749865

Research Article

Solving an Inverse Sturm-Liouville Problem by a Lie-Group Method

Chein-Shan Liu^{1, 2}

1Department of Mechanical and Mechatronic Engineering, Taiwan Ocean University, Keelung 20224, Taiwan

2Department of Harbor and River Engineering, Taiwan Ocean University, Keelung 20224, Taiwan

Correspondence should be addressed to Chein-Shan Liu,[email protected] Received 8 September 2007; Revised 21 December 2007; Accepted 29 January 2008 Recommended by Colin Rogers

Solving an inverse Sturm-Liouville problem requires a mathematical process to determine unknown function in the Sturm-Liouville operator from given data in addition to the boundary values. In this paper, we identify a Sturm-Liouville potential function by using the data of one eigenfunction and its corresponding eigenvalue, and identify a spatial-dependent unknown function of a Sturm- Liouville differential operator. The method we employ is to transform the inverse Sturm-Liouville problem into a parameter identification problem of a heat conduction equation. Then a Lie-group estimation method is developed to estimate the coefficients in a system of ordinary differential equations discretized from the heat conduction equation. Numerical tests confirm the accuracy and efficiency of present approach. Definite and random disturbances are also considered when com- paring the present method with that by using a technique of numerical differentiation.

Copyrightq2008 Chein-Shan Liu. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

The problem to describe the interaction between colliding particles is a fundamental one in the physics of particle, where the identification of Schr ¨odinger operator is utmost important. It is one sort of the inverse Sturm-Liouville problems which have various versions. Among them, the best known one is studied by Gel’fand and Levitan1, in which the potential function is uniquely determined by spectral function. McLaughlin2has given an analytical method to treat this type of inverse problems.

There were many works to develop algorithms for solving the inverse Sturm-Liouville problem of reconstructing potential function from eigenvalues3,4, which is known as the inverse spectral problem or inverse eigenvalue problem5. On the other hand, McLaughlin 6first noted that it is possible to obtain the potential function and boundary conditions using

only the set of nodal points. This interesting problem has soon been known as the inverse nodal problem7–10.

Numerical methods often transform the inverse Sturm-Liouville problem into an inverse eigenvalue problem of a certain matrix11. However, many of these discretizations into a ma- trix form have higher eigenvalues significantly differening from those of true eigenvalues. As a consequence, the inverse algorithms based on these discretizations require careful implemen- tation3,4.

In this study, the data of spectral function is chosen in order to identify a spatial- dependent potential function; hence, the present inverse Sturm-Liouville problem is less diffi- cult than those considered in3,4,6–10.

First, we transform the inverse Sturm-Liouville problem into a parameter identification problem governed by a parabolic type partial differential equationPDE. Then, a one-step group-preserving schemeGPSfor a semidiscretization of that PDE is established, which can be used to derive a closed-form solution of the estimated potential function at discretized spa- tial points. This type approach is first time appeared in the literature.

Let us consider a second-order ordinary differential equation ODE describing the Sturm-Liouville boundary value problem:

d dx

pxdy

dx

qx λrx

y Fx inx₀≤x≤x_f, 1.1

y A0 atx x0, 1.2

y B0 atx xf. 1.3

The direct problem for the given conditions in1.2and1.3and the given functionspx, qx,rx, andFxis to find the solutionyxof the second-order boundary value problem BVP. Specifically, whenFx 0, we have a Sturm-Liouville problem to determine the eigen- valueλand eigenfunctionyx.

The present inverse problem of Sturm-Liouville is to estimateqxby using the informa- tion of one eigenfunctionyxand its corresponding eigenvalueλ, and for the Sturm-Liouville differential operator is to estimatepxby using the data ofyxwhenqx rx 0.

For the case whenpxfis known andqx rx 0 in1.1, we propose a noniterative method to calculatepxat discretized spatial points. This problem could also be solved by the iterative method given by Keung and Zou12for the elliptic problem∇ ·p∇u F. Some of the numerical examples in Keung and Zou12involve Sturm-Liouville problems, but the method proposed here requires less computation for these problems.

For the case ofqx rx 0 from1.1, it follows directly that pxyx p

x0

y x0

_x

x0

Fsds. 1.4

Ifyx,px0andyx0are available, the above equation simply gives the unknown param- eterpxby dividing both sides by yx. However, becauseyx is usually not given in a closed-form and is given discretizedly under a perturbation by noise, we require a numerical technique to findyx. As mentioned by Li13, several techniques were developed to con- struct useful difference formulas for numerical derivativesNDs. In addition to the references in13, we also mention the book by Shu14. Among the many NDs, we only employ the

method by Ahn et al.15to compare it with our new method for numerical examples given inSection 6. Ahn et al.15have used a Volterra integral equation of the second kind to derive the following numerical derivative of a functionfxunder noise denoted byf_δx:

f_δx −1

α²exp −x α

_x

0

exp s α

f_δsdsf_δx

α , 1.5

whereαis a regularized parameter andf_δxis a numerical derivative off_δx.

Lie-group is a differentiable manifold, endowed a group structure that is compatible with the underlying topology of manifold. The main purpose of Lie-group solver is for pro- viding a better algorithm that retains the orbit generated from numerical solution on the man- ifold which associated with the Lie-group16,17. The retention of Lie-group structure under discretization is vital in the recovery of qualitatively correct behavior in the minimization of numerical error18,19.

Liu20has extended the GPS developed in19for ODEs to solve the BVPs, and the numerical results reveal that the GPS is a rather promising method to effectively solve the two- point BVPs. In that construction of Lie-group method for the calculations of BVPs, Liu20has introduced the idea of one-step GPS by utilizing the closure property of Lie-group, and hence, the new shooting method has been named the Lie-group shooting method.

It should be stressed that the one-step property of Lie-group is usually not shared by other numerical methods because those methods do not belong to the Lie-group type. This important property has been used by Liu21to establish a one-step estimation method to es- timate the temperature-dependent heat conductivity, and then extended to estimate heat con- ductivity and heat capacity 22–24. Its group structure gives the Lie-group method a great advantage over other numerical methods. It is a powerful technique to solve the inverse prob- lem of parameter identification.

This paper is arranged as follows. We introduce a novel approach of an inverse Sturm- Liouville problem inSection 2by transforming it into an identification problem of a parabolic type PDE, and then discretizing the PDE into a system of ODEs at discretized spatial points.

InSection 3, we give a brief sketch of the GPS for ODEs for a self-content reason. Due to its good property of Lie-group, we will propose a one-step GPS which can be used to identify the parameters appeared in the PDE. The resulting algebraic equation is derived inSection 4when we apply the one-step GPS to identifyqx. We demonstrate how the Lie-group theory can help us to solve the parameter estimation equation in a closed-form. InSection 5, we turn our attention to the estimation ofpxwhich leads again to a closed-form solution of the parameter pxat discretized spatial points. InSection 6, several numerical examples are examined to test the Lie-group estimation methodLGEM. Finally, we give conclusions inSection 7.

2. A novel approach

2.1. Transformation into a PDE

In the solution of linear PDE, a common technique is the separation of variables from which the PDE is transformed into ODEs. We may reverse this process by considering

ux, t 1tyx, 2.1

such that1.1–1.3are changed to

∂ux, t

∂t

∂

∂x

px∂ux, t

∂x

qx λrx

ux, t hx, t inx0≤x≤xf, 0< t≤T,

2.2 u

x0, t

A01t, 2.3

u x_f, t

B₀1t, 2.4

ux,0 yx, 2.5

wherehx, t yx−1tFx, and the last initial condition follows from2.1directly.

Equation2.2is a heat conduction equation, where we are attempting to estimatepx orqxunder a given sourcehx, t.

2.2. Semidiscretization

The semidiscrete procedure of PDE produces a coupled system of ODEs. For the one- dimensional heat conduction2.2, we adopt the numerical method of line to discretize the spatial coordinatexby

∂ux, t

∂x _{x x}

i x0iΔx

ui1t−uit

Δx , 2.6

∂²ux, t

∂x² _{x x}

i x0iΔx

ui1t−2uit ui−1t

Δx² , 2.7 whereΔx xf−x0/n1is a uniform discretization spacing length, anduit ux0iΔx, t for a simple notation. Such that2.2can be approximated by

˙

uit pi

Δx²

ui1t−2uit ui−1t

p_iui1t−uit Δx

qiλri

uit hit, i 1, . . . , n,

2.8

wherepi pxi,p_i pxi,qi qxi,ri rxi, andhit yi−1tFiwithyi yxiand F_i Fxi.

Wheni 1, the termu₀tis determined by boundary condition2.3withu₀t A₀1 t. Similarly, wheni n, the termun1tis determined by boundary condition 2.4 with un1t B01t. The next step is to advance the solution from a given initial condition to a desired timeT. However,2.8has totallyncoupled linear ODEs for thenvariablesu_it, i 1, . . . , n, which can be numerically integrated to obtainuiT.

In this section, we have transformed the inverse Sturm-Liouville problem in1.1into an inverse parameter identified problem for the PDE in2.2, and finally to an estimation ofnco- efficientsq_iorp_iin then-dimensional linear ODEs system. The data required in the estimation are the discretization ofyxat discretized spatial points, that is,yi yxi.

3. GPS for differential equations system 3.1. Group-preserving scheme

Upon letting u u1, . . . , u_n^Tand denoting f the right-hand side of2.8, we can write it as a vector form:

˙u fu, t, u∈Rⁿ, t∈R. 3.1 Liu19has embedded3.1into an augmented dynamical system, which is concerned with not only the evolution of state variables but also the evolution of the magnitude of the state variables vector:

d dt

u u

⎡

⎢⎢

⎢⎣

0n×n fu, t u f^Tu, t

u 0

⎤

⎥⎥

⎥⎦ u

u

. 3.2

Equation 3.2gives us a Minkowskian structure of the augmented state variables of X : u^T,u^Tto satisfy the cone condition:

X^TgX 0, 3.3

where

g

I_n 0n×1

01×n −1

3.4 is a Minkowski metric, I_nis the identity matrix of ordern, and the superscript T stands for the transpose. In terms ofu,u,3.3becomes

X^TgX u·u− u² u²− u² 0, 3.5 where the dot between twon-dimensional vectors denotes their Euclidean inner product. The cone condition is thus the most natural constraint that we can impose on the dynamical system 3.2.

Consequently, we have ann1-dimensional augmented system:

X˙ AX 3.6

with a constraint3.3, where

A :

⎡

⎢⎢

⎢⎣

0n×n fu, t u f^Tu, t

u 0

⎤

⎥⎥

⎥⎦, 3.7

satisfying that

A^TggA 0 3.8

is a Lie algebrason,1of the proper orthochronous Lorentz groupSO_on,1.

Although the dimension of the new system is raised one more, it has been shown that the new system has an advantage to permit the group-preserving schemeGPSgiven as follows 19:

X1 GX, 3.9

G^TgG g, 3.10

det G 1, 3.11

G⁰₀>0, 3.12

whereG⁰₀is the 00th component of G, Xdenotes the numerical value of X at the discrete time t, and G∈SOon,1is the group value of G at a timet. If Gsatisfies the properties in 3.10–3.12, then Xsatisfies the cone condition in3.3.

The Lie-group G can be generated from A∈son,1by an exponential mapping,

G exp

ΔtA

⎡

⎢⎢

⎢⎢

⎣

Ina−1

f² ff^T bf f bf^T

f a

⎤

⎥⎥

⎥⎥

⎦, 3.13

where

a : cosh Δtf u

,

b : sinh Δtf u

.

3.14

Substituting3.13for Ginto3.9, we obtain

u1 uηf, 3.15

u₁ au b

ff·u, 3.16 where

η : buf a−1

f·u

f² 3.17 is an adaptive factor. From f·u≥ −fu, we can prove that

η≥

1−exp −Δtf

u u

f >0 ∀Δt >0. 3.18 This scheme is group properties preserved for allΔt >0.

3.2. One-step GPS

Applying scheme3.15on2.8, we can compute the heat conduction equation by the GPS.

Assume that the total timeTis divided byKsteps, that is, the time step size we use in the GPS isΔt T/K.

Starting from an initial augmented condition X0 X0, we may want to calculate the value XTat a desired timet T. By3.9, we can obtain that

X_T G_KΔt· · ·G₁ΔtX0, 3.19 where XTapproximates the real XTwithin a certain accuracy depending onΔt. However, let us recall that each Gi, i 1, . . . , K, is an element of the Lie-groupSOon,1, and by the closure property of Lie-group, G_KΔt· · ·G₁Δtis also a Lie-group denoted by GT. Hence, we have

XT GTX0. 3.20

This is a one-step transformation from X0to XT.

Usually, it is very hard to find an exact solution of GT; however, a numerical one may be obtained approximately without any difficulty. The most simple method to calculate GT is given by

GT

⎡

⎢⎢

⎢⎢

⎢⎣

I_na−1

f0² f₀f^T₀ bf₀ f₀ bf^T₀

f0 a

⎤

⎥⎥

⎥⎥

⎥⎦, 3.21

where

a: cosh Tf0 u₀

,

b: sinh Tf₀ u0

.

3.22

Then from3.15and3.16, we obtain a one-step GPS:

u_T u₀ηf₀, 3.23

uT au0bf0·u0

f₀ , 3.24 where

η a−1f0·u₀bu₀f₀

f₀² . 3.25

4. Identifyingqxby the LGEM

In this section, we will start to estimate the potential functionqx. By using the one-step GPS, we also suppose that the initial value ofux,0 yxis given and its corresponding eigen- value is known.

Applying the one-step GPS in3.23on2.8from timet 0 to timet T, we obtain a nonlinear equation forq_i:

u^T_i u⁰_i ηpi

Δx²

u⁰_i1−2u⁰_i u⁰_i−1

ηp_iu⁰_i1−u⁰_i Δx η

q_iλr_i

u⁰_i ηh_i0. 4.1 It is not difficult to rewrite4.1as

q_i 1 u⁰_i

u^T_i −u⁰_i η − pi

Δx²

u⁰_i1−2u⁰_i u⁰_i−1

− p_i Δx

u⁰_i1−u⁰_i

−λr_iu⁰_i −h_i0

, 4.2

noting thatηin the above is not a constant but a nonlinear function ofq_ias shown by3.25.

Therefore, in this stage, we cannot calculateqiby a simple equation. However, we will prove below thatηis fully determined byu⁰_i andu^T_i.

In order to solveq_i, let us return to3.23:

f₀ 1 η

u_T−u₀

. 4.3

Substituting it for f0into3.24, we obtain uT

u₀ a b

uT−u0

·u0

u_T−u₀u₀, 4.4 where

a: cosh Tu_T−u₀ ηu₀

, 4.5

b: sinh Tu_T−u₀ ηu0

. 4.6

Let

cosθ:

u_T−u₀

·u₀

uT−u0u0, 4.7 S: Tu_T−u₀

u0 , 4.8 and from4.4–4.6, it follows that

uT

u₀ cosh S η

cosθsinh S η

. 4.9

Upon defining

Z: exp S η

, 4.10

and from4.9, we obtain a quadratic equation forZ:

1cosθZ²−2u_T

u0 Z1−cosθ 0. 4.11 The solution is found to be

Z uT/u0±uT/u0²−

1−cos²θ

1cosθ if±cosθ >0; 4.12

and from4.10, we obtain a closed-form solution ofη:

η Tu_T−u₀

u0lnZ . 4.13 Up to here, we must point out that for a givenT,ηis fully determined by u0 and uT which are supposed to be known. Therefore, the original nonlinear equation4.2becomes a linear equation forq_i.

By using2.1, we have

u^T_i 1Tu⁰_i 1Tyi, 4.14

and thus the vector uTis proportional to u0with a multiplier 1T larger than 1. Under this condition, we have cosθ 1, andZis given by

Z u_T

u0 1T, 4.15 and hence from4.13, we have

η T²

ln1T. 4.16

Inserting4.14and4.16into4.2, we obtain a very simple formula to estimateq_iby qi 1

yi

yiln1T

T − pi

Δx²

yi1−2yiyi−1

− p_i Δx

yi1−yi

−λriyi−yiFi

. 4.17 This solution is in a closed-form forq_i.

In the above, we have mentioned thatηis a nonlinear function ofqi; however, by viewing 4.7,4.12, and4.13, it is known thatηis fully determined by u0and uT. Furthermore, by using4.14ηbecomes a constant given by4.16. This point is very important for our closed- form solution of parameter. The key points rely on the construction of the method by using the one-step GPS for the estimation of parameter, and the full use of then1 equations3.23and 3.24. To distinguish the present method by a joint use of the one-step GPS and the closed-form solution with the aid of3.24, we may call the new method a Lie-group estimation method LGEM.

5. Applying the LGEM to estimatepx

In this section, we will derive a simple linear equations system to solve the coefficientsp_i, i 1, . . . , n. However, for simplicity, we assume thatqx rx 0 in this section.

A similar finite difference as that used in2.6foruxcan be used forpxin2.8. In doing so, we can obtain a system of ODEs foruwithtas an independent variable:

˙

uit pi1−p_i Δx

ui1t−uit

Δx piui1t−2uit ui−1t

Δx² hit. 5.1

The known initial condition is given by u⁰_i y

xi

, i 1, . . . , n, 5.2

which is obtained from2.1by a discretization.

Applying the same idea of LGEM on5.1, we can obtain a closed-form formula to esti- matep_i:

pi Δx² u⁰_i −u⁰_i−1

u⁰_i1−u⁰_i

Δx² pi1hi0−1 η

u^T_i −u⁰_i

, 5.3

and, moreover, by using the data ofu^T_i given by4.14and4.16forη, we can derive a much simple equation forp_i:

p_i Δx² y_i−y_i−1

y_i1−y_i

Δx² pi1y_i−F_i−y_iln1T T

. 5.4

This will be called a closed-form estimation method. The above equation can be used sequen- tially to findpi, i n, . . . ,1,if we knowpn1a priori. Here,pn1is the right-end boundary value ofpx, and is supposed to be known for simplicity.

However, we can develop another estimation method through iterations. The numer- ical procedures for estimating pi are described as follows. We assume an initial value ofpi, for example, pi 1. Substituting it into 5.1, we can apply the GPS to integrate it from t 0 to t T through T/Δt steps. Then, we obtain u^T_i. Inserting it into 5.3, we can calculate a new pi, which is then compared with the old pi. If the difference of these two sets ofpi is smaller than a given criterion, then we stop the iteration and the finalpi is ob- tained.

The processes are summarized as follows:

igive an initialp_i 1;

iiforj 1,2. . ., we repeat the following calculations; calculateu^T_i by using the GPS in 3.15to integrate5.1fromt 0 tot T, where f is a vector form of the right-hand side of5.1;

iiiinsert the above calculatedu^T_i denoted byu^T_ijtogether withu⁰_i given by5.2into

p^j_i Δx² u⁰_i −u⁰_i−1

u⁰_i1−u⁰_i

Δx² p^j_i1− 1 η^j

u^T_ij−u⁰_i h_i0

, 5.5

where η^j is calculated from4.13by substituting u^j_T for each step. Ifp^j_i converges according to a given stopping criterion:

n i 1

p^j1_i −p^j_i²

< ε², 5.6

then stop, otherwise, go to stepii.

Basically, the present method has repeatedly used the time direction integration of5.1 to obtain the final time data to adjustpi, which will be called an iterative estimation method.

6. Numerical examples

Example 6.1. For a first example, we consider an inverse Sturm-Liouville problem to identify the potential function in the Schr ¨odinger equation:

yx

β−α²x²

yx 0,

y−∞ y∞ 0. 6.1

Here,β 2k1αis the eigenvalue, and

y_kx H_kxexp −αx² 2

6.2

is the eigenfunction, where the Hermite polynomials fork 0,1,2,3,4 are given by H₀x 1,

H₁x 2x, H₂x −24x², H₃x −12x8x³, H4x 12−48x²16x⁴.

6.3

In general,H_kx −1^ke^x2d^ke^−x2/dx^k.

In order to recover the potential functionqxfrom the given eigenvalue and eigenfunc- tion, we apply4.17on this problem by takingp 1,p 0,r 1, andF 0. We also take x₀ −5 andx_f 5 and letΔx 10/300,α 1, andT 0.0001.

The estimation errors ofqxare shown in Figures1aand1b, fork 1 andy1x andk 3 and y3xas the inputs on4.17. From these two figures, it can be seen that the estimations ofqxare quite accurate. However, near the boundaries, the errors are increased.

In order to avoid the boundary effect on the estimation ofqx, we can extend the range ofx into a larger one.

4E-2

2E-2

0E0

−5 0 5

x

Estimationerrorofqx

a

6E-2 5E-2 4E-2 3E-2 2E-2 1E-2 0E0

−5 0 5

x

Estimationerrorofqx

b

Figure 1:Example 6.1of identifying the Schr ¨odinger equationaplotting the estimation error by using the second spectral function, whilebby using the fourth spectral function.

Example 6.2. For a second test example, we consider the Sturm-Liouville problem with x⁻¹yx λ1x⁻³yx 0,

y1 ye 0. 6.4

The eigenvalue isλk k²π², k ∈N, and the eigenfunction isyk axsinkπlnx, whereais an arbitrary nonzero constant fixed to bea 100.

InFigure 2, we compare the estimation errors by considering the noisy data withy_σ axsinπlnx σRi, whereRiare random numbers between−1 and1. TheL² error for σ 0 is about 0.327 while that forσ 0.01 is about 0.373.

Example 6.3. In this example, we estimatepxby settingqx 0 andrx 0. Let us use the following example to demonstrate the process inSection 5. This example is given by

px x−3²,

Fx 6x−3², hx, t yx−1tFx −x−3²56t. 6.5 Under the boundary conditions

u0, t 91t, u1, t 41t, 6.6

and the initial condition

ux,0 x−3², 6.7

0.25

0.2

0.15

0.1

0.05

0

0.8 1.2 1.6 2 2.4 2.8

x

Estimationerrorofqx

Error withσ 0 Error withσ 0.01

Figure 2:Example 6.2of identifyingqx, the estimation errors were plotted.

the exact solution is given by

ux, t x−3²1t. 6.8

We apply the LGEM on this identification ofpx, where we have fixed Δx 1/30, Δt 5×10⁻⁵, andT 0.01. Under the stopping criterion with 10⁻³, the process is conver- gent within 34 iterations. InFigure 3a, we plot the tentativep_ifor the first iteration, the fifth iteration, the tenth iteration, and the fifteenth iteration, the last of which is already close to the exact solution. The numerical solutions ofpiare close to the exact ones with theL²error about 0.0156, and the maximum relative error about 4×10⁻³as shown inFigure 3b.

Example 6.4. The following example has been calculated by Keung and Zou12 using the augmented Lagrange method:

px 32x²−2 sin2πx, Fx 2π

4x−4πcos2πx

cos2πx−4π²

32x²−2 sin2πx

sin2πx. 6.9

The exact solution is given by

ux, t 1tsin2πx. 6.10

We first apply the iterative LGEM on this identification ofpx, where we have fixed Δx 2/60,Δt 10⁻⁵, andT 0.003. Under the stopping criterion with 10⁻³, the process is convergent within 37 iterations. InFigure 4, we compare the estimated solution with exact solution. For this estimation, we have a maximum error with 0.377 and anL²error with 1.62.

9 8 7 6 5 4 3

15 10

5 1

0 0.2 0.4 0.6 0.8 1

Estimated Exact

x

px

a 4E-3

2E-3

0E0

x

0 0.2 0.4 0.6 0.8 1

Relativeerrorofpx

b

Figure 3:Example 6.3by using an iterative method:acomparing estimated and exactpx, andbplot- ting the relative error of estimation.

On the other hand, we also apply5.4on this estimation by usingΔx 2.2/150 andT 0.001, of which the result is shown inFigure 4by the dashed-dotted line. For the later estimation, we have a maximum error with 0.59 and anL²error with 1.402. It is slightly more accurate than the iterative method; however, there appear four kinks near the extremal points ofpx, where px 0.

Example 6.5. For this example, we consider a simple exactyx sinπxunder a noise given byyδ sinπx δcos3πx. The following results are used:

px x,

Fx πcosπx−π²xsinπx. 6.11

For the integration in1.5, we have employed the trapezoidal rule to calculatey_δxat dis- cretized spatial points. Then inserting them into1.4, we can obtainpxby the ND method underδ 0 andδ 0.02, whose errors are shown inFigure 5a. The best parameterαis fixed to be 0.95 forδ 0 and 0.85 forδ 0.02, not 0.0085 as that used by Ahn et al.15.

On the other hand, we also apply5.4on this estimation ofpxby usingΔx 0.01 and T 0.01, of which the errors underδ 0 andδ 0.02 are shown inFigure 5a. It can be seen that the maximum errors of LGEM are much smaller than that obtained by the ND method.

Especially, whenδ 0 the LGEM with all its absolute error is smaller than that by the ND method.

12 10 8 6 4 2

0

0 0.4 0.8 1.2 1.6 2

Exact

Estimated by closed-form method Estimated by iterative method

x

px

Figure 4: InExample 6.4we compare exactpxwith the estimated ones by using an iterative method and a closed-form method.

For both methods applied in this example, we also consider the random noise distur- bance given by y_σ sinπxi σRi, where Ri are random numbers between −1 and 1. In Figure 5b, we compare the numerical errors by applying the LGEM and ND for this case under a noise with σ 0.001. For the ND, the best parameter of α is α 0.5.

Obviously, the absolute error of LGEM is smaller than that given by the ND method. This also shows that the method of LGEM can be against the random and also definite distur- bances.

Example 6.6. Let us consider the following linear BVP:

xyx yx Fx, y1 0, y2 2, 6.12

whereFx x^3/2cosxhas no closed-form integral of_x

1Fsds, such thatyxhas also no closed-form solution. In this case, we applied the Lie-group shooting method 20to calcu- latey.

In 1.4, we require to know both px0 and yx0. This is a great drawback of the ND method. Because the calculated data ofyare discretized, we may approximateyx0by y1−y₀/Δx. Then, inserting the calculated data ofyinto1.4, we can obtainpxby the ND method underσ 0 andσ 0.001, whose errors, as shown inFigure 6, are very large. The best parameterαis fixed to be 0.45 forσ 0 and 0.43 forσ 0.001.

Then, we apply5.4on this estimation ofpxby usingΔx 0.01 andT 0.01 of which the errors underσ 0 andσ 0.02 are shown inFigure 6. The errors of LGEM are much smaller than that obtained by the ND method with a ratio 10⁻⁴of these two maximum errors for the un-noised case and 6.7×10⁻³for the noised case.

1E0

1E-1

1E-2

1E-3

0 0.1 0.2 0.3 0.4 0.5

x Error of ND withδ 0.02 Error of ND withδ 0 Error of LGEM withδ 0.02 Error of LGEM withδ 0

Errors

a

1E-1

1E-2

1E-3

1E-4

0 0.1 0.2 0.3 0.4 0.5

x Error of ND withσ 0.001 Error of LGEM withσ 0.001

Errors

b

Figure 5:Example 6.5comparing the estimation errors ofpx:adefinite disturbances withδ 0, 0.02, andbrandom disturbance withσ 0.001.

7. Conclusions

In order to estimate the potential function under a given spectral function and its correspond- ing eigenvalue, we have employed the LGEM to derive an algebraic equation and solved it in a closed form. We transformed the inverse Sturm-Liouville problem into a parameter identifi- cation problem for a parabolic type PDE, and then established a one-step GPS for the semidis- cretization of that PDE. We also established an iterative method to estimate the unknown coef- ficient in a second-order Sturm-Liouville operator.

1 1.2 1.4 1.6 1.8 2 1E2

1E1 1E0 1E-1

1E-2

1E-3

1E-4

x Error of ND withσ 0.001 Error of ND withσ 0 Error of LGEM withσ 0.001 Error of LGEM withσ 0

Errors

Figure 6:Example 6.6comparing the estimation errors ofpxunder random disturbances withσ 0 and σ 0.001.

Numerical examples were worked out, which show that the new LGEM is applicable for the estimations of unknown functions. When disturbances are exerted on the input data, we also verified that the present approach can be against them very well. The tested case shows that the LGEM is superior than the ND method examined here. Through this study, it can be concluded that the new estimation method is accurate, effective, and stable. Its numerical implementation is very simple and the computational speed is very fast.

References

1 I. M. Gel’fand and B. M. Levitan, “On the determination of a differential equation from its spectral function,” American Mathematical Society Translations, vol. 1, pp. 253–304, 1955.

2 J. R. McLaughlin, “Analytical methods for recovering coefficients in differential equations from spec- tral data,” SIAM Review, vol. 28, no. 1, pp. 53–72, 1986.

3 A. L. Andrew, “Numerical solution of inverse Sturm-Liouville problems,” The Australian & New Zealand Industrial and Applied Mathematics Journal, vol. 45, pp. C326–C337, 2004.

4 J. Paine, “A numerical method for the inverse Sturm-Liouville problem,” SIAM Journal on Scientific and Statistical Computing, vol. 5, no. 1, pp. 149–156, 1984.

5 M. T. Chu, “Inverse eigenvalue problems,” SIAM Review, vol. 40, no. 1, pp. 1–39, 1998.

6 J. R. McLaughlin, “Inverse spectral theory using nodal points as data—a uniqueness result,” Journal of Differential Equations, vol. 73, no. 2, pp. 354–362, 1988.

7 P. J. Browne and B. D. Sleeman, “Inverse nodal problems for Sturm-Liouville equations with eigenpa- rameter dependent boundary conditions,” Inverse Problems, vol. 12, no. 4, pp. 377–381, 1996.

8 Y.-H. Cheng, C. K. Law, and J. Tsay, “Remarks on a new inverse nodal problem,” Journal of Mathematical Analysis and Applications, vol. 248, no. 1, pp. 145–155, 2000.

9 O. H. Hald and J. R. McLaughlin, “Solutions of inverse nodal problems,” Inverse Problems, vol. 5, no. 3, pp. 307–347, 1989.

10 X.-F. Yang, “A solution of the inverse nodal problem,” Inverse Problems, vol. 13, no. 1, pp. 203–213, 1997.

11 D. Boley and G. H. Golub, “A survey of matrix inverse eigenvalue problems,” Inverse Problems, vol. 3, no. 4, pp. 595–622, 1987.

12 Y. L. Keung and J. Zou, “An efficient linear solver for nonlinear parameter identification problems,”

SIAM Journal on Scientific Computing, vol. 22, no. 5, pp. 1511–1526, 2000.

13 J. Li, “General explicit difference formulas for numerical differentiation,” Journal of Computational and Applied Mathematics, vol. 183, no. 1, pp. 29–52, 2005.

14 C. Shu, Differential Quadrature and Its Application in Engineering, Springer, London, UK, 2000.

15 S. Ahn, U. J. Choi, and A. G. Ramm, “A scheme for stable numerical differentiation,” Journal of Com- putational and Applied Mathematics, vol. 186, no. 2, pp. 325–334, 2006.

16 A. Iserles and A. Zanna, “Preserving algebraic invariants with Runge-Kutta methods,” Journal of Com- putational and Applied Mathematics, vol. 125, no. 1-2, pp. 69–81, 2000.

17 C.-S. Liu, “New integrating methods for time-varying linear systems and Lie-group computations,”

Computer Modeling in Engineering & Sciences, vol. 20, no. 3, pp. 157–175, 2007.

18 A. Iserles, H. Z. Munthe-Kaas, S. P. Nørsett, and A. Zanna, “Lie-group methods,” Acta Numerica, vol. 9, pp. 215–365, 2000.

19 C.-S. Liu, “Cone of non-linear dynamical system and group preserving schemes,” International Journal of Non-Linear Mechanics, vol. 36, no. 7, pp. 1047–1068, 2001.

20 C.-S. Liu, “The Lie-group shooting method for nonlinear two-point boundary value problems exhibit- ing multiple solutions,” Computer Modeling in Engineering & Sciences, vol. 13, no. 2, pp. 149–163, 2006.

21 C.-S. Liu, “One-step GPS for the estimation of temperature-dependent thermal conductivity,” Interna- tional Journal of Heat and Mass Transfer, vol. 49, no. 17-18, pp. 3084–3093, 2006.

22 C.-S. Liu, “An efficient simultaneous estimation of temperature-dependent thermophysical proper- ties,” Computer Modeling in Engineering & Sciences, vol. 14, no. 2, pp. 77–90, 2006.

23 C.-S. Liu, “Identification of temperature-dependent thermophysical properties in a partial differential equation subject to extra final measurement data,” Numerical Methods for Partial Differential Equations, vol. 23, no. 5, pp. 1083–1109, 2007.

24 C.-S. Liu, L.-W. Liu, and H.-K. Hong, “Highly accurate computation of spatial-dependent heat con- ductivity and heat capacity in inverse thermal problem,” Computer Modeling in Engineering & Sciences, vol. 17, no. 1, pp. 1–18, 2007.