15th Annual Conference on Applied Mathematics, Univ. of Central Oklahoma, Electronic Journal of Differential Equations, Conference 02, 1999, pp. 125–131. ISSN: 1072-6691. URL: http://ejde.math.swt.edu or http://ejde.math.unt.edu ftp ejde.math.swt.edu (login: ftp)

Operational WKB solution to the initial/final-value problem for Beechem-Haas

equations ^{∗}

Valentino Anthony Simpao

Abstract

Operational and Wantzel-Kramers-Brillouin (WKB) methods are utilized to obtain an analytical approximate solution of the Beechem-Haas equation for the initial/final-value problem with system parameters of radial profile. The Beechem-Haas equation herein is expanded from its traditional form, which determines its solutions in a spherical coordinate system as functions of the radial and time coordinates only, to one which includes the angular coordinates as well, thus providing a natural modeling framework for polarization phenomena.

1 Introduction

Consider the Beechem-Haas equation and its initial/boundary conditions, as given in [1],

∂_{t}\overline{N}(r, t) = ( ∑_{j=1}^{n} \frac{α_{j}}{τ_{j}} (1 + (R_{0}/r)^{6}))\overline{N}(r, t) + \frac{1}{N_{0}(r)}∂_{r}[N_{0}(r)D(r)∂_{r}\overline{N}(r, t)] ,

\overline{N}(r, t)|_{t=0} = 1 , \overline{N}_{t=0}(r, t) = 1 ,  (1)

∂_{r}\overline{N}(r, t)|_{r=rmin} = 0 , ∂_{r}\overline{N}(r, t)|_{r=rmax} = 0 ,

where 0 < r_{min} ≤ r ≤ r_{max}, t ≥ 0, and α_{j}, τ_{j}, R_{0} are positive real numbers. Unless otherwise directed, all variables, parameters, and domains are real.

In the present work, the full 3-dimensional form of the spatial derivatives in spherical coordinates shall be applied to (1), resulting in an enhanced form of the Beechem Haas equation, which will allow for the inclusion of angular components to the solution in a natural way. Particularly, this allows (1) to be extended to include the nontrivial effect of angular coupling to the radial component of the density function (i.e., polarization phenomena). Finally, the solution of the initial/final-value problem for the enhanced version of the Beechem-Haas equation will be addressed herein. This is basically the “on/off” control problem in the time-domain for the system. As a result of solving this problem, one then determines the system function necessary to produce a required output state, given a prescribed input state. Full technical details of the terms in (1), other than those supplied herein, are unnecessary for the present work. Therefore, the interested reader is encouraged to examine the original source [1]. The primary purpose of the present work is to simply supply a derivation of the results named above; applications and deeper studies shall appear elsewhere.

∗1991 Mathematics Subject Classifications: 92C45, 34E20, 58J35.

Key words and phrases: operational methods, WKB, Analytic solution.

\begin{matrix}c \\ ◯ \end{matrix}2000 Southwest Texas State University and University of North Texas.

Published January 21, 2000.

2 Main result

With these notions in place, the enhanced Beechem-Haas equation equation with attendant generalized auxiliary conditions becomes

∂_{t}\overline{N}_{E}(r, θ, ϕ, t) = ( ∑_{j=1}^{n} \frac{α_{j}}{τ_{j}} (1 + (R_{0}/r)^{6}))\overline{N}_{E}(r, θ, ϕ, t)
+\frac{4πr^{2}}{N_{0}(r)}∇ · [\frac{N_{0}(r)}{4πr^{2}}D(r)∇\overline{N}_{E}(r, θ, ϕ, t)] ,  (2)
\overline{N}_{E}(r, θ, ϕ, t)|_{t=0} = \overline{N}_{E_{t=0}}(r) , \overline{N}_{E}(r, θ, ϕ, t)|_{t=T} = \overline{N}_{E_{t=T}} (r)
∂_{r}[\overline{N}_{E}(r, θ, ϕ, t)|_{r=rmin} = 0 , ∂_{r}[\overline{N}_{E}(r, θ, ϕ, t)|_{r=rmax} = 0 ,

where 0 < r_{min} ≤ r ≤ r_{max}, 0 ≤ θ ≤ π, 0 ≤ ϕ ≤ 2π, t ≥ 0, and α_{j}, τ_{j}, R_{0} are positive real numbers. Expanding the above expression,

\frac{1}{D(r)}∂_{t}\overline{N}_{E}(r, θ, ϕ, t)

= \frac{1}{D(r)}( ∑_{j=1}^{n} \frac{α_{j}}{τ_{j}} (1 + (R_{0}/r)^{6}))\overline{N}_{E}(r, θ, ϕ, t)  (3)

+∇[ln (\frac{N_{0}(r)D(r)}{4πr^{2}})] · ∇\overline{N}_{E}(r, θ, ϕ, t) + Δ\overline{N}_{E}(r, θ, ϕ, t) .

Recall that the finite Fourier transform and its inverse [2] are defined as

f̂(n) = F_{T} [f(t)] = ∫_{0}^{T} f(t)e^{−inωt} dt ,

f(t) = F_{T}^{−1}[f̂(n)] = \frac{1}{T} (f̂(0) + ∑_{n=1}^{∞}[f̂(n)e^{inωt} + f̂(−n)e^{−inωt}])

Applying the finite Fourier transform on the time domain of (3), we obtain

\frac{1}{D(r)}[\overline{N}_{E_{t=T}} (r) − \overline{N}_{E_{t=0}}(r)]

= \frac{1}{D(r)}( ∑_{j=1}^{n} \frac{α_{j}}{τ_{j}} (1 + (R_{0}/r)^{6}) − inω)N̂_{E}(r, θ, ϕ, n)  (4)

+∇[ln (\frac{N_{0}(r)D(r)}{4πr^{2}})] · ∇N̂_{E}(r, θ, ϕ, n) + ΔN̂_{E}(r, θ, ϕ, n) .

Reducing this to canonical form,

\frac{1}{D(r)}[\overline{N}_{E_{t=T}} (r) − \overline{N}_{E_{t=0}}(r)]\sqrt{\frac{4πr^{2}}{N_{0}(r)D(r)}}

= ΔÑ_{E}(r, θ, ϕ, n) + (\frac{1}{D(r)}( ∑_{j=1}^{n} \frac{α_{j}}{τ_{j}} (1 + (R_{0}/r)^{6}) − inω)  (5)

−\frac{1}{2}Δ[ln (\frac{N_{0}(r)D(r)}{4πr^{2}})] − \frac{1}{4}|∇[ln (\frac{N_{0}(r)D(r)}{4πr^{2}})]|^{2})Ñ_{E}(r, θ, ϕ, n) ,

where

Ñ_{E}(r, θ, ϕ, n) = N̂_{E}(r, θ, ϕ, n)\sqrt{\frac{N_{0}(r)D(r)}{4πr^{2}}} ,

∂_{r}[Ñ_{E}(r, θ, ϕ, n)\sqrt{\frac{4πr^{2}}{N_{0}(r)D(r)}}]|_{r=rmin} = 0 ,

∂_{r}[Ñ_{E}(r, θ, ϕ, n)\sqrt{\frac{4πr^{2}}{N_{0}(r)D(r)}}]|_{r=rmax} = 0 .

Now we consider the homogeneous equation corresponding to (5),

0 = ΔÑ_{E_{h}}(r, θ, ϕ, n) + (\frac{1}{D(r)}( ∑_{j=1}^{n} \frac{α_{j}}{τ_{j}} (1 + (R_{0}/r)^{6}) − inω)  (6)

−\frac{1}{2}Δ[ln (\frac{N_{0}(r)D(r)}{4πr^{2}})] − \frac{1}{4}|∇[ln (\frac{N_{0}(r)D(r)}{4πr^{2}})]|^{2})Ñ_{E_{h}}(r, θ, ϕ, n) .

Let Ñ_{E_{h}}(r, θ, ϕ, n) = R_{c}(r, n)Θ(θ)Φ(ϕ). Then via separation of variables,

\frac{1}{r^{2}R_{c}(r, n)}D_{r}(r^{2}D_{r}R_{c}(r, n))

+\frac{1}{r^{2} sin(θ)Θ(θ)}D_{θ}(sin(θ)D_{θ}Θ(θ)) + \frac{1}{r^{2} sin^{2}(θ)Φ(θ)}D_{θ}^{2}Φ(θ)  (7)
= \frac{1}{D(r)}( ∑_{j=1}^{n} \frac{α_{j}}{τ_{j}} (1 + (R_{0}/r)^{6}) − inω)

−\frac{1}{2}Δ[ln (\frac{N_{0}(r)D(r)}{4πr^{2}})] − \frac{1}{4}|∇[ln (\frac{N_{0}(r)D(r)}{4πr^{2}})]|^{2} .

Hence

D_{r}^{2}(rR_{c}(r, n)) − (\frac{λ}{r^{2}} + \frac{1}{D(r)}( ∑_{j=1}^{n} \frac{α_{j}}{τ_{j}} (1 + (R_{0}/r)^{6}) − inω)  (8)

−\frac{1}{2}Δ[ln (\frac{N_{0}(r)D(r)}{4πr^{2}})] − \frac{1}{4}|∇[ln (\frac{N_{0}(r)D(r)}{4πr^{2}})]|^{2})rR_{c}(r, n) = 0 ,

\frac{1}{Φ(ϕ)}D_{ϕ}^{2}Φ(ϕ) = −m^{2} ,

\frac{1}{sin θ} D_{θ}(sin(θ)D_{θ}Θ(θ)) − \frac{m^{2}}{sin^{2} θ} Θ(θ) + λΘ(θ) = 0 ,

where herein λ, m are arbitrary real parameters and n is a positive integer.

Regarding the solutions of the above angular component equations, they are most generally expressed as

Θ(θ)Φ(ϕ) = [c_{ϕ_{1}}e^{imϕ} + c_{ϕ_{2}}e^{−imϕ}]  (9)

×(c_{θ_{1}} \frac{1}{Γ(1 − m)}(\frac{1 + cos θ}{1 − cos θ})^{m/2} + c_{θ_{2}} \frac{π cot(mπ)}{2Γ(1 − m)}(\frac{1 + cos θ}{1 − cos θ})^{m/2}

−c_{θ_{2}} \frac{πΓ(m + λ + 1)}{2 sin(mπ)Γ(λ − m + 1)Γ(1 + m)}(\frac{1 + cos θ}{1 − cos θ})^{−m/2})

×F_{2,1}(−λ, λ; 1 − m; \frac{1 − cos θ}{2} ) ,

where c_{ϕ_{1}}, c_{ϕ_{2}}, c_{θ_{1}}, c_{θ_{2}} are arbitrary constants. From (8) it is clear that the radial ordinary differential equation couples to the angular ordinary differential equation only when λ is non-zero. When λ ≠ 0, (9) can take different forms, depending on the values of λ and m, which may be arbitrary real numbers. Indeed for λ and m integers, the angular dependence may be specified in terms of a finite superposition of spherical harmonics Y_{λ}^{m}(θ, ϕ), rather than via the Gauss hyper-geometric functions F_{2,1} used in in (9).

However, even if m is an integer and λ is not so restricted, the angular solutions will still [in general, by necessity] require the Gauss hyper-geometric formulation [3]. With an eye toward future results, the phenomena modeled by (8) may or may not require λ and m to be integer, depending upon the possible presence of more general non-homogeneous auxiliary conditions. For economy of notation result (9) may be denoted as Λ(θ, λ : ϕ, m) in general. However, since the present homogeneous auxiliary conditions on the solutions are just the standard ones for boundedness and single-valuedness, it is here sufficient to particularly designate the angular components as spherical harmonics Y_{λ}^{m}(θ, ϕ) with λ and m integers. Hence

Ñ_{Eh}(r, θ, ϕ, n) = R_{c}(r, n)Θ(θ)Φ(ϕ) = R_{c}(r, n)Λ(θ, λ : ϕ, m) = R_{c}(r, n)Y_{λ}^{m}(θ, ϕ) .

The bulk of the remaining analysis concerns the radial equation (8). Now applying the WKB admissibility criterion [4] to (8),

Φ(r, n) = −\frac{λ}{r^{2}} − \frac{1}{D(r)}( ∑_{j=1}^{n} \frac{α_{j}}{τ_{j}} (1 + (R_{0}/r)^{6}) − inω)  (11)

+\frac{1}{2}Δ[ln (\frac{N_{0}(r)D(r)}{4πr^{2}})] + \frac{1}{4}|∇[ln (\frac{N_{0}(r)D(r)}{4πr^{2}})]|^{2} ,

\frac{|D_{r}Φ(r, n)/2|}{|Φ(r, n)|^{3/2}} ≪ 1 .

Hence the WKB results are accurate when the above criterion is met by the coefficient over the given (r, n)-domain. Recalling (10), there obtains the linearly independent WKB analytical solutions for the homogeneous radial equation in (8):

R_{c_{1}}(r, n) ≈ R_{c_{1}WKB}(r, n)
= \{ \begin{matrix}(−Φ(r, n))^{−1/4} cos(∫ (−Φ(r, n))^{1/2} dr) & if Φ(r, n) < 0, \\ (Φ(r, n))^{−1/4} exp(∫ (Φ(r, n))^{1/2} dr) & if 0 < Φ(r, n). \end{matrix}
R_{c_{2}}(r, n) ≈ R_{c_{2}WKB}(r, n)  (12)
= \{ \begin{matrix}(−Φ(r, n))^{−1/4} sin(∫ (−Φ(r, n))^{1/2} dr) & if Φ(r, n) < 0, \\ (Φ(r, n))^{−1/4} exp(− ∫ (Φ(r, n))^{1/2} dr) & if 0 < Φ(r, n). \end{matrix}

Hence, the general solution of the radial part of the full non-homogeneous equation (5) may be supplied via (12) as

c_{1}(n)R_{c_{1}}(r, n) + c_{2}(n)R_{c_{2}}(r, n)  (13)

+R_{c_{2}}(r, n) ∫ \frac{\overline{N}_{E_{t=T}} (r) − \overline{N}_{E_{t=0}}(r)}{D(r)w̃(r, n)}R_{c_{1}}(r, n)\sqrt{\frac{4πr^{2}}{N_{0}(r)D(r)}} dr

−R_{c_{1}}(r, n) ∫ \frac{\overline{N}_{E_{t=T}} (r) − \overline{N}_{E_{t=0}}(r)}{D(r)w̃(r, n)}R_{c_{2}}(r, n)\sqrt{\frac{4πr^{2}}{N_{0}(r)D(r)}} dr ,

where

w̃(r, n) = R_{c_{1}}(r, n)D_{r}R_{c_{2}}(r, n) − R_{c_{2}}(r, n)D_{r}R_{c_{1}}(r, n) .

By ansatz (10), (13), and (5), the time-transformed version (4) of equation (12) has general solution

N̂_{E}(r, θ, ϕ, n) = \sqrt{\frac{4πr^{2}}{N_{0}(r)D(r)}}Ñ_{E}(r, θ, ϕ, n)

= \sqrt{\frac{4πr^{2}}{N_{0}(r)D(r)}}Y_{λ}^{m}(θϕ)(c_{1}(n)R_{c_{1}}(r, n) + c_{2}(n)R_{c_{2}}(r, n)  (14)

+R_{c_{2}}(r, n) ∫ \frac{\overline{N}_{E_{t=T}} (r) − \overline{N}_{E_{t=0}}(r)}{D(r)w̃(r, n)}R_{c_{1}}(r, n)\sqrt{\frac{4πr^{2}}{N_{0}(r)D(r)}} dr

−R_{c_{1}}(r, n) ∫ \frac{\overline{N}_{E_{t=T}} (r) − \overline{N}_{E_{t=0}}(r)}{D(r)w̃(r, n)}R_{c_{2}}(r, n)\sqrt{\frac{4πr^{2}}{N_{0}(r)D(r)}} dr) ,

Now consider the boundary conditions of (4), which are the time-transformed boundary conditions of (1),

∂_{r}N̂_{E}(r, θ, ϕ, n)|_{r=rmin} = 0 , ∂_{r}N̂_{E}(r, θ, ϕ, n)|_{r=rmax} = 0 .  (15)

Without loss of generality, we may ignore the angular factor Y_{λ}^{m}(θϕ) and consider only the radial terms in (14). Given the particular boundary conditions (15), the formula (14) may be used to calculate the appropriate particular solution matching the boundary conditions in the (r, n)-domain. In general, there will result a 2x2 matrix equation for the coefficients c_{1}(n) and c_{2}(n).

( \begin{matrix}P_{11} & P_{21} \\ P_{21} & P_{22} \end{matrix} ) ( \begin{matrix}c_{1}(n) \\ c_{2}(n) \end{matrix} ) = ( \begin{matrix}Q_{1} \\ Q_{2} \end{matrix} ) ,  (16)

where

P_{11} = ∂_{t}(\sqrt{\frac{4πr^{2}}{N_{0}(r)D(r)}}R_{c_{1}}(r, n))|_{r=rmin} ,

P_{12} = ∂_{t}(\sqrt{\frac{4πr^{2}}{N_{0}(r)D(r)}}R_{c_{2}}(r, n))|_{r=rmin} ,

P_{21} = ∂_{t}(\sqrt{\frac{4πr^{2}}{N_{0}(r)D(r)}}R_{c_{1}}(r, n))|_{r=rmax} ,

P_{22} = ∂_{t}(\sqrt{\frac{4πr^{2}}{N_{0}(r)D(r)}}R_{c_{2}}(r, n))|_{r=rmax} ,

Q_{1} = −∂_{r}(R_{c_{2}}(r, n) ∫ \frac{\overline{N}_{E_{t=T}} (r) − \overline{N}_{E_{t=0}}(r)}{D(r)w̃(r, n)}R_{c_{1}}(r, n)\sqrt{\frac{4πr^{2}}{N_{0}(r)D(r)}} dr

−R_{c_{1}}(r, n) ∫ \frac{\overline{N}_{E_{t=T}} (r) − \overline{N}_{E_{t=0}}(r)}{D(r)w̃(r, n)}R_{c_{2}}(r, n)\sqrt{\frac{4πr^{2}}{N_{0}(r)D(r)}} dr)|_{r=rmin} ,

Q_{2} = −∂_{r}(R_{c_{2}}(r, n) ∫ \frac{\overline{N}_{E_{t=T}} (r) − \overline{N}_{E_{t=0}}(r)}{D(r)w̃(r, n)}R_{c_{1}}(r, n)\sqrt{\frac{4πr^{2}}{N_{0}(r)D(r)}} dr

−R_{c_{1}}(r, n) ∫ \frac{\overline{N}_{E_{t=T}} (r) − \overline{N}_{E_{t=0}}(r)}{D(r)w̃(r, n)}R_{c_{2}}(r, n)\sqrt{\frac{4πr^{2}}{N_{0}(r)D(r)}} dr)|_{r=rmax} .

Then by Cramer’s Rule, the functions c_{1}(n) and c_{1}(n) have values

c_{1}(n) = \frac{P_{22}Q_{1} − P_{21}Q_{2}}{P_{11}P_{22} − P_{21}P_{12}} , c_{2}(n) = \frac{P_{11}Q_{2} − P_{12}Q_{1}}{P_{11}P_{22} − P_{21}P_{12}} .  (17)

Hence the cefficients in (17) are the functions of n necessary for the given time-transformed boundary conditions of (15) to be met in (14). Finally, substituting the coefficients in (17) into (14), and applying the inverse finite Fourier transform to N̂_{E} in (14) yields the operational WKB analytical solution to the boundary and initial/final value problem (2) in the (r, θ, ϕ, t)-domain.

\overline{N}_{E}(r, θ, ϕ, t) = F_{T}^{−1}[N̂_{E}(r, θ, ϕ, n)]

= \frac{1}{T} (N̂_{E}(r, θ, ϕ, 0) + ∑_{n=1}^{∞}[N̂_{E}(r, θ, ϕ, n)e^{inωt} + N̂_{E}(r, θ, ϕ, −n)e^{−inωt}]) ,

with N̂_{E}(r, θ, ϕ, n) given by (14).

References

[1] J. Beechem, E. Haas, Simultaneous determination of intramolecular distance distributions and conformational dynamics by global analysis of energy trnasfer measurements, Biophys J., vol. 55, 1989, pp. 1225-1236.

[2] W. Kaplan, Operational Methods for Linear Systems, Addison-Wesley 1962.

[3] I. S. Gradshteyn, I. M. Ryzhik, Table of Integrals, Series, and Products, Academic Press, p. 1013.

[4] P. Morse, W. Feshbach, Methods of Theoretical Physics, vol. 2, McGraw-Hill.

Valentino Anthony Simpao Mathematical Consultant Services 108 Hopkinsville St. Greenville, Kentucky 42345 USA e-mail: mcs007@muhlon.com Tel.: 502-338-5543