About this Journal Submit a Manuscript Table of Contents
Journal of Applied Mathematics
Volume 2013 (2013), Article ID 685137, 21 pages
http://dx.doi.org/10.1155/2013/685137
Research Article

Numerical Approximation of Higher-Order Solutions of the Quadratic Nonlinear Stochastic Oscillatory Equation Using WHEP Technique

1Engineering Mathematics & Physics Department, Faculty of Engineering, Cairo University, Giza, Egypt
2Mathematics Department, Faculty of Science, Northern Borders University, Arar, Saudi Arabia

Received 21 March 2013; Accepted 4 July 2013

Academic Editor: Turgut Öziş

Copyright © 2013 Mohamed A. El-Beltagy and Amnah S. Al-Johani. 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.

Abstract

This paper introduces higher-order solutions of the stochastic nonlinear differential equations with the Wiener-Hermite expansion and perturbation (WHEP) technique. The technique is used to study the quadratic nonlinear stochastic oscillatory equation with different orders, different number of corrections, and different strengths of the nonlinear term. The equivalent deterministic equations are derived up to third order and fourth correction. A model numerical integral solver is developed to solve the resulting set of equations. The numerical solver is tested and validated and then used in simulating the stochastic quadratic nonlinear oscillatory motion with different parameters. The solution ensemble average and variance are computed and compared in all cases. The current work extends the use of WHEP technique in solving stochastic nonlinear differential equations.

1. Introduction

Analysis of the response of linear and nonlinear systems subjected to random excitations is of considerable interest to the fields of mechanical and structural engineering [1]. Stochastic differential equations based on the white noise process provide a powerful tool for dynamically modeling complex and uncertain aspects. In many practical situations, it may be appropriate to assume that the nonlinear term affecting the phenomena under study is small enough; then its intensity is controlled by means of a small parameter, say [2].

According to [3], the solution of stochastic partial differential equations (SPDEs) using Wiener-Hermite expansion (WHE) has the advantage of converting the problem to a system of deterministic equations that can be solved efficiently using the standard deterministic numerical methods. The main statistics, such as the mean, covariance, and higher-order statistical moments, can be calculated by simple formulae involving only the deterministic Wiener-Hermite coefficients. In WHE approach, there is no randomness directly involved in the computations. One does not have to rely on pseudorandom number generators, and there is no need to solve the stochastic PDEs repeatedly for many realizations. Instead, the deterministic system is solved only once.

The application of the WHE [410] aims at finding a truncated series solution to the solution process of a stochastic differential equation. The truncated series is composed of two major parts: the first is the Gaussian part which consists of the first two terms, while the rest of the series constitutes the non-Gaussian part. In nonlinear cases, there exist always difficulties in solving the resultant set of deterministic integrodifferential equations got from the applications of a set of comprehensive averages on the stochastic integrodifferential equation obtained after the direct application of WHE. Many authors introduced different methods to face these obstacles. Among them, the WHE with perturbation (WHEP) technique [4] was introduced using the perturbation technique to solve perturbed nonlinear problems.

The WHE was originally started and developed by Wiener in 1938 and 1958 [11]. Wiener constructed orthonormal random bases for expanding homogeneous chaos depending on white noise and used it to study problems in statistical mechanics. Cameron and Martin [12] developed a more explicit and intuitive formulation for WHE (now known as Wiener chaos expansion, WCE). Their development is based on an explicit discretization of the white noise process through its Fourier expansion, which was missed in Wiener’s original formalism. This approach is much easier to understand and more convenient to use and hence replaced Wiener’s original formulation. Since Cameron and Martin’s work, WHE has become a useful tool in stochastic analysis involving white noise (Brownian motion) [3]. We will denote it by Imamura formulation [13]. Another formulation was suggested and applied by Imamura and his coworkers [13, 14]. They have developed a theory of turbulence involving a truncated WHE of the velocity field. The randomness is taken up by a white-noise function associated, in the original version of the theory, with the initial state of the flow. The mechanical problem then reduces to a set of coupled integrodifferential equations for deterministic kernels. In [1], the WHE (Imamura formulation [13]) was used to compute the nonstationary random vibration of a Duffing oscillator which has cubic nonlinearity under white-noise excitation. Solutions up to second order are obtained by solving the equivalent deterministic system by an iterative scheme. El-Tawil and his coworkers [410, 15, 16] used the WHEP technique to solve a perturbed nonlinear stochastic diffusion equation and the quadratic and cubic nonlinear stochastic oscillatory equation with first-order approximations.

In [17], the analysis of nonlinear random vibration has been studied using several methods, such as equivalent linearization method [18], stochastic averaging method [19], the WHE approach with nonstationary excitations [1], the WHEP technique [15], eigenfunction expansions [20], and the method of detailed balance [21]. All of the above methods are applied and used for nonlinear random oscillations of real systems subjected to random nonstationary (or stationary) excitations.

According to [5, 6], quadrate oscillation arises through many applied models in applied sciences and engineering when studying oscillatory systems [22]. These systems can be exposed to a lot of uncertainties through the external forces, the damping coefficient, the frequency, and/or the initial or boundary conditions. These input uncertainties cause the output solution process to be also uncertain. For most of the cases, getting the probability density function (p.d.f.) of the solution process may be impossible. So, developing approximate techniques (through which approximate statistical moments can be obtained) is an important and necessary work. There are many techniques which can be used to obtain statistical moments of such problems. The main goal of this paper is to introduce higher-order WHEP solutions and to suggest a numerical solver suitable to handle the equivalent deterministic system.

In [16], the WHEP technique is generalized to th nonlinearity, general order of WHE, and general number of corrections. Also, the extension to handle white noise in more than one variable and general nonlinearities are outlined. The generalized algorithm is implemented and linked to MathML [23] script language to print out the resulting equivalent deterministic system.

In the current work, the WHE formulation suggested by Meecham and his coworkers (Imamura formulation) is used to solve the stochastic nonlinear differential models of the form with the proper set of initial conditions which will be assumed to be deterministic. The differential operator is a general linear operator. The nonlinearity is introduced as losses of degree strengthened by a deterministic small parameter (). For the quadratic nonlinearity, will be equal 2. The uncertainty is introduced through white noise scaled by a deterministic envelope function . The function is a deterministic forcing function. Theorem 1 will be used in the derivation of the WHEP technique.

Theorem 1. The solution of (1), if exists, is a power series in ; that is, .

Proof. Using Picard’s method which generates a sequence of approximations that converge to the solution; that is, the solution is obtained as where the th approximation is computed as apply the inverse operator, , on both sides to get Let which is the solution at .
Then the Picard’s th approximation is computed as Now, we need to prove that the th approximation is a power series in ; that is, it can be written as where is the number of terms in the series. Using the mathematical induction, we need to prove that a power series solution will be obtained at and at provided that the th approximation is a power series solution.
At , the Picard’s 1st approximation will be , which is a power series in . Now we need to prove that is a power series in given that is a power series in . The Picard’s th approximation is computed as or, after substituting the power series of , we get The second term in the right hand side can be expanded using the multinomial theorem to get where and the counter runs over all the combinations of the positive integers such that . In a more simplified form, we can write where . Let ; then the expansion takes the form Substitute in Picard’s th approximation to get or which is a power series in . This completes the proof. We note that the previous theorem is applied also in case of deterministic force term; that is, . Also, the theorem is applied if the unknown function is a function of more than one variable [16].

This paper is organized as follows. In Section 2, the WHEP technique is reviewed and the generalized WHEP derivation steps are outlined. The deterministic set of equations equivalent to the stochastic nonlinear oscillatory equation are tabulated in Section 3. The analytical and numerical solutions of the oscillatory equation are described in Section 4. The simulations up to third order and fourth correction are shown in Section 5.

2. WHEP Technique

As a consequence of the completeness of the Wiener-Hermite set [1], any arbitrary stochastic process can be expanded in terms of the Weiner-Hermite polynomial set, and this expansion converges to the original stochastic process with probability one.

The solution function can be expanded in terms of Wiener-Hermite functionals as [4] or after eliminating the parameters, for the sake of brevity, we get where and is a -dimensional integral over the variables . The first term in the expansion (14) is the nonrandom part or ensemble mean of the function. The first two terms represent the normally distributed (Gaussian) part of the solution. Higher-order terms in the expansion depart more and more from the Gaussian form [16]. The Gaussian approximation is usually a bad approximation for nonlinear problems, especially when high-order statistics are concerned [18].

The components are called the deterministic kernels of the WHE of . They are functions of time and space variables and fully account for the time dependence of as well as for its statistical properties [3]. is a random output of a triple probability space , where is a sample space, is a -algebra associated with , and is a probability measure. For simplicity, will be dropped later on.

The functional is the th order Wiener-Hermite time-independent functional. The WH functionals form a complete set [1], and they satisfy the following recurrence relation for : with and : the white noise. By construction, the Wiener-Hermite functionals are symmetric in their arguments and are statistically orthonormal with respect to the weighting function ; that is, The average of almost all Wiener-Hermite functionals vanishes, particularly, The expectation and variance of the solution will be The WHE method can be elementarily used in solving stochastic differential equations by expanding the solution as well as the stochastic input processes via the WHE. The resultant equation is more complex than the original one due to being a stochastic integrodifferential equation. Taking a set of ensemble averages together with using the statistical properties of the WHE functionals, a set of deterministic integrodifferential equations are obtained in the deterministic kernels , . To obtain approximate solutions of these deterministic kernels, one can use perturbation theory in the case of having a perturbed system depending on a small parameter . Expanding the kernels as a power series of , another set of simpler iterative equations in the kernel series components are obtained. This is the main algorithm of the WHEP algorithm [4]. The technique was successfully applied to several nonlinear stochastic equations; see for example [4, 5, 10].

The WHEP technique for general nonlinear exponent (), general order (), and general number of corrections () follows the following steps [16].(1)Truncate the expansion (14) to contain only () kernels , ; that is, .(2)Substitute into the stochastic oscillatory equation (1).(3)Use the multinomial theorem to expand the quadratic nonlinear term in (1), , .(4)Multiply by , , and then apply the ensemble average. This will lead to equations in the kernels , .(5)For each kernel , , apply the perturbation technique up to corrections; that is, . (6)Equate the coefficients of , , in both sides to get equations for each kernel , .

This will lead to the following equations: where And the expectations are computed as It was explained in [16] how to get in terms of the Dirac delta functions and then use them to reduce the integrals that appear in . The summation means that all variations , , , that satisfy the equality are selected. This can be done be a searching technique. For these variations, the factors will be multiplied by each other to get ; that is, . The Kronecker delta functions are defined as The counter , in the summation in the right hand side of (20), runs over all of the combinations of the positive integers such that .

Equations (19) and (20) can always be solved using the proper sequence. The first equations of (19) are solved independently to get , ; then they are used to compute the other components in (20). For , the component is obtained by solving with the original initial and boundary conditions. For , the component is obtained by solving with zero initial and boundary conditions. The other components , , will be zeros due to zero right hand side and zero initial and boundary conditions. Equation (20) specifies the solution sequence to be followed. The component is evaluated in terms of the previously computed components , , . This means that the 1st corrections for all kernels, ,  , are solved firstly and then the 2nd corrections for all kernels, , up to the th corrections for all kernels , .

These results are consistent with the known results obtained using WHE. In WHE, higher-order kernels are driven by lower-order kernels, and at the bottom, the Gaussian kernels are driven by the random forcing directly. So, the lower-order kernels are usually dominant in magnitude [3].

The statistical properties of the solution will be calculated as If , then it will be convergent if [16] for . This means that should obey an upper bound condition after which divergence is obtained.

3. The Equivalent Deterministic Equations

Apply the previous WHEP algorithm to get the following systems of equations of the quadratic () nonlinear oscillatory equation and first-order () Gaussian approximation and different number of corrections (). The initial conditions are assumed deterministic and hence only the zero-order and zero-correction kernel equation will have the initial conditions ( and ). Other kernels equations will have zero initial conditions.

For the quadratic nonlinear oscillatory stochastic equation, the application of the WHEP technique will result in the following set of equations (Tables 1, 2, and 3). The equations will be written for the first, second, and third orders. For a certain correction level, the deterministic system will include also the equations from previous levels.

tab1
Table 1: The equivalent deterministic system for first-order approximation with different number of corrections (NC).
tab2
Table 2: The equivalent deterministic systems for second-order approximation with different number of corrections (NC).
tab3
Table 3: The equivalent deterministic system for third-order approximation with different number of corrections (NC).

In case of zero initial conditions, and ; that is, , and we will have the following reduced system of equations:

4. Oscillatory Equation and the Numerical Solver

For the linear oscillatory equation, the linear operator will be This means that the linear oscillatory equation can be written as . The parameter is the undamped angular frequency of the oscillator and is the damping ratio. The initial conditions will be considered as and . In case of zero initial conditions, the exact solution that can be obtained using different methods such as the theory of linear differential equations or the Laplace transform will be the convolution where with , assuming underdamping ().

For , the solution will be , which results in The numerical solution can be obtained for a model equation and then used for all kernels in the proper sequence. The model equation in this case will take the form where and are assumed constant. We can use any difference scheme, but as we are working with a white noise, the Dirac delta function is expected to appear in the equations of some kernels. So, an integral numerical scheme such as FEM or FVM will be more suitable as the integration of the Dirac delta function is easier to be handled. In our case, the finite volume method (FVM) will be considered. The time axis where will be discretized into equal intervals of size . The interval extending from to is taken as the control volume. This technique in one dimension will be equivalent to the trapezoidal integration rule. The second-order linear equation (31) will be decomposed into two simultaneous first-order equations. This can be done by substituting , and then we will have Integrate (32) along the control volume to get Approximate the integral with the trapezoidal rule which is of accuracy proportional to to get Also, integrate (32) along the control volume and substitute with from (34) to get where . If the Dirac delta function appears in the right hand side, a special treatment for is considered. In this case when and when . Equation (35) will be solved at each node to get and then (34) computed to get .

The numerical solution can be validated with the exact solution in case of (see (30)). With , , , , and , the comparison in Figure 1 shows that the numerical solution has sufficient accuracy with a relative error of 0.03% at . The overall convergence rate of the developed scheme was tested, and it has a convergence order of 1.5 as shown in Figure 2.

685137.fig.001
Figure 1: Comparison between the exact and the numerical solutions at .
685137.fig.002
Figure 2: The convergence order of the developed numerical solver.

5. Results

The following output is simulated using the previous developed numerical solver. The solution (34) of the model equation (31) is used to get all kernels with the proper right hand side. The mean response and the response variance are then calculated from the kernels using (24).

Figures 3 and 4 show the first-order response mean and variance for different values of . The simulations are done for the case of zero initial conditions, zero deterministic excitation, and unit envelope function multiplied by the white noise. The angular frequency is and the damping ratio is . For the response mean, the 3rd and 4th corrections are coincident. For the variance, the 2nd and 3rd corrections are coincident and also the 4th and 5th corrections.

fig3
Figure 3: First-order response mean and variance for (a) , (b) , and (c) .
fig4
Figure 4: First-order response mean and variance for (a) and (b) .

As it is shown in the figures, the nonlinearity strength greatly affects the amplitudes of the mean and variance. It should not be increased after a certain value to obtain a convergent solution. This value depends on the different parameters of the problem. As the nonlinearity strength increases, we need higher number of corrections. For and at seconds, the ratio between the response mean of each correction and the proceeding one is around 0.98, which is greater than . This means that the condition of convergence is satisfied and hence the solution converges in this case. Also, the convergence condition is satisfied for and . In case , the ratio between the 4th and 5th corrections is 0.6 and between the 2nd and 3rd corrections is 0.58. Both ratios are lesser than 0.7, which means that we will have a divergent solution for . Also, will produce a divergent solution.

Figures 5 and 6 show the first-order response mean and variance using the same parameters in Figures 3 and 4, but the envelope function is taken as . As it is clear from the figures, the effect of the nonlinearity strength on the response variance is negligible. The multiplication of attenuates the effect of the white noise and it becomes negligible with the time increase. The variance vanishes with the time and the solution becomes nearly deterministic. Also, the response variance is nearly invariant with the number of corrections.

fig5
Figure 5: The first-order response mean and variance for (a) , (b) , and (c) . Case of exponential function multiplied by the white noise.
fig6
Figure 6: The first-order response mean and variance for (a) and (b) . Case of exponential function multiplied by the white noise.

Figures 7 and 8 show the second-order response mean and variance using the same parameters in the first-order simulations, Figures 3 and 4. Figure 9 shows the second-order response mean and variance for longer time interval, seconds. This was done to ensure that seconds is a sufficient interval for the response mean and variance to reach their steady state values. For the second-order approximation, the response variance is computed practically as Figures 10 and 11 show the third-order response mean and variance with the same parameters used in the first- and second-orders simulations. For the third-order approximation, the response variance is computed practically, for time and memory saving, as Table 4 shows a comparison between the computing times for different levels of correction for the third-order approximation. The computing times for the response mean and variance are also shown. The table displays also the order of integrals in each level. The time step was 0.2 seconds and the computations are done on Intel Core i5, 2.4 GHz, machine with W7, 32 bits. We can note that around 90% of the computing time is consumed in level 4 (4th correction level). This is due to the higher orders of integrals computed in this level. As the level of correction increases, the order of integral increases and hence the computational time will also increase.

tab4
Table 4: The computing time for each level and for the response mean and variance compared with the total time.
fig7
Figure 7: The second-order response mean and variance for (a) , (b) , and (c) .
fig8
Figure 8: The second-order response mean and variance for (a) and (b) .
fig9
Figure 9: The second-order response mean and variance for and time interval of 20 seconds.
fig10
Figure 10: The third-order response mean and variance for (a) , (b) , and (c) .
fig11
Figure 11: The third-order response mean and variance for (a) and (b) .

Figures 12 and 13 show a comparison, at , between the response mean and variance for different orders and different number of corrections. As it was described earlier, the solution is convergent at . As it is shown in the figures, the response variance is more sensitive to the approximation order than the response mean. The variance converges as the approximation order increases. This means that, for a convergent solution, higher-order approximations are required for accurate prediction of the stochastic response.

fig12
Figure 12: Comparison between the response mean and variance for the first, second, and third orders with number of corrections 1 and 2, .
fig13
Figure 13: Comparison between the response mean and variance for the first, second, and third orders with numbers of corrections 3 and 4, .

Figures 14 and 15 show a comparison, at , between the response mean and variance for different orders and different number of corrections. In this case, the solution is divergent. The response variance diverges as the approximation order increases and also as the correction level increases.

fig14
Figure 14: Comparison between the response mean and variance for the first, second, and third orders with numbers of corrections 1 and 2, .
fig15
Figure 15: Comparison between the response mean and variance for the first, second, and third orders with numbers of corrections 3 and 4, .

It is worth to note that the WHEP technique used in the current work can be extended to solve stochastic PDEs with white noise in multiple dimensions and of different colors as described in [16].

6. Summary and Conclusions

The mean response of the quadratic nonlinear oscillatory system subjected to nonstationary random excitation was investigated using WHEP technique. The equivalent deterministic equations are derived up to third order. The solution is approximated up to fifth correction for the first order and up to fourth correction for the second and third orders. Numerical integral solution of the equivalent deterministic set of equations was applied using the FVM. The numerical treatment is validated after comparing the results with the analytical solution. The numerical solver is utilized in simulating the mean and variance of the nonlinear stochastic oscillatory motion with higher order, higher number of corrections, and different strengths of the nonlinear term. The values of the nonlinearity strength required for convergent solution are estimated. It was found that the numerical solution is efficient, and higher-order approximations are required for accurate prediction of the stochastic response.

References

  1. A. Jahedi and G. Ahmadi, “Application of Wiener-Hermite expansion to nonstationary random vibration of a Duffing oscillator,” Transactions of the ASME. Journal of Applied Mechanics, vol. 50, no. 2, pp. 436–442, 1983. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  2. J. C. Cortes, J. V. Romero, M. D. Rosello, and R. J. Villanueva, “Applying the Wiener-Hermite random technique to study the evolution of excess weight population in the region of Valencia (Spain),” American Journal of Computational Mathematics, vol. 2, no. 4, pp. 274–281, 2012. View at Publisher · View at Google Scholar
  3. W. Lue, Wiener chaos expansion and numerical solutions of stochastic partial differential equations [Ph.D. thesis], California Institute of Technology, Pasadena, Calif, USA, 2006.
  4. M. A. El-Tawil, “The application of WHEP technique on stochastic partial differential equations,” International Journal of Differential Equations and Applications, vol. 7, no. 3, pp. 325–337, 2003. View at Zentralblatt MATH · View at MathSciNet
  5. M. A. El-Tawil and A. S. Al-Johani, “Approximate solution of a mixed nonlinear stochastic oscillator,” Computers & Mathematics with Applications, vol. 58, no. 11-12, pp. 2236–2259, 2009. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  6. M. A. El-Tawil and A. S. El-Johani, “On solutions of stochastic oscillatory quadratic nonlinear equations using different techniques, a comparison study,” Journal of Physics, vol. 96, no. 1, 2008.
  7. M. A. El-Tawil and A. F. Fareed, “Solution of stochastic cubic and quintic nonlinear diffusion equation using WHEP, Pickard and HPM methods,” Open Journal of Discrete Mathematics, vol. 1, no. 1, pp. 6–21, 2011. View at Publisher · View at Google Scholar · View at MathSciNet
  8. M. A. El-Tawil and A. El-Shikhipy, “Approximations for some statistical moments of the solution process of stochastic Navier-Stokes equation using WHEP technique,” Applied Mathematics & Information Sciences, vol. 6, no. 3, pp. 1095–1100, 2012. View at MathSciNet
  9. M. A. El-Tawil and A. A. El-Shekhipy, “Statistical analysis of the stochastic solution processes of 1-D stochastic navier-stokes equation using WHEP technique,” Applied Mathematical Modelling, vol. 37, no. 8, pp. 5756–5773, 2013. View at Publisher · View at Google Scholar
  10. A. S. El-Johani, “Comparisons between WHEP and homotopy perturbation techniques in solving stochastic cubic oscillatory problems,” in Proceedings of the AIP Conference, vol. 1148, pp. 743–752, 2010. View at Publisher · View at Google Scholar
  11. N. Wiener, Nonlinear Problems in Random Theory, MIT Press, John Wiley, Cambridge, Mass, USA, 1958. View at MathSciNet
  12. R. H. Cameron and W. T. Martin, “The orthogonal development of non-linear functionals in series of Fourier-Hermite functionals,” Annals of Mathematics, vol. 48, pp. 385–392, 1947. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  13. T. Imamura, W. C. Meecham, and A. Siegel, “Symbolic calculus of the Wiener process and Wiener-Hermite functionals,” Journal of Mathematical Physics, vol. 6, pp. 695–706, 1965. View at Publisher · View at Google Scholar · View at MathSciNet
  14. W. C. Meecham and D. T. Jeng, “Use of the Wiener-Hermite expansion for nearly normal turbulence,” Journal of Fluid Mechanics, vol. 32, pp. 225–235, 1968.
  15. E. F. Abdel-Gawad and M. A. El-Tawil, “General stochastic oscillatory systems,” Applied Mathematical Modelling, vol. 17, no. 6, pp. 329–335, 1993.
  16. M. A. El-Beltagy and M. A. El-Tawil, “Toward a solution of a class of non-linear stochastic perturbed PDEs using automated WHEP algorithm,” Applied Mathematical Modelling, vol. 37, no. 12-13, pp. 7174–7192, 2013.
  17. X. Yong, X. Wei, and G. Mahmoud, “On a complex duffing system with random excitation,” Chaos, Solitons and Fractals, vol. 35, no. 1, pp. 126–132, 2008. View at Publisher · View at Google Scholar
  18. P. Spanos, “Stochastic linearization in structural dynamics,” Applied Mechanics Reviews, vol. 34, pp. 1–8, 1980.
  19. W. Q. Zhu, “Recent developments and applications of the stochastic averaging method in random vibration,” Applied Mechanics Reviews, vol. 49, no. 10, pp. 72–80, 1996. View at Publisher · View at Google Scholar
  20. J. Atkinson, “Eigenfunction expansions for randomly excited nonlinear systems,” Journal of Sound and Vibration, vol. 30, no. 2, pp. 153–172, 1973. View at Publisher · View at Google Scholar
  21. A. Bezen and F. Klebaner, “Stationary solutions and stability of second order random differential equations,” Physica A, vol. 233, no. 3-4, pp. 809–823, 1996. View at Publisher · View at Google Scholar
  22. A. H. Nayfeh, Problems in Perturbation, John Wiley & Sons, New York, NY, USA, 1993. View at MathSciNet
  23. MathML, http://www.w3.org/Math/.