27. Stochastic Differential Equations¶
PDF pages 1167–1202
Chapter 27 Stochastic Differential Equations Now we will consider the numerical solution to stochastic differential equations (SDEs) of the form
(27.1) where dW(t) is the Wiener process as usual (see Chapter 17). We will stick to considering only It¯o SDEs for simplicity. Because of the unusual and singular nature of dW(t), the methods that apply to ODEs (Chapter 24) do not work well here, and we will have to develop new methods. We will also see that accurately solving SDEs is much more difficult than ODEs. However, we will then have the formalism to put stochastic Schrödinger and master equations (Chapters 18 and 19) on the computer.1 27.1 Stochastic Euler Method The simplest numerical method we considered in solving ODEs was the Euler method (24.2). This applies
(27.2) where again the solution is evolved in finite time steps of ∆t, and we are using the condensed notation
linear approximation to the stochastic term, to arrive at the stochastic Euler method (often called the Euler–Maruyama method):2
(27.3) (stochastic Euler method) We are defining the time increment as before, and the Wiener increment ∆Wn is defined in an analogous way: ∆tn := Z tn+1 tn dt′ ∆Wn := Z tn+1 tn dW(t′). (27.4) (time increments) Of course, ∆tn ≡∆t is independent of n, so we need not indicate explicit time dependence. The key, however, is that the approximation works when \alpha and \beta vary slowly over the interval [tn, tn+1], and thus 1One of the best and most complete references on this subject is Peter E. Kloeden and Eckhard Platen, Numerical Solution of Stochastic Differential Equations, 3rd ed. (Springer, 2000). 2Gisirō Maruyama, ‘‘Continuous Markov Processes and Stochastic Equations,’’ Rendiconti del Circolo Matematico di Palermo 4, 48 (1955).
27.1.1 Truncation Error¶
Chapter 27. Stochastic Differential Equations when we try to compute the solution
Z tn+1 tn
Z tn+1 tn dW(t′) \beta(y(t′), t′), (27.5) we can, to first approximation, treat \alpha and \beta as constants and pull them out of their respective integrals, which precisely yields Eq. (27.3). 27.1.1 Truncation Error Recall the It¯o chain rule (17.193) for a function f(y), where y(t) satisfies the SDE (27.1):
2f ′′(y) \beta2(y, t) dt + f ′(y) \beta(y, t) dW. (27.6) Integrating this expression from t′ = t0 to t,
Z t t0
2f ′′(y) \beta2(y, t′) dt′ + Z t t0 f ′(y) \beta(y, t′) dW(t′). (27.7) Now we can let f(y) −\rightarrow \alpha(y, t), f(y) −\rightarrow \beta(y, t), and t0 −\rightarrow tn, and then put the two resulting expressions in Eq. (27.5) to obtain
- Z tn+1 tn dt′ Z t′ tn dt′′
+ Z tn+1 tn dt′ Z t′ tn dW(t′′) \alpha′(y, t′′) \beta(y, t′′) + Z tn+1 tn dW(t′) Z t′ tn dt′′
+ Z tn+1 tn dW(t′) Z t′ tn dW(t′′) \beta′(y, t′′) \beta(y, t′′). (27.8) The primes here are equivalent to the partial derivatives \partial y. Again, we can pretend that the \alpha and \beta functions are constant over the short time interval ∆t. Then we have simple integrals of the form Z tn+1 tn dt′ Z t′ tn dt′′; Z tn+1 tn dt′ Z t′ tn dW(t′′); Z tn+1 tn dW(t′) Z t′ tn dt′′; Z tn+1 tn dW(t′) Z t′ tn dW(t′′), (27.9) to deal with. The first integral is just ∆t2/2, or for our purposes, simply O(∆t2). Similarly, if we proceed with counting each ∆W on average as equivalent to ∆t1/2, the second and third integrals are O(∆t3/2), and the last integral is O(∆t1). Clearly, any error in treating \alpha and \beta as constants will result in higher-order errors, so the truncation error here is O(∆t). Thus, we may again write the stochastic Euler method as
(27.10) (stochastic Euler method) The local truncation error here is, as we see, much worse than the deterministic Euler method (27.2). To estimate the global truncation error, again suppose we use this method to evolve the solution out to a fixed
the stochastic term. The dominant errors will thus be random and uncorrelated, so we consider the error in a random walk of N steps of mean size of order t/N, which will scale as \sqrt
\sqrt
Thus the global error of the stochastic Euler scheme converges very badly, as O(∆t1/2). This is indicative of the fact that stochastic differential equations are much more difficult to solve than the deterministic counterparts. In any case, as in the ODE case, we call this method an O(∆t1/2) method, because the local truncation is correct to this order, or equivalently, the global error is O(∆t1/2).
27.2.1 Multiplicative vs. Additive Noise¶
27.2 Milstein Method 27.2 Milstein Method To construct the next, relatively simple, higher-order method, we note that the only term we neglected in Eq. (27.8) that was O(∆t) was the last one, involving the double integral Z tn+1 tn dW(t′) Z t′ tn
Z tn+1 tn dW(t′) [W(t′) −W(tn)] = 1 Z tn+1 tn d[W(t′)]2 −dt
−W(tn) Z tn+1 tn dW(t′) = 1 W 2(tn+1) −W 2(tn) −∆t 2 −W(tn)∆W(tn) = 1 2 [W(tn+1) + W(tn)] ∆W(tn) −∆t 2 −W(tn)∆W(tn) = 1 2 [W(tn+1) −W(tn)] ∆W(tn) −∆t = 1 [∆W(tn)]2 −∆t
, (27.11)
Milstein method3 is to keep this correction term, keeping in mind that the factor we just worked out is multiplied by \beta(y, tn) \beta′(y, tn), which we are treating as constant to this order of approximation:
2\beta(yn, tn) \beta′(yn, tn) ∆W 2 n −∆t
(Milstein method) (27.12) Again, the prime refers to partial differentiation with respect to yn. Clearly the local truncation error is O(∆t3/2), since we have ignored the two integrals in Eqs. (27.9) of the same order. By the same argument as for the Euler method, the global error is a factor of ∆t1/2 worse, or O(∆t), which is the same global error as the deterministic Euler method. 27.2.1 Multiplicative vs. Additive Noise Note that the correction term in the Milstein method (27.12) is of the form \beta\beta′, and thus is only required for multiplicative noise (where \beta is a function of y). For additive noise (\beta independent of y), the correction term vanishes, and the stochastic Euler method has O(∆t) global error. In general, we expect additive noise to be easier to solve numerically than multiplicative noise, and one strategy to improve the accuracy of numerical methods is to use the Lamperti transform (Section 17.7.4.2) to change a multiplicative-noise process to an additive process, assuming that the transform can be inverted to give the original solution. 27.3 Stochastic Taylor Expansion What we have written down in Eq. (27.8) is something like a Tayler expansion for yn+1 in terms of yn. Again, treating \alpha and \beta as constants, essentially what we have is a Taylor expansion, neglecting O(∆t3/2) terms. The full Taylor expansion is given by iterating the above procedure, for example using Eq. (27.7) to replace the functions \alpha and \beta evaluated at time t′′. This procedure obviously introduces triple stochastic integrals; the next iteration introduces quadruple integrals, and so forth. 3G. N. Mil’shtein, ‘‘Approximate Integration of Stochastic Differential Equations,’’ Theory of Probability and its Applications 19, 557 (1974).
27.3.1 Single and Double Integrals¶
Chapter 27. Stochastic Differential Equations 27.3.1 Single and Double Integrals The first thing that is clear is that we will need to employ a hierarchy of stochastic integrals. The ones we have considered so far are the single integrals
Z tn+1 tn dt′ = ∆t
Z tn+1 tn
(27.13) as well as the double integrals
Z tn+1 tn dt′ Z t′ tn dt′′ = ∆t2
Z tn+1 tn dt′ Z t′ tn dW(t′′)
Z tn+1 tn dW(t′) Z t′ tn
Z tn+1 tn dW(t′) Z t′ tn
n −∆t 2 . (27.14) The first and the last integral we worked out before, and the third integral follows from
Z tn+1 tn dW(t′) (t′ −tn) = Z tn+1 tn dW(t′) t′ −tn∆Wn
tn − Z tn+1 tn dt′ W(t′) −tn∆Wn
Z tn+1 tn dt′ [W(t′) −W(tn)] −∆t W(tn) −tn∆Wn
(27.15) where we integrated by parts in the third step. Note that all the double integrals can be reduced to expressions in terms of single integrals and (I10)n; however, this last double integral is irreducible in the sense that it cannot be written only in terms of single integrals. We can thus characterize it more completely. It is clearly Gaussian, as from its definition it is a sum over independent, Gaussian random variables. It has mean, variance, and covariance with ∆Wn given by
(I10) 2 n
= ∆t3
2 . (27.16)
27.3.2 Multiple It¯o Integrals¶
27.3 Stochastic Taylor Expansion The mean is obvious, as (I10)n is again the sum over independent, Gaussian random variables of zero mean. The variance and covariance we leave as an exercise (Problem 27.1). For the purposes of simulation, given two independent, standard-normal random numbers z1 and z2 (i.e., variance 1 and mean 0), it is not hard to verify that the linear combinations ∆Wn = \sqrt ∆t z1
z1 + z2 \sqrt (27.17) have the correct statistics for the two desired quantities. However, we will return to more useful strategies for computing these numbers below. 27.3.2 Multiple It¯o Integrals The notation here generalizes readily to further multiple integrals. In the case of I0, the zero subscript indicates a simple integration of dt over the time interval; for I1, the unit subscript indicates a simple integration of dW(t) instead. For the double-integrals, the subscripts have the same meaning, but the result of integrating according to the first index is integrated again according to the second integral; thus, I10 means to integrate dW, and then integrate the result with a the differential dt. In the general case, with Ij1j2j3...jn, where j\alpha \in {0, 1}, we again integrate dt or dW according to the value of j1, then integrate the result over dt or dW according to the value of j2, and so on. Clearly, when we iterate the procedure leading to Eq. (27.8) to generate the stochastic Taylor expansion, we are generating on the very next iteration triple integrals, such as I110. Roughly speaking, when counting order we should count each 1 as a 1/2 order, while a 0 is a full extra order. Thus, the Taylor expansion to order ∆t involves I0, I1, and I11. To order ∆t3/2, the expansion also involves I10, I01, and I111; while to order ∆t2, the expansion also involves I00, I110, I101, I011, and I1111. It is not hard to see that the stochastic Taylor expansion is much more complicated than the regular Taylor expansion: the regular Taylor expansion only involves I0, I00, I000, and so on. In terms of actual calculations, there is the additional complication that at higher orders, new random quantities are introduced that cannot be expressed entirely in terms of lower-order quantities. Again, these integrals may be written in terms of other integrals, although this becomes complicated where many indices are involved; however, it can be shown that4
n X i=0 Ij1...ji−1jji+1...jn + n X i=1 j \deltajji∆W Ij1...ji−10ji+1...jn, (27.18) where recall that j and all ji are either 0 or 1. It is precisely the existence of extra terms that causes traditional numerical methods for ODEs to fail in general for SDEs. In certain special cases (such as additive noise), ordinary methods may provide better performance (since certain terms in the stochastic Taylor expansion will vanish). However, in general different methods must be developed to handle SDEs at ‘‘high’’ order.5 4Kloeden and Platen, op. cit., Proposition 5.2.3, p. 170. 5see Kevin Burrage, Pamela Burrage, Desmond J. Higham, Peter E. Kloeden, and Eckhard Platen, ‘‘Comment on ’Numerical methods for stochastic differential equations,’ ’’ Physical Review E 74, 068701 (2006) (doi: 10.1103/PhysRevE.74.068701), which is a comment on the use of standard Runge–Kutta methods for SDEs by Joshua Wilkie, ‘‘Numerical methods for stochastic differential equations,’’ Physical Review E 70, 017701 (2004) (doi: 10.1103/PhysRevE.70.017701).
27.3.4 Multiple Wiener Processes¶
Chapter 27. Stochastic Differential Equations 27.3.3 Expression for the Taylor Expansion As we have indicated thus far, the stochastic Taylor expansion is quite complicated. We will simply state the result here for the It¯o–Taylor expansion6 for the solution y(t) to the scalar SDE (27.1), which is
\infty X n=1 X j1...jn=0 Ij1...jnfj1...jn, (27.19) (It¯o–Taylor expansion) where the multiple stochastic integrals Ij1...jn are defined over the time interval from t0 to t, and the coefficient functions are defined recursively in terms of the functions with one index removed: fj1...jn = Lj1fj2...jn, (27.20)
by
2\beta2\partial 2 y
(27.21) Thus, for example,
2\beta2\partial 2 y \beta. (27.22) Clearly, these coefficient functions become much more complicated as the order increases. However, it should also be reasonably clear how this expansion comes about by the iterated application of the It¯o integral forumula. This is one case where Stratonovich calculus is fairly nice: the corresponding Stratonovich formulae for the multiple integrals and coefficients are simplified, for example, by the absence of the final term in the expression for L0. 27.3.4 Multiple Wiener Processes As you can imagine, a vector SDE driven by multiple Wiener processes dWk is way more complicated than the scalar case. In particular, the multiple integrals must be generalized to indices beyond 0 and 1, so that there are many more terms in the Taylor expansion, and the differential operators must be similarly generalized. In any case, some extensions of the methods we show below to vector Wiener processes are given by Kloeden and Platen7 27.4 Stochastic Runge–Kutta Methods The Milstein method (27.12) is the simplest example of ‘‘Taylor’’ methods that explicitly cancel the higher- order terms in the It¯o–Tayler expansion. The problem with these methods is that they require not only specifying the functions \alpha and \beta in the SDE (27.1), but also their derivatives. In the Milstein case, we have to specify \beta′, but for higher-order methods, more and more derivatives must be specified. This is not a problem in principle, but merely for convenience of the user. Instead, we can try an approach following the ODE case of Runge–Kutta methods, where multiple function evaluations with different arguments are used to estimate the information that would otherwise 6Kloeden and Platen, op. cit., Section 5.5, p. 181. The It¯o–Taylor expansion given in this reference is much more rigorous than what we have written here; here, we have written what amounts to an infinite number of iterations of the It¯o integral formula. Truncations of the expansion should in principe be done carefully. 7Kloeden and Platen, op. cit.
27.4 Stochastic Runge–Kutta Methods be supplied by the derivatives. Again, though, these are not the same Runge–Kutta methods as in the deterministic case, but different methods designed to cancel error terms in the It¯o–Taylor expansion above. We will simply quote two It¯o methods here.8 An O(∆t) method (i.e., a stochastic Runge–Kutta method of comparable accuracy to the Milstein method) is
\sqrt ∆t
\sqrt ∆t
n −∆t]. (stochastic Runge–Kutta method, order 1) (27.23) An O(∆t3/2) method is
\sqrt ∆t
\sqrt ∆t
\sqrt ∆t
- 1
- \sqrt ∆t
∆W 2 n −∆t +
+
∆W 2 n −∆t ∆Wn. (stochastic Runge–Kutta method, order 3/2) (27.24) At the time of writing, there seems to have been no true O(∆t2) methods that have been reported for the general (multiplicative noise) case. To see the relative performance of the schemes we have presented thus far, we test them on the sample problem9
(27.25) where a and b are constants. This is a multiplicative SDE with the analytic solution
(1 + y0) exp[−2at + 2bW(t)] + 1 −y0). (27.26) For this calculation we took a = b = 1 and y0 = 0. The base step size was ∆t = 0.00015, with the other runs having the same base step size reduced by powers of 2, down to ∆t = 0.00015/215. The solutions were
trajectory in each case. Every calculation was performed on the same Brownian path W(t) (see below). The following plot compares the accuracy of the stochastic Euler method (27.3), the O(∆t) Runge–Kutta method (27.23), and the O(∆t3/2) Runge–Kutta method (27.24). 8Kloeden and Platen, op. cit.; P. E. Kloeden and E. Platen, ‘‘Higher-Order Implicit Strong Numerical Schemes for Stochastic Differential Equations,’’ Journal of Statistical Physics 66, 283 (1992) (doi: 10.1007/BF01060070). 9K. Burrage and P. M. Burrage, ‘‘High strong order explicit Runge–Kutta methods for stochastic ordinary differential equations,’’ Applied Numerical Mathematics 22, 81 (1996); Peter E. Kloeden and Eckhard Platen, Numerical Solution of Stochastic Differential Equations, 3rd ed. (Springer, 2000).
Chapter 27. Stochastic Differential Equations -8.6 -3.6 Dt -8 -7 -6 -5 -4 -3 -13 rms error -12 -11 -10 -9 -8 -7 -6 -5 -4 order 0.5 order 1.0 order 1.5 The scaling behavior in each case is clear, though again the convergence is slow compared to ODEs, and the step sizes required to achieve high accuracy are quite small. Note that in the O(∆t3/2) case, rounding dominates the error at the smallest step sizes. In terms of cpu time, the higher order methods are again superior in this problem. -2.8 cpu time -1 -2 -2.5 -13 rms error -12 -11 -10 -9 -8 -7 -6 -5 -4 -3 order 0.5 order 1.0 order 1.5 In this case, each method takes a similar amount of time for a given time step ∆t. This is because the test SDE is very simple, and the cpu effort is dominated by random-number generation and construction of the Brownian path. For very complicated SDEs, there will obviously be more of a difference in run times between the different methods. However, it is clear that higher-order methods still provide substantial benefit despite the added complexity. 27.5 Implicit Schemes As we discussed in Section 24.2.2, implicit methods have some advantages in certain classes of deterministic ODEs. It is therefore natural to try to construct similar implicit methods for SDEs. Unfortunately, we can quickly run into problems. Consider the simple SDE dy = −ay dt −by dW (27.27)
27.5 Implicit Schemes with constant coefficients a and b. The stochastic Euler method for this SDE reads yn+1 = yn −ayn ∆t −byn ∆Wn. (27.28) The obvious generalization of the impicit Euler method for deterministic ODEs is to evaluate the functions on the right-hand side at time tn+1 instead of t: yn+1 = yn −ayn+1 ∆t −byn+1 ∆Wn. (27.29) Note that ∆Wn is unchanged, because it must be chosen according to the convention for It¯o SDEs. Solving this equation for yn+1,
yn
. (27.30) Clearly there is a problem here: ∆Wn is a stochastic variable that can take on essentially any value with nonzero probability, and typically takes on values of the order of \sqrt ∆t. For small ∆t it may thus happen that the denominator can come very close to zero, when ∆Wn takes on values near −(a/b)∆t. Thus we lose the stability properties we gained in the deterministic case, due just to the stochastic nature of the Wiener process. The problem is apparently just with the stochastic part of the equation. There is no problem, for example, if we choose a hybrid scheme where the deterministic step is taken implicitly, but the stochastic parts are treated explicitly, as in yn+1 = yn −ayn+1 ∆t −byn ∆Wn (27.31) for the example problem, which becomes
yn 1 + a ∆t −byn ∆Wn. (27.32) which has no special problems with small denominators. Again, we will simply quote two implicit It¯o methods here,10 corresponding to the explicit O(∆t) and O(∆t3/2) methods above. An O(∆t) method (i.e., a stochastic Runge–Kutta method of comparable accuracy to the Milstein method) is
\sqrt ∆t
- \sqrt ∆t
n −∆t]. (implicit stochastic Runge–Kutta method, order 1) (27.33) Here, we have taken a simple average of \alpha(yn, tn) and \alpha(yn+1, tn+1). Of course, any weighted average could 10Kloeden and Platen, op. cit.; P. E. Kloeden and E. Platen, ‘‘Higher-Order Implicit Strong Numerical Schemes for Stochastic Differential Equations,’’ Journal of Statistical Physics 66, 283 (1992) (doi: 10.1007/BF01060070).
Chapter 27. Stochastic Differential Equations be taken, interpolating between explicit and fully implicit. Similarly, an O(∆t3/2) method is
\sqrt ∆t
\sqrt ∆t
- \sqrt ∆t
∆W 2 n −∆t +
- \sqrt ∆t
(I10)n −∆Wn ∆t +
∆W 2 n −∆t ∆Wn. (stochastic Runge–Kutta method, order 3/2) (27.34) This form is specific to the choice of a half ‘‘degree of implicitness’’ (that is, an average of \alpha at the present and advanced times). When comparing these methods to the corresponding explicit methods in the last section on the example problem (27.25), the performance is about the same (with a slight advantage in the order 1 case) for a given time step, and the cpu time is again about the same since the calculation is dominated by the random-number generation, not by the implementation of the finite-difference formulae or by the functional iteration in the implicit schemes. Of course, similar tricks can be done with any of the schemes we have presented so far: the deterministic step can be taken with high-order deterministic methods (such as fourth-order Runge–Kutta), so long as the stochastic parts are treated according to SDE-adapted methods as presented here. The formulae here and in the last section have been implemented in a publicly-available Fortran 90 module.11 This module includes the facilities described below to generate consistent Brownian paths with different step sizes, and a sample code implementing the test problem above is included. 27.6 Strong and Weak Convergence Up until now, we have been considering a particular kind of convergence of solutions of SDEs, strong convergence. This means we are considering pathwise convergence of solutions. If ˜y(t; ∆t) is a finite- difference approximation to the true solution y(t) to the SDE (27.1), then the method for generating ˜y(t; ∆t) is of strong order \gamma if at fixed t,
(27.35) (strong convergence condition) The order here again refers to the scaling behavior of the global error. Another less-demanding convergence criterion refers to convergence of ensemble means. The same numerical method is said to have weak order \delta if for every polynomial g(y),
(27.36) (weak convergence condition) at fixed time t. Strong convergence at order \gamma implies weak convergence of at least the same order. For example, we argued that the stochastic Euler method has strong order 1/2. This is due to an error term proportional to I11, which has zero mean; when considering expectation values, this error term is wiped out, and the stochastic Euler scheme is actually of weak order 1.12 The Milstein method turns out to have 11SDERK90, available online at http://atomoptics.uoregon.edu/~dsteck/computer.html. 12G. N. Milshtein, ‘‘A Method of Second-Order Accuracy Integration of Stochastic Differential Equations,’’ Theory of Prob- ability and its Applications 23, 396 (1978) (doi: 10.1137/1123045).
27.7 Consistent Brownian Paths both strong and weak order 1.13 Since we are normally thinking about solutions to SDEs as individually interesting objects (modeling individual realizations of a continuous measurement process, for example), we will generally be concerned with the more difficult case of strong convergence. 27.7 Consistent Brownian Paths In considering the strong convergence of the solution to an SDE, it is in general easy to check the convergence where the exact solution is known, as in the test problem above. However, this is obviously more difficult in the generic case where the exact solution is not known. In deciding whether or not to accept a particular numerical solution, in the deterministic case you would usually just run the integration again, but with a different step size (say half the original step size). If the two solutions match to within some tolerance, then you can accept the solution. But this is trickier with SDEs. Suppose that you use a random number generator to generate a sequence of Wiener increments ∆W0, ∆W1, ∆W2, . . . , (27.37) where
∆W 2 n
= ∆t. To cut the step size in half and rerun things, first of all, the same set of random numbers must be used, otherwise the two runs will not in general be well correlated. But even if we use the same set of random numbers to generate the new Wiener increments, ∆W (1/2) , ∆W (1/2) , ∆W (1/2) , ∆W (1/2) , ∆W (1/2) , ∆W (1/2) , . . . , (27.38)
n
random number in the two cases generated ∆W0 and ∆W (1/2) , while the second random number generated ∆W1 and ∆W (1/2) . However, ∆W1 and ∆W (1/2) don’t correspond to the same absolute time. In fact, for consistency, what we require is
;
;
, (27.39) and so on. These conditions allow both sequences to correspond to time integrals (I1)n of the same ideal, particular realization of the Wiener process W(t), which we will refer to as a particular Brownian path. Thus, what we require is a procedure for constructing sequences of Wiener increments with different step sizes, but corresponding to different paths. One convenient method for doing this is to start with the sequence of Wiener increments ∆Wn on the coarser time step ∆t, and then to refine it consistently onto the finer time grid of step ∆t/2.14 That is, given a Gaussian Wiener increment ∆Wn with variance ∆t, we will use another Gaussian random variable z of unit variance (i.e., standard normal) to generate two new, uncorrelated Gaussian random numbers ∆W (1) n and ∆W (2) n , each of variance ∆t/2 and satisfying ∆W (1) n + ∆W (2) n = ∆Wn. These conditions are clearly satisfied if we choose15 ∆W (1) n = 1 2∆Wn + \sqrt ∆t z ∆W (2) n = 1 2∆Wn − \sqrt ∆t z. (27.40) (refinement of Brownian path) This procedure may obviously be iterated to obtain consistent Brownian paths with yet smaller time steps (by powers of two), creating a tree structure of Wiener increments, with different levels corresponding to different step sizes. 13G. N. Milshtein, op. cit. 14J. G. Gaines and T. J. Lyons, ‘‘Variable Step Size Control in the Numerical Solution of Stochastic Differential Equations,’’ SIAM Journal on Applied Mathematics 57, 1455 (1997) (doi: 10.1137/S0036139995286515). 15Paul Lévy, Processus Stochastiques et Mouvement Brownien (Gauthier–Villars, 1948).
27.7.1 Consistent Multiple It¯o Integrals¶
Chapter 27. Stochastic Differential Equations ∆W0 ∆W 1/2 ∆W 1/4 ∆W 1/4 ∆W 1/2 ∆W 1/4 ∆W 1/4 ∆W1 ∆W 1/2 ∆W 1/4 ∆W 1/4 ∆W 1/2 ∆W 1/4 ∆W 1/4 Descending one level on this Brownian tree involves generating a new random number for each Wiener increment on the coarser level, and then combining them as in Eqs. (27.40) to generate the Wiener increments on the finer level. This procedure may be continued indefinitely to consistently realize Brownian paths with arbitrarily fine steps. This procedure allows you, for example, to run one simulation with one time step ∆t and another ∆t/2 on the same Brownian path, and then compare the two runs to estimate the numerical error. Or you can run many different step sizes and examine the convergence behavior of the solution. Consistency of the Brownian path is also an essential first step in implementing adaptive time-step methods for SDEs, which are much more difficult than for ODEs. (The stochastic Euler method, for example, as a basis for an adaptive method can lead to convergence to a wrong solution!16) 27.7.1 Consistent Multiple It¯o Integrals One remaining point to address is that for the order 3/2 Runge–Kutta methods, we also need to refine (I10)n consistently onto finer time grids. One conceptually simple way to generate the (I10)n in the first place is to generate the Wiener increments on a finer scale than needed for the integration. That is, if we intend to generate a finite-difference solution to an SDE with time step ∆t, Wiener increments ∆Wn, and multiple integrals (I10)n, then we should start by generating finer Wiener increments \deltaWn for step size \deltat < ∆t.
compute ∆Wn by the simple sum, ∆Wn = N\delta X j=1 \deltaWj, (27.41) where N\delta is the number of temporal substeps per integration step (20 or 2, as we mentioned above). This expression is exact, and hasn’t gotten us anything because we could have just picked ∆Wn directly. The point is that we may also compute (I10)n via the deterministic Euler method
Z tn+1 tn
N\delta X j=2 j X k=1 \deltaWk ! . (27.42)
However, (I10)n itself is O(∆t3/2), and so accounting for the coefficient of the global error term, the error
the order 1.5 methods (27.24) and (27.34). Furthermore, from the last section we know how to refine the increments \deltaWn onto a finer time grid, and thus we can consistently generate the (I10)n also with finer time steps on the same Brownian path. 16Gaines and Lyons, op. cit.
27.7 Consistent Brownian Paths 27.7.1.1 Lévy Areas A better approximation to the integral (I10)n is to compute the Lévy area17
Z tn+1 tn dW (j)(t′) Z t′ t dW (i)(t′′) − Z tn+1 tn dW (i)(t′) Z t′ t dW (j)(t′′) ! = 1 2 [(Iij)n −(Iji)n] , (Lévy area) (27.43) where in the notation here, dW 0 ≡dt and dW 1 ≡dW. We can also approximate this by the (stochastic) Euler scheme (which is equivalent to the Milstein scheme for these integrals),
N\delta X p=2 \deltaW (j) p p−1 X q=1 \deltaW (i) q − N\delta X p=2 \deltaW (i) p p−1 X q=1 \deltaW (j) q ! . (27.44) In particular, the Lévy area that we want to approximate is
N\delta X p=2 \deltat p−1 X q=1 \deltaWq − N\delta X p=2 \deltaWp p−1 X q=1 \deltat !
N\delta X p=2 p−1 X q=1 \deltaWq − N\delta X p=2 \deltaWp(p −1) ! . (27.45) Then since
2∆t ∆Wn, (27.46) where we used Eqs. (27.14), we can compute (I10)n based on this approximation to the Lévy area. To see why the Lévy area is better despite the extra complication of the formula, we can compute the variance of the numerical approximant (27.42) to (I10)n: DD (ˆI10) 2 n EE
(N\delta −1)2
\deltaW 2
\deltaW 2
\deltaW 2
\deltaW 2 N\delta = \deltat3
N\delta−1 X j=1 j2 = \deltat3 N 3 \delta + N 2 \delta −11N\delta + 1 = ∆t3 " 1 + 3 \deltat ∆t −11 \deltat ∆t 2 + 3 \deltat ∆t 3# . (27.47)
approximation gets variance of the generated (I10)n’s wrong by an O(\deltat/∆t) bias. 17Gaines and Lyons, op. cit.
Chapter 27. Stochastic Differential Equations For the Lévy area, we can also compute the variance of the numerical approximant: DD ( ˆA10) 2 n EE = \deltat2 (N\delta −1)2
\deltaW 2
\deltaW 2
\deltaW 2
\deltaW 2 N\delta = \deltat3
N\delta−1 X j=1 [j −(N\delta −j)]2 = \deltat3 N 3 \delta −4N\delta + 1 = ∆t3 " 1 −4 \deltat ∆t 2 + 3 \deltat ∆t 3# . (27.48) The true variance of the Lévy area is
(A10) 2 n
=
(I10) 2 n
- ∆t2
∆W 2 n
−2∆t
3 + ∆t2 4 ∆t−2∆t ∆t2 = 1 12∆t3, (27.49) and thus the approximant biases the variance by an error of O[(\deltat/∆t)2]. Since we can compute the integral (I10)n directly from the Lévy area by adding a quantity that can be computed exactly, the bias to the variance of (I10)n by this method is only O[(\deltat/∆t)2], which is one order better than the previous method. 27.7.1.2 Direct Refinement Another method for directly refining ∆Wn and (I10)n onto a finer temporal grid, in the spirit of Eqs. (27.40), is to combine them directly with two independent, standard-normal random numbers z1 and z2, according to the linear transformation18 ∆W (1) n (I10)(1) n ∆W (2) n (I10)(2) n = − \sqrt ∆t −1 2∆t ∆t3/2 \sqrt −∆t \sqrt ∆t −3 2∆t −∆t3/2 \sqrt ∆t3/2 ∆t −1 z1 z2 ∆Wn (I10)n . (refinement of Brownian path) (27.50) Here, the ∆W (1) n and ∆W (2) n are the Wiener increments for the two subintervals of duration ∆t/2, and (I10)(1) n and (I10)(2) n are the corresponding double integrals. This method obviously has the advantage of fewer extra random numbers generated and fewer arithmetic operations for each refinement, compared to the Lévy-area method above. This method also has no bias in terms of the variance of the refined variables. To see where this method comes from, consider the inverse problem, which is that if you have a Wiener path defined on a fine time grid, how do you define the Wiener path on a coarser subgrid? The transformation for the coarse increments ∆W are easy to compute by summing the finer increments, but the transformations for the multiple integrals are more complicated. Consider the multiple integrals over the time interval [0, t], which is divided into N subintervals of time ∆= t/N. Then the multiple integrals 18Kevin Burrage, Pamela Burrage, and Taketomo Mitsui, ‘‘Numerical solutions of stochastic differential equations — imple- mentations and stability issues,’’ Journal of Computational and Applied Mathematics 125, 171 (2000) (doi: 10.1016/S0377- 0427(00)00467-2). Actually, these authors give a more general transformation to two time intervals of nonequal duration.
27.8.1 Uniform Distribution¶
27.8 Random Numbers on the full interval in terms of the fine-step integrals may be written (Problem 27.4)
N−1 X n=0 ∆W(n∆t) (I10)t 0 = N−1 X n=0 h (I10)(n+1)∆t n∆t + (N −1 −n)∆t ∆W(n∆t) i (I100)t 0 = N−1 X n=0 h (I100)(n+1)∆t n∆t
n∆t + ∆t2 2 (N −1 −n)2∆W(n∆t) i (I1000)t 0 = N−1 X n=0 (I1000)(n+1)∆t n∆t
n∆t + ∆t2 2 (N −1 −n)2(I10)(n+1)∆t n∆t + ∆t3 3! (N −1 −n)3∆W(n∆t) (I11)t 0 = N−1 X n=0 h (I11)(n+1)∆t n∆t + ∆W(n∆t) W(n∆t) i (I110)t 0 = N−1 X n=0 (I110)(n+1)∆t n∆t + ∆t(N −1 −n) h (I11)(n+1)∆t n∆t + ∆W(n∆t) W(n∆t) i
n∆t W(n∆t) . (coarser Brownian path) (27.51) Here, the explicit time ranges are specified for each multiple integral for clarity, so that for example (I10)t2 t1 := Z t2 t1 dt Z t t1 dW(t′). (27.52) Then the transformations (27.50) follow from writing each refined variable as a linear combination of the coarse variables and the extra random variables, and enforcing consistency of Eqs. (27.51) as well as the proper variances and covariances of the refined variables. 27.8 Random Numbers The ability to generate quality ‘‘random’’ numbers on a computer is obviously an important part of simulating stochastic systems. We will thus spend a bit of time surveying some useful techniques towards these goals. 27.8.1 Uniform Distribution The workhorse of stochastic simulations is a random-number generator that generates numbers uniformly distributed on the interval [0, 1). These numbers can then be transformed to whatever distribution needed by various methods. Most compilers and canned software packages already include such a generator, so why bother thinking about it? First, you might not want to trust just any random-number generator that was handed to you, since there have historically been problems found with random-number generators, some subtle.19 Second, it is useful to have a portable random-number generator, which can generate the exact same sequence of ‘‘random’’ numbers under any compiler/architecture, which is impossible using ‘‘built- in’’ generators. Doing this generally requires implementing the generator in integer arithmetic to avoid floating-point rounding variations on different architectures. This is very handy when porting your code to 19Alan M. Ferrenberg, D. P. Landau, and Y. Joanna Wong, ‘‘Monte Carlo simulations: Hidden errors from ‘good’ random number generators,’’ Physical Review Letters 69, 3382 (1992) (doi: 10.1103/PhysRevLett.69.3382). The authors found that certain Monte-Carlo methods combined with certain random-number generators (which otherwise passed standard statistical tests) produced clear systematic errors in the fourth or fifth decimal places of calculated quantities in the 2D Ising model.
Chapter 27. Stochastic Differential Equations a new platform: the easiest test is just to run it and verify that you get essentially the same output for a particular run. This won’t happen if the random numbers are different, however. Finally, having several known methods available allows you to switch methods to make sure you aren’t seeing any artifacts due to a particular method. So just what is a random-number algorithm on a computer? Obviously, it is a deterministic algorithm, and ultimately, with the finite precision of any computer, the algorithm will be periodic. The idea is to come up with an algorithm that (1) produces the correct (uniform) distribution, (2) has a period much larger than the number of pseudorandom numbers needed for simulations, and (3) shows no detectable correlations over sets of pseudorandom numbers large enough for useful simulations. The algorithm must generally be chosen carefully to meet these criteria.20 Useful algorithms may be ‘‘seeded’’ with an initial value to produce the same set of random numbers each time. Also, it is desireable that different seeds correspond to different random-number sequences (that is, they should start the generator with initial conditions such that the sequences do not overlap for many iterations. Batteries of statistical tests are available to check the quality of random-number generators, such as Marsaglia’s ‘‘Diehard Battery’’21 and the battery of tests described by Knuth.22 We point out three algorithms here that meet all these criteria, plus some methods for improving the random numbers further.23 27.8.1.1 L’Ecuyer’s Multiple Recursive Generator A class of computationally very simple methods go by the name of linear congruential generators (LCGs),24 and implement the recurrence
(27.53) (linear conguential generator) in integer arithmetic. These algorithms, implemented with real numbers, would clearly be simple chaotic systems, where the state xn is stretched by a, shifted by c, and then folded by m. Of course, in integer arithmetic, the sequence is periodic, and the constants here must be chosen carefully to give decent per-
their simplicity, but they can have problems, such as a period that may be short in the context of modern processors, and they can have some problems with statistical tests25 A multiple recursive generator (MRG) improves on LCGs at the expense of added complexity by expanding the depth of the recurrence. An MRG of order k has the form
(multiple recursive generator) (27.54) L’Ecuyer’s combined multiple recursive generator26 uses two MRGs, where the first generator hsa coefficients
a11 = 0 a12 = 63 308 a13 = −183 326 b1 = 0 (27.55) 20For a particularly amusing discussion, see Donald E. Knuth, The Art of Computer Programming, Volume 2: Seminumerical Algorithms, 3rd ed. (Addison Wesley, 1998), Chapter 3. 21George Marsaglia and Wai Wan Tsang, ‘‘Some difficult-to-pass tests of randomness,’’ Journal of Statistical Software 7, No. 3, 1 (2002). 22Donald E. Knuth, op. cit., Section 3.3, p. 41. 23All of the algorithms described here are implemented in the Fortran 90 module/library RANDOM_PL, available online at http://atomoptics.uoregon.edu/~dsteck/computer.html. Also implemented are the shuffling and combination algorithms for any combination of the three generators, as well as facilities for checkpointing, generating vectors of numbers, generating Gaussian numbers, and running multiple equivalent (but differently seeded) generators in parallel. These algorithms are also implemented in the SDERK module mentioned above. 24see Donald E. Knuth, op. cit., Sections 3.2 and 3.3 for a detailed discussion of LCGs. 25Donald E. Knuth, op. cit., Sections 3.3.3 and 3.3.4. 26Pierre L’Ecuyer, ‘‘Combined Multiple Recursive Random Number Generators,’’ Operations Research 44, 816 (1996), specif- ically the generator in Example 4.
27.8 Random Numbers while the second has m2 = 2 145 483 479 a21 = 86 098 a22 = 0 a23 = −539 608 b2 = 0. (27.56) Some of the coeffiecients are zero as a compromise between speed and quality. The idea then is to run both generators simultaneously and combine their outputs via
n + x(2) n ) mod m1. (27.57) Since the first random-number generator is uniform over the positive range [0, m1−1], we can think of it being a uniform distribution on the circle. Adding the second random number corresponds to a random rotation of the circle, which doesn’t affect its distribution. This combination further breaks up any correlations or problems that might occur with each individual generator. This combined generator has a period of about 2185, and outputs integers in the range [0, 231 −2], which can then be divided by 231 to produce uniform random numbers on [0, 1). This generator can be seeded in an elegant way that guarantees that different seeds will produce sequences that don’t overlap those of other seeds for a maximally long time. First, six initial numbers are needed to start the recurrence, so we can just pick them to be some fixed but otherwise arbitrary, ‘‘random’’ numbers. For concreteness, we can choose x(1) = 1 852 689 663 x(1) = 1 962 642 687 x(1) = 580 869 375 x(2) = 2 039 711 750 x(2) = 1 671 394 257 x(2) = 879 888 250. (27.58) Now the idea is that with the MRGs, we can efficiently skip ahead to anywhere in the random sequence via a divide-and-conquer algorithm.27 To understand this, first note that the MRGs can be implemented as the linear transformation xn+1 xn xn−1 = a1 a2 a3 xn xn−1 xn−2 mod m =: A xn xn−1 xn−2 mod m (27.59) on the three-dimensional state of the generator. Then advancing the generator forward n times is equivalent to instead applying the matrix (An) mod m Then since the period is just slightly more than 2184, and if we seed with a single 32-bit integer, there are 232 possible seed values. Thus, we can start different seeds 2184−32 or 5.7 \times 1045 iterations apart. Thus if s is the seed value, we must apply the matrices (A2184−32)s mod m to the initial state of the random number generators, where s is possibly of the order of 232. This seems like an absurdly large amount of matrix multiplication just to get started! The trick behind the divide and conquer algorithm is to compute the large power of the matrix A recursively, according to
A
(AAn−1) mod m (n > 1, n odd) (An/2An/2) mod m (n > 1, n even). (divide-and-conquer algorithm) (27.60) 27This algorithm is described in depth by Pierre L’Ecuyer, ‘‘Random Numbers for Simulation,’’ Communications of the ACM 33, 85 (1990) (doi: 10.1145/84537.84555).
Chapter 27. Stochastic Differential Equations In this way, An may be computed in only O(log2 n) matrix multiplications. Specifically, (A2184−32) mod m may be computed in about 152 matrix multiplications. Then the matrix (A2184−32)s mod m needed to compute the seed can be computed in at most about 32 extra matrix multiplications, which is certainly quite feasible. 27.8.1.2 Knuth’s Lagged-Fibonacci Generator The second generator we mention is Knuth’s subtractive, lagged-fibonacci generator.28 This generator im- plements the integer recurrence
(lagged–Fibonacci generator) (27.61) which is obviously a multiple recursive generator. However, to improve the random numbers here, 1009 numbers can be generated at a time, but only the first 100 numbers used (the rest to be discarded). Knuth provides an algorithm for initializing the first block of values so that the recurrence can take over; this is done in such a way that each seed in the range [0, 230 −3 = 1 073 741 821] gives a distinct sequence for roughly at least the first 270 numbers (but take away a factor of 10 to account for discarding numbers). This generator has a period of 229 \cdot (2100 −1), or 1038 numbers, which is quite good (again, take away a factor of 10 to account for discarding numbers). The generator gives a 32-bit integer output in the range [0, 231 −2] (note that the least significant bit is always 0, so the number is always odd), which can then be divided by 231 to produce a uniform random number on [0, 1). 27.8.1.3 Mersenne Twister The Mersenne Twister generator of Matsumoto and Nishimura is rather more complicated than the other algorithms, being of much higher dimension. Thus we do not describe the algorithm, but refer to the original reference, where C source code is included.29 The reason we mention it is that it has an incredible period of 219937 −1, and so by this standard is an extremely good generator, while still being fairly computationally efficient. The state of this generator is represented by 624 32-bit integers, so the same number of initial values are required. This method can be seeded by a single number simply by using it as the first integer, and using a simple LCG method to seed the rest. The period is so long that it is exceedingly unlikely that two seeds will produce overlapping sequences on any reasonable scale. 27.8.1.4 Randomizing Random Numbers There are a few methods to try to improve the quality of random-number generators like the ones above. The first method is to ‘‘shuffle’’ or ‘‘scramble’’ the output of one random-number generator, based on the output of another.30 The idea is to maintain a ‘‘buffer’’ table of, say, 64 output numbers filled by the primary generator. Then a random number is chosen from the second generator, and based on this, one of the numbers in the buffer table is chosen to be the output number, and it is replaced by the primary generator. This procedure extends the period of the shuffled sequence to be the least common multiple of the two parent sequences, which for the above methods is to say incredibly long. An alternative is to to a similar shuffle of a random sequence using its own randomness to choose which element of the shuffle table to output. This is the Bays–Durham shuffle31 Then for example, the last output random number is used to select which element of the scramble table to output. Even though only one generator is involved, the method still improves the ‘‘randomness’’ of the parent generator. Both of these shuffling methods break up short-range correlations that may exist in the methods above. 28Donald E. Knuth, op. cit., p. 186, Eq. (2) and the following discussion. 29Makoto Matsumoto and Takuji Nishimura, ‘‘Mersenne Twister: A 623-Dimensionally Equidistributed Uniform Pseudo-Ran- dom Number Generator,’’ ACM Transactions on Modeling and Computer Simulation 8, 3 (1998) (doi: 10.1145/272991.272995). 30Donald E. Knuth, op. cit., p. 33, Algorithm M. 31Carter Bays and S. D. Durham, ‘‘Improving a Poor Random Number Generator,’’ ACM Transactions on Mathematical Software 2, 59 (1976) (doi: 10.1145/355666.355670); Donald E. Knuth, op. cit., p. 34, Algorithm B.
27.8.2 Gaussian Distribution¶
27.8 Random Numbers The other method for combining two random-number generators is the subtraction we mentioned above for the L’Ecuyer algorithm. Combining in a second generator can’t make the first one worse, but may help if it has some subtle defects. Combining three random number generators via subtraction and scrambling gives the ultimate security for the paranoid computational scientist. 27.8.2 Gaussian Distribution To pick random deviates from a standard normal distribution, or a Gaussian distribution with zero mean and unit variance, with probability density
\sqrt
(27.62) it turns out to be convenient to do so by generating two uniform deviates at a time, and then transform them to Gaussian deviates. The first algorithm we mention for doing this is the Box–Muller method,32 To understand it, we first write the joint distribution for two standard normal deviates as
(27.63)
(27.64) via the usual polar transformation. This is still separable, as we may regard the polar distribution to be the product of
2\pi ,
(27.65) where the first distribution is obviously uniform over [0, 2\pi). For the second, we may set
e−r2/2 , (27.66) then it is consistent to identify u = 1 − e−r2/2 , (27.67) or solving for r, r = p −2 log(1 −u). (27.68) Of course, 1 −u has the same distribution as u, so we may as well use it instead: r = p −2 log u. (27.69) What this means is that if we take u to be a uniform deviate and use this transformation to find r, then
function. So, to summarize the Box–Muller algorithm, we choose U1 and U2 to be uniform deviates on [0, 1), and set U = U1 and \Theta = 2\piU2. Then r = \sqrt−2 log u, and since x = r cos \theta and y = r sin \theta, then the two transformed numbers Z1 = p −2 log(1 −U1) cos(2\piU2), Z2 = p −2 log(1 −U1) sin(2\piU2), (Box–Muller transformation) (27.70) 32G. E. P. Box, Mervin E. Muller, ‘‘A Note on the Generation of Random Normal Deviates,’’ Annals of Mathematical Statistics 29, 610 (1958).
27.8.3 Other Distributions¶
Chapter 27. Stochastic Differential Equations are two independent, standard-normal deviates. Note here that we were careful to compute the logarithm of 1 −U1 rather than U1 here—if U1 is uniform on [0, 1) (as is the case for most random-number generators on computers) it is best to avoid the (rare but problematic) possibility of computing log(0), whereas log(1) is no problem at all. An variation on the Box–Muller method is the polar Marsaglia method.33 Suppose w1 and w2 define the coordinates for a uniformly distributed random number inside the unit circle. Then let R2 = W 2 1 + W 2 2 . From the usual transformation to polar coordinates, the probability density for R is 2R (the area element
identify the probability density for R2 as being uniform on [0, 1). Thus, we may transform R to the radius
Z1 = W1 r −2 log R2 R2 , Z2 = W2 r −2 log R2 R2 . (27.71) (polar Marsaglia method) The advantage here is that the sin and cos functions are eliminated, and they may potentially be compu- tationally expensive. Now the remaining issue is, how to pick W1 and W2? Simply pick the usual uniform
1 + W 2 2 . If R2 < 1, then proceed with the transformation; otherwise pick new deviates U1 and U2 and try again until you end up with a pair inside the unit circle. The probability on one attempt of succeeding by ending up in the unit circle is \pi/4 \approx 78.5%; on average, the
times to successfully generate the two normal numbers. In principle, the polar Marsaglia method should be faster than the Box–Muller method because the sin and cos functions can be expensive to evaluate. However, this is tempered by the presence of a division and the generation of extra random numbers. To test this, 108 normal deviates were generated on a relatively modern computer (2.16 GHz Intel Core Duo) in Fortran 90. Using Knuth’s lagged Fibonacci generator, the polar Marsaglia method is indeed faster by about 10-15%. However, using a more complicated generator (the L’Ecuyer combined MRG subtractively mixed with the Mersenne twister, then scrambled by Knuth’s generator), the two methods were about equally fast. Thus, either method is acceptable, and the Box–Muller method may even be preferable in that the generator will advance by a known number of iterations when generating normal deviates. 27.8.3 Other Distributions While not directly useful in simulations of the type of SDEs that we considered at the start of the chapter, the methods so far for generating random numbers are useful in generating random deviates with other distributions useful in quantum optics, so it is worth considering these briefly here. 27.8.3.1 Angular Distributions One useful class of random-deviate distributions are angular distributions, as for choosing random spontaneous- emission directions. A simple idea here is to simply use the rejection method, as we mentioned above in the polar Marsaglia method to convert a uniform variate on the unit box to a uniform variate on the unit circle. In that case, we simply rejected any deviate that fell outside the desired circle. Thus, suppose we have a distribution function f(\theta, \phi) that we wish to simulate. Our procedure will be a ‘‘double rejection’’ as follows: 1. Choose three uniform deviates U1, U2, and U3 on [0, 1). 2. Convert them to obtain a uniform deviate in the unit cube (that is, from [−1, 1) in all directions) by setting X = 2U1 −1, Y = 2U2 −1, and Z = 2U3 −1. 33G. Marsaglia, ‘‘Improving the polar method for generating a pair of random variables,’’ Boeing Scientific Research Labora- tory report D1-82-0203 (1962); Donald E. Knuth, op. cit., Algorithm P (p. 122).
27.8 Random Numbers 3. Obtain a uniform deviate in the unit sphere by rejection: Compute R2 = X2 + Y 2 + Z2 and continue
\sqrt R2); go back to step 1 otherwise. 4. Convert to spherical coordinates by computing \Theta = tan−1 \sqrt X2 + Y 2 Z \Phi = tan−1 Y X . (27.72) Note that for computing \Phi, the arctangent must be defined such that it can distinguish between
This is handled in Fortran and other languages by the atan2 function. 5. Obtain a deviate with the desired angular distribution by testing to see if U4 \le f(\Theta, \Phi)
(27.73) where U4 is another uniform deviate on [0, 1). To save from running the random-number generator again, we can alternately use R3 in place of U4, since it is likewise a uniform deviate on [0, 1), and it is independent of the angles. We normalize by the maximum possible value of f(\theta, \phi) to maximize the efficiency of the method. If the condition is true, we’re done; if not, start over at step 1. For example, if we were generating a dipole distribution
(27.74) in step 5 we would check to see if R3 \le sin2 \Theta, (27.75) and if so, keep the generated angles. In this example, the success probability for each triple of uniform deviates U1,2,3 to generate a uniform deviate in the unit sphere is \pi/6 \approx 52%. Once this happens, the probability of generating a deviate with the dipole distribution pattern is 2/3 \approx 67%, for a combined success probability of \pi/9 \approx 34.9%. The first probability is given by the volume ratio of the unit sphere and cube, while the second is given by the volume under the ‘‘unwrapped’’ surface sin3 \theta (i.e., not the volume inside
Of course, the normalized values X/R, Y /R, and Z/R give the projections of the unit vector pointing in the (\Theta, \Phi) direction along the x, y, and z directions. A slightly more efficient method generates the angle values more directly. The normalized probability density for equally-likely emission in any direction is
4\pi , (27.76) with normalization convention Z \pi d\theta Z 2\pi
(27.77) That is, the \phi variable is already uniform, so we can handle it by picking a uniform deviate u2 and rescaling: \phi = 2\piu2. Thus, we need only consider the \theta variable, with normalized density
2 sin \theta. (27.78)
this is equivalent to equating the cumulative probability functions. For the uniform deviate, the cumulative distribution is simply u1, while for \theta, the cumulative distribution is
2 −1
(27.79)
Chapter 27. Stochastic Differential Equations Equating the cumulative distributions, we obtain \theta = sin−1\sqrtu1. Thus, to summarize, we should choose two uniform deviates U1 and U2, and then compute the angles via \Theta = sin−1p U1,
(transformation from uniform to angular variables) (27.80) This saves one uniform random deviate compared to the rejection method (plus any rejected deviates), although it requires evaluating the more complicated arcsine function. 27.8.3.2 Exponential Distributions Another simple example is the exponential distribution,
(27.81) defined on the nonnegative half of the real line. The cumulative density is an exponential rise
(27.82) and the inverse of this function gives x in terms of a uniform deviate U: X = −1 \gamma log(1 −U). (27.83) (exponential-deviate generation) Of course, the result is always a positive real number. Also, U can be substituted for U in the logarithm, but if the random-number generator gives U \in [0, 1), the above form guards against the possibility of crashing the code by computing log(0). 27.8.3.3 Power-Law Distributions Suppose we have a power-law probability density defined between a and b:
0 < a \le x \le b. (27.84) We have chosen the bounds ‘‘safely’’ here, but in some cases they may be negative, too. Here, \eta is a normalization factor, given by
\gamma + 1
log b −log a
(27.85) Note that we can safely take a −\rightarrow 0 if \gamma > −1, otherwise we need a > 0. Similarly, we can safely take b −\rightarrow \inftyif \gamma < −1, otherwise we need b < \infty. If \gamma = 1, we strictly need 0 < a < b < \infty. First, in the case \gamma̸ = −1, we have a cumulative distribution
(27.86) Setting this equal to a cumulative distribution u for a uniform deviate, and solving, we find the algorithm X =
=
(power-law-deviate generation) (27.87) for a power-law deviate X in terms of a uniform deviate U. However, if \gamma = −1, then we have the cumulative distribution
log b −log a
(27.88)
27.8 Random Numbers This leads to the transformation X = exp log a + (log b −log a)U = exp (1 −U) log a + U log b
(power-law-deviate generation) (27.89) In both cases, the transformation is a weighted average of the endpoints (according to the value U), within a nonlinear transformation (power-law or logarithm). One context in which the power-law distribution comes up in the Monte-Carlo evaluation of an integral of the form I = Z \infty dx x\gammaf(x) (\gamma < −1), (27.90) where the restriction on \gamma is just to simplify the following analysis. Clearly, this integral only converges if f(x) converges to a finite constant as x −\rightarrow \infty, and if f(x) converges quickly enough to zero as x −\rightarrow 0. For the power-law distribution to do us any good here we must have a lower bound a > 0. We will thus
to the integral—we can always find some appropriate boundary that satisfies this condition to any desired precision). Then I = Z \infty xmin dx x\gammaf(x) (\gamma < −1), (27.91) and thus we can take advantage of the power-law distribution to turn this into an ensemble average,
f(x)
min \gamma + 1
f(x)
x (\gamma < −1), (27.92) where from Eq. (27.87), we generate the values of x according to
(\gamma < −1). (27.93) Note again that we could replace 1 −U by U to simplify the formula here, but it is safer if U is defined on
27.8.3.4 Gamma Distribution A more complicated distribution, something of a combination of power-law and exponential distributions, is the gamma distribution
(27.94) named thusly due to the gamma function in the normalization factor. In the case of \gamma = 1, this is simply an exponential distribution, as in Eq. (27.81). One useful interpretation of the gamma distribution comes from the property that if X1 ∼f(x; \gamma1, \beta) (i.e., X1 is gamma distributed with parameters \gamma1 and \beta) and X2 ∼f(x; \gamma2, \beta), then (Problem 27.2)
(27.95) Thus, for example, if \gamma is an integer, the gamma-distributed deviate X can be interpreted as the sum of \gamma, independent exponential deviates. This also means that as \gamma becomes large, f(x; \gamma, \beta) should converge to a Gaussian distribution via the central limit theorem. It is also fairly simple to choose a gamma-distributed deviate if \gamma = a is a relatively small integer, since the combination of exponential deviates X = −1 \beta a X j=1
\beta log a Y j=1 (1 −Uj) (gamma-distributed deviate, integer \gamma) (27.96)
Chapter 27. Stochastic Differential Equations works, for a uniform deviates Uj. The latter form is more efficient on a computer, since it avoids the calculation of multiple (expensive) logarithms. Another useful special case of the gamma distribution occurs for \gamma = 1/2:
r \beta
(x > 0, \beta > 0). (27.97) Changing variables via x = z2/2 and dx = z dz = \sqrt 2x dz leads to
r 2\beta
(z > 0, \beta > 0). (27.98) This is a Gaussian distribution (or at least the positive half) with variance 1/\beta. Thus, it is also easy to pick a gamma-distributed deviate for any half -integer value of \gamma. That is, if
2 Z2 −1 \beta a X j=1
2 Z2 −1 \beta log a Y j=1 (1 −Uj) , (gamma-distributed deviate, half-integer \gamma) (27.99) where a is a positive integer, Z is a standard-normal deviate, and the Uj are uniform deviates, then X is
By the same coordinate change, the \gamma = 3/2 case,
\sqrt\pi
(x > 0, \beta > 0), (27.100) becomes
r \beta3
(z > 0, \beta > 0). (27.101) Notice that this is the same as a Gaussian distribution of variance 1/\beta for three variables x1, x2, and x3, but transformed into spherical coordinates by z2 = x 2 1 + x 2 2 + x 2 3 . Thus, for \gamma = 3/2, it is possible to create a deviate by summing the squares of three normal deviates (with appropriate scaling), or to add the square of a normal deviate to an exponential deviate as in Eq. (27.99). For more general \gamma (i.e., neither integer nor half-integer), it is in general necessary to use more com- plicated methods, based on various forms of rejection sampling.34 34Luc Devroye, Non-Uniform Random Variate Generation (Springer, 1986), Section IX.3 (ISBN: 0387963057), http://www. nrbook.com/devroye/.
27.9 Exercises 27.9 Exercises Problem 27.1 Let us define ∆W and I10 over the same time interval ∆t as usual as ∆W := Z ∆t dW(t) I10 := Z ∆t dt Z t dW(t′). (27.102) Show that the variance of I10 and its covariance with ∆W are given by Eq. (27.16),
I 2
= ∆t3
I10 ∆W
= ∆t2 (27.103) as follows. View I10, defined over a time step ∆t, as the solution of the differential equation
(27.104)
becomes exact as N −\rightarrow \infty). Problem 27.2 Prove the addition property of the gamma distribution [Eq. (27.95)]: if X1 ∼f(x; \gamma1, \beta) (i.e., X1 is gamma distributed with parameters \gamma1 and \beta) and X2 ∼f(x; \gamma2, \beta), where the gamma probability density is
(27.105) then the sum is also gamma-distributed:
(27.106) Problem 27.3 A useful probability density for sampling the T integral in worldline path integrals is
\Gamma D+m e−2d2/T
(27.107) Show that this can be transformed to a gamma distribution, and describe the generation of deviates T for D = 4 and integer m \ge 0. Problem 27.4 Derive the equations (27.51) relating the multiple integrals on the full time interval [0, t] to the multiple
Chapter 27. Stochastic Differential Equations integrals defined on subintervals ∆t = t/N:
N−1 X n=0 ∆W(n∆t) (I10)t 0 = N−1 X n=0 h (I10)(n+1)∆t n∆t + (N −1 −n)∆t ∆W(n∆t) i (I100)t 0 = N−1 X n=0 h (I100)(n+1)∆t n∆t
n∆t + ∆t2 2 (N −1 −n)2∆W(n∆t) i (I1000)t 0 = N−1 X n=0 (I1000)(n+1)∆t n∆t
n∆t + ∆t2 2 (N −1 −n)2(I10)(n+1)∆t n∆t + ∆t3 3! (N −1 −n)3∆W(n∆t) (I11)t 0 = N−1 X n=0 h (I11)(n+1)∆t n∆t + ∆W(n∆t) W(n∆t) i (I110)t 0 = N−1 X n=0 (I110)(n+1)∆t n∆t + ∆t(N −1 −n) h (I11)(n+1)∆t n∆t + ∆W(n∆t) W(n∆t) i
n∆t W(n∆t) . (27.108)
Index A gauge, 462 E gauge, 462 P function, 129–130 Q function, 129 T1 and T2, 178 \sigma-algebra, 869 It¯o–Taylor expansion, 1169–1172 It¯o calculus, 948–953 Abraham–Lorentz force, 39, 451 absolute-value process, 776–777, 802 absorption, 98 absorption coefficient, 57, 102–103 absorption image, 42 absorption oscillator strength, 35–36 ac Stark shift, 162, 272, 273, 397 due to thermal field, 508 action functional, 403, 910, 914 expansion, 914, 915 action principle, 403–404 for Schrödinger equation, 452 optics, 451 Adams–Bashforth method, 1131 Adams–Moulton methods, 1131 adaptive measurement, 863 adaptive time stepping, 1130 adiabatic approximation, 104, 191, 230–231, 234, 236, 271–273, 528–530 for momentum, 734–737 via master equation, 230–231 adiabatic passage, 162–166, 177 adiabatic rapid passage, 166 adiabatic theorem, 163 advanced Green function, 654, 655 advection equation, 239 affine parameter, 993 Aharonov–Bohm effect, 465–467 aliasing, 1139–1140 ancilla, 143 angular momentum, 301–315 anomalous diffusion, 705 antibunching, 76, 834 resonance fluorescence, 209–211 anticipating calculus, 953–954 antistandard-ordered distribution, 130 arcsine law, 775–777 associated Legendre functions, 412 atom two-level, 151–161 atom interferometry, 281–283 atom optics, 39–49 dipole force, 39–45 Doppler temperature, 47–49 magneto-optical trap (MOT), 49 optical molasses, 46–49 radiation pressure, 45–49 atom–mirror interaction, 50–53, 57–58, 507, 509, 637–638 atom–photon molecule, 490–491 atomic clock, 171–174 Autler–Townes doublet, 217–221, 256 autocorrelation function, 64 auxiliary functions to sine and cosine integrals, 537 averaging operator, 958 avoided crossing, 161 backward differentiation formulae, 1131 backward Euler method, 1125 Baker–Campbell–Hausdorff expansion, 124, 908, 1149 bare states, 161 Bays–Durham shuffle, 1184 beam splitter relation among coefficients, 512, 687 Bernoulli numbers, 631 Bessel process, 807–811 Bessel’s equation, 413, 455 big-O notation, 1121–1123 Blackman pulse, 277–278 bleaching, 103 Bloch sphere, 166–177 Bloch vector, 168–178 Bloch–Siegert shift, 42, 251–252 Boltzmann statistics, 99 Borel \sigma-algebra, 870 Born approximation, 145 Born rule, 869 Born–Markov approximation, 145 bound charge density, 579, 1010
Index bound current density, 580 boundary crossing, 748–754, 758–761 Box–Muller method, 1185–1186 Bragg scattering, 232–235 rate, 235 branch prediction, 1115–1116 Breit–Rabi formula, 397 Brownian area, 732 Brownian Bridge, 827–828 Brownian bridge, 737–807, 913, 925, 927, 930, 931, 961, 980, 982, 983, 985, 986 boundary crossing, 750–754, 758–761, 985, 986 escape, 754–758 escape probability, 983 local time, 797–807 sojourn time, 777–807 Brownian bridges, 988 Brownian excursion, 738 Brownian meander, 738 Brownian motion, 719–720, 925 reflected, 776–777, 802 skew, 802–807 Brownian path consistent refinement, 1177–1181 bunching, 73–76 cache, 1111–1112 line, 1111 thrashing, 1112 Cameron–Martin transformation, 949 Cameron–Martin–Girsanov transformation, 949–950 canonical cosymplectic form, 135–137 canonical quantization, 409 canonical transformation, 136–137, 463 generating function, 463 Casimir effect path integral, 965–989 Casimir energy, 638–647 temperature dependence, 646–647 Casimir–Polder potential, 42, 464, 533–557, 563–578, 596–635 atom–atom potentials, 616–621 comparison to classical dipole, 543–544 dielectric and general media, 596–635 excited levels, 632–635 general planar interface, 601–614 Lifshitz expression, 608 path integral, 965–989 rarefied dielectric, 563–573, 610 spherical cavity, 614–616 temperature dependence, 621–632, 988–989 three atoms, 620–621 Cauchy integral formula, 589, 601, 628 Cauchy principal value, 289, 499, 542, 590, 650 Cauchy principle value, 543 Cauchy probability distribution, 704, 710 Cauchy’s integral formula, 635 caustics, 921 cavity damping, 492–493 cavity decay, 146, 491–492, 509–519, 531 decay rate, 516–519 nonexponential, 682 cavity driving field, 491 cavity photon number measurement of, 492–493 cavity QED coupling constant, 486 center-of-mass coordinates, 475–477 central limit theorem, 697–706 central-limit theorem, 711 Chapman–Kolmogorov equation, 717, 908 characteristic function, 702 chronological operator, 112 Church of the Larger Hilbert Space, 142–143 Clebsch–Gordan coefficients, 304–309 cloning, 131–132 closed set, 870 CNOT gate, 132 coarse-graining approximation, 145 coherence, 63–76 length, 70–72 time, 70–72 coherent population trapping, 273, 283–295 velocity-selective, 286–287 coherent state, 128, 187, 491 collapses and revivals, 487–489 collisions dephasing due to, 178 commutator with space-dependent mass, 934 composite systems, 130–131 conditional probability density, 716 Condon–Shortley phase convention, 307 conductivity, 59 conjugate momentum, 404 constitutive relations, 1009 continuous function, 1122 controlled-NOT gate, 132 convective derivative, 457 convergence definition, 1121–1123 convolution, 698–706, 709–710 of box functions, 699 of two Gaussians, 700–701 with \delta-function, 698–699 convolution kernel, 699
Index convolution theorem, 699–700, 1133 cosine integral Ci(x), 537, 574–575 cosymplectic form, 135–137 Coulomb gauge, 401, 422–423, 459 counters, hardware performance, 1119 counting process, 728, 833–835 coupled-mode theory, 511–519 critical sampling, 1134–1135 cross section absorption, 35, 101 natural, on resonance, 101, 102, 192 scattering, 195–197 crossover resonance, 225–228 D line, 397 dark state, 285 decay rate modified by macroscopic body, 50–53, 635–638 decay rate, natural, 271 degree of coherence, 94 degree of first-order temporal, 694 first-order, 66 first-order quantum, 686 first-order temporal, 66 second-order quantum, 686 second-order temporal, 72–76 delta function chain rule, 706, 828–829 derivative of, 794, 828–829 longitudinal, 423–425 transformation for surface constraint, 706–710 transverse, 423–425 density matrix, 109–110 path integral, 928–930 density operator, 107–110 purity, 110, 255 derivative process, 820 detector efficiency, 839–840, 853–855, 878–879 diamagnetic interaction, 473, 475 Dicke superradiance, 298–299 dielectric constant, 581 diffusion spontaneous emission, 247 diffusion coefficient, 239, 698 diffusion equation, 239 diffusion process, 697 diffusion rate discontinuity, 802–807 dipole approximation, 151, 460 dipole force, 39–45, 237–238, 254 dressed-atom picture, 256–258 dipole interaction Hamiltonian, 151–157, 457, 460–465, 485–486, 495 dipole matrix element relation to decay rate, 192 dipole moment electric, 474 magnetic, 474 dipole operator, 152–153, 229 relation to Bloch vector, 170 direct product, 130–131 direction cosines, 425–426 Dirichlet problem Laplace equation, 758–760 Poisson equation, 760–761 discrete Fourier Transform, 1136–1147 discrete spectrum, 488 dispersion, 580–581 divergence theorem, 706 divergent sequence, 1121 divide-and-conquer algorithm, 1183–1184 domination of sequences, 1121 Doppler broadening, 57, 100, 174 Doppler limit, 250–251 Doppler shift, 259, 273, 280 Doppler temperature, 47–49 Doppleron resonance, 224 dressed states, 161–166, 254, 489 drift coefficient, 239 Drude–Lorentz model, 59 dynamical tunneling, 234 effective potential, 150 effective, non-Hermitian Hamiltonian, 185 Einstein A and B coefficients, 98–100, 104–105 Einstein A coefficient, 178 Einstein rate equations, 98–105 EIT, 287–290 electric dipole moment, 474 electric displacement, 435, 579 electric displacement field, 462 electric flux density, 435, 579 electric octupole moment, 474 electric quadrupole moment, 474 electromagnetic field commutators, 425–430 gauge freedom, 431–434 Hamiltonian structure, 402–405, 431–434 mass coupling, 455–456 quantization, 406–420 electromagnetically induced transparency, 287–290 as position measurement, 899–904 energy spectral density, 64 entanglement, 131–140, 142–143 environment, 143, 144
Index error tolerance, numerical, 1130 escape probability, 754–758, 930 Euler method, 1123–1125, 1167 implicit, 1123–1125 stochastic, 1167–1168 Euler’s constant, 537 Euler–Cromer method, 1162–1163, 1165 Euler–Lagrange equation, 403, 405, 451–452, 458 Euler–Maclaurin summation formula, 631 Euler–Maruyama method, 1167–1168 Euler–Mascheroni constant, 575, 671 exchange interaction, 141–142 exponential integral E1(x), 575, 667, 681 Fabry–Perot cavity, 261, 296, 490, 517, 519, 755 Fano profile, 286 far-off resonance trap (FORT), 42 Fermat’s principle, 451 Fermi’s Golden Rule, 509–510, 531 Fermi’s golden rule, 500–502, 635, 637 Feynman–Kac formula, 761–770, 923, 927–928 Feynman–Vernon–Hellwarth representation, 169 fidelity, 117 finesse, 517 finite-difference method, 1123 finite-difference operator, 916 first variation, 402, 406, 914 first-passage time Brownian bridge, 752–754 Wiener path, 750 fixed-point arithmetic, 1109 fixed-point iteration, 1124 floating-point arithmetic, 1109–1110 fluctuation–dissipation relation, 48, 267, 522, 621–626 fluorescence, resonance, 194–211 Fock state, 187 Fock—Schwinger proper time, 975 Fokker–Planck equation, 238–240, 258–259, 716–718, discontinuity in, 802–807 multidimensional, 718, 826 forbidden transition, 343 forward Euler method, 1125 Fourier Transform discrete, 1136–1147 Fourier transform, 1133–1147 free charge density, 1010 free spectral range, 517 free-electron gas, 58–60 Fresnel reflection coefficients, 564, 603 FTIR spectroscopy, 67, 94 functional, 402, 406 functional derivative, 402–406, 452, 914 second, 914 functional determinant, 915, 971 as eigenvalue product, 973–975 functional integral, 907–962 fundamental charge, 152 Göppert-Mayer transformation, 462–463 gain coefficient, 105 gamma distribution, 1189 gauge transformation, 461, 463–467 Gaussian fluctuations complex, 74–76 exponentiated, 95 Gel’fand–Yaglom method, 916–918 general Legendre equation, 412 generalized Rabi frequency, 159 generalized susceptibility, 592–596 generating function for canonical transformation, 463 generator for SDE, 717 geometric Brownian motion, 826 geometric optics, 451 Girsanov transformation, 949 Glauber–Sudarshan P function, 129–130 Green function, 648–649, 653–655 classical harmonic oscillator, 648–649 for Schrödinger equation, 653–655 Green tensor, 563, 581–586 for free space, 583–586 Kramers–Kronig relations, 596 planar interface, 603–604 related to field correlation, 595, 596 scattering part, 601 Hamilton’s equations, 404, 405 Hamilton’s principle, 403 Hamiltonian, 404–406 Hanbury-Brown–Twiss experiment, 75–76, 685, 688 hardware performance counters, 1119 harmonic function, 759 harmonic oscillator, 31–36 coherent state, 187 damped, 33–36, 150 Fock state, 187 Green function, 648–649 quantum, 186–191 quantum damped, 187–191, 261, 266–267 heat bath, 144 Heisenberg picture, 112–114 Hellmann–Feynman theorem, 640 Helmholtz equation, 407 Green function, 681, 1002–1003
Index scalar vs. vector solutions, 414–418, 454 Helmholtz theorem, 421–422 Hermite–Gaussian functions, 121 heterodyne detection, 851–853, 859–863 heterodyne spectroscopy, 79–81 Hilbert transform, 289 homodyne detection, 840–851, 858–859, 861–863, 900 electromagnetically induced transparency, 899–900 Hong–Ou–Mandel dip, 694 Hong–Ou-Mandel effect, 688–694 Husimi distribution, 129 upper bound, 149 hydrogen maser, 174 hyperfine structure, 551–553 IEEE 754 standard, 1109–1110 impedance of vacuum, 40 implicit Euler method, 1123–1125 stochastic, 1174–1175 indistinguishability, 140–142 inhomogeneous broadening, 174 inlining, procedure, 1117 input–output theory, 519–531 integer arithmetic, computer, 1109 integrated Brownian motion, 732 intensity of Poisson process, 726 interaction picture, 114, 515 for atom–field interaction, 156–157 interference between partially coherent sources, 72 coherence, 63–76 two-photon, 688–694 visibility, 69–70 interferometer Michelson, 66–67, 94 Itō calculus, 712–726 ensemble average, 714–715 Itō rule, 713–714 Jaynes–Cummings model, 485–491, 515 as model for spontaneous emission, 487–489 Johnson noise, 624–626 joint probability density, 717 Jordan–Pauli commutators, 426–428 jump process, 728 Kirkwood distribution, 130 Klein–Gordon equation, 998–1002 Knuth lagged Fibonacci generator, 1184 Kolmogorov backward equation, 718, 769–770 Kolmogorov forward equation, 718 Kramers–Heisenberg formula, 596–600 Green tensor, 599–600 polarizability, 596–599 Kramers–Kronig relations, 289, 589–592, 648–649 Green tensor, 596 Kraus operator, 868 L’Ecuyer multiple recursive generator, 1182–1184 Lévy flight, 280 Lévy’s arcsine law, 775–777 ladder operators, 302 ladder structure, 232 lagged Fibonacci generator, 1184 Lagrangian, 403–406 for string, 451–452 Laguerre–Gaussian functions, 121 Lamb dip, 221–228 Lamb shift, 42, 464, 465, 499, 508–509, 533, 536, 557–563, 663, 982 Lamb–Dicke effect, 507 lamp shift, 162 Lamperti transform, 747–748 temporal scaling, 748 Landau–Zener crossing, 162–166, 254–255, 298 Laplace equation, 422, 758–760 Laplace transform, 180–181, 656, 681, 772, 1047–1050 Laplace–Beltrami operator, 935, 937 Laplace–Beltrami ordering, 933 Laplacian as spherical average, 759–760 spherical coordinates, 411 Larmor precession, 168 laser Lorenz–Haken model, 261–266 laser medium three-level, 104–105 Lebesgue measure, 870 level-shift operator, 659 Levi-Civita symbol, 301, 395, 582 light thermal, 74–76 likelihood-ratio estimator for derivative, 813 limit from above, 1122 from below, 1122 of function, 1122 of sequence, 1121 Lindblad form, 185, 866 Lindblad superoperator, 185, 187 interference represented by, 292 line shape, 100 natural, 101
Index linear congruential generator, 1182 linewidth measurement of, 78–93 Liouville distribution, 127 Liouville–von Neumann equation, 108, 180 Liouvillian superoperator, 108, 180 little-o notation, 1121–1123 local interactions, 465–467 local realism, 131 local time, 793–807 long-wavelength approximation, 151, 460 longitudinal decay rate, 178 loop unrolling, 1114–1115 Lorentz force, 457–458 Hamiltonian, 458 Lagrangian, 458 Lorentz model, 33–36, 186–189 damping coefficient, 38–39 Lorentzian absorption, 35 line shape, 100, 101, 665 lineshape, 68 noise spectrum, 78, 80–85 Lorenz gauge, 401 Lorenz–Haken model, 261–266 low-pass filter, 720 Lévy area, 1178–1180 magic wavelength, 174, 252–253 magnetic dipole moment, 474 magnetic flux density, 435, 579 magnetic-dipole transition, 397 magnetizability, 1018 magnetization field atomic, 468–469, 473, 480 magneto-optical trap (MOT), 49 Malliavin derivative, 819 Markov approximation, 145, 521 Markovian evolution, 717 martingale, 715 master equation, 143–148, 185, 187, 508 atom in thermal field, 508 Born–Markov, 143–148, 508 Lindblad form, 866 stochastic, 833–863, 865–868, 875–904 unconditioned, 833 matrix elements of r and p, 463–464 of atom–field interaction Hamiltonian, 464 Matsubara frequencies, 628 maximum-value process, 776 Maxwell equations, 401–402, 579–580 measurable set, 870 measure, 869–871 positive, 870 measurement record, 835 mechanical effects of light, 228–253 Mellin transform, 1047–1049 Mersenne Twister, 1184 methane-stabilized He–Ne laser, 225 method of images for diffusion, 805 Michelson interferometer, 66–67, 94 midpoint method, 1126 Milstein method, 1168–1169 minimal coupling, 458–460 Hamiltonian, 460 replacement, 460 mixed state, 107 Mollow triplet, 201–209, 255–256 nonclassical correlations, 208–209 moment-generating function, 1048–1049 moments inverse, 1048–1049 momentum-shift operator, 232, 272 Moore machine, 1115 Moyal bracket, 127–128, 150 multiple recursive generator, 1182–1184 multipole expansion, 473–475, 583 of localized current, 583–584 multipole fields, vector, 418–420 multistep methods, 1131 Naimark’s theorem, 872 natural line width, 102 Neumark’s theorem, 872 Newton’s Second Law, 404 NIST-7 cesium-beam clock, 174, 255 NIST-F1 cesium fountain clock, 174 no-cloning theorem, 131–132 nonexponential decay, 503–506, 666–675, 679–680 cavity, 682 normal modes of atom–cavity system, 489 normal-mode decomposition, 406, 972 octupole moment electric, 474 Ohm’s law, 625 one-photon Rabi frequency, 486, 487, 496 open quantum systems, 142–148 open set, 870 optical Bloch equations, 177–186, 189, 191–194, 230, 275, 500 Torrey’s solutions, 181–185, 255
Index optical lattice accelerating, 259 optical molasses, 46–49 optical theorem, 43–44 optics geometric, 451 Ornstein–Uhlenbeck process, 718–720 correlation function, 718–719 low-pass filter and, 720 oscillator strength, 35–36 oscillatory sequence, 1121 out-of-order execution, 1114 Paley–Wiener theorem, 506 parabolic-cylinder functions, 165 parallel programming, 1118 parametric downconversion, 688 parity, 152–153 Parseval’s theorem, 72 partition function, 1002–1003 path integral, 922–931, 936 path integral, 907–963, 965–989 It¯o, 941–944, 946–950 ensemble-average normalization, 923–927, 960–962 midpoint, 938–941, 946–950 moment relations, 944–946 Monte–Carlo method, 922–933, 960–962 operator-ordering issues, 933–962 partition function, 922–931 phase-space, 910–911 prepoint, 941–944, 946–950 product ordering, 963 Stratonovich, 938–941, 946–950 Pauli operators, 166 Peres–Horodecki criterion, 132–140 perimittivity imaginary frequencies, 591–592 permittivity, 581, 586–592 imaginary part as loss, 587–589 phase diffusion, 826–827 phase estimation, 861–863 phase lag of circuit response, 650 phase modulation, 280, 296 phase noise, 76–95 phase space, 115 phase-space action, 404 photon, 156 photon blockade, 489–490 photon echo, 174–177 photon scattering rate, 42–45 photon-recoil momentum, 272 pipeline, cpu, 1113–1116 Planck blackbody distribution, 100 plasma frequency, 58 plasma model, 58 Poincaré cone condition, 758 Poisson bracket, 127, 128, 135–137, 934 Poisson distribution, 726–728, 826 Poisson equation, 421, 760–761 Poisson process, 726–731, 826, 839, 843 inhomogeneous, 728 intensity, 726 white-noise limit, 728–729 Poisson sum rule, 931 polar Marsaglia method, 1186 polarizability, 32–33, 255 diagonal in principle coordinates, 595 related to dipole correlation, 594, 595 scalar, 595 static, 542 polarization field, 33, 579, 581 atomic, 461, 467–468, 470, 484 coupling to electric field, 60–61 effective sources, 579–580 polarization sum, 425–426 polarizer, 507–508 pole approximation, 660, 662, 663 position measurement electromagnetically induced transparency, 899–904 probe excitation, 895–899 resonance fluorescence, 885–895 positive map, 865–866 positive measure, 870 positive partial transpose (PPT), 132–140 positive semidefinite, 138 positive semidefinite operator, 133 positive-operator-valued measure, 868–875 imprecise measurement, 873–875 spontaneous emission, 872–873 POVM, 868–875 imprecise measurement, 873–875 spontaneous emission, 872–873 power broadening, 197 power spectral density, 65 one- vs. two-sided, 68–69 Power–Zienau transformation, 460–465, 467, 469–484, 544–548 power-equivalent width, 71 predictor–corrector methods, 1131 probability density for combinations of random variables, 709–710 probability measure, 871 probe absorption driven, two-level atom, 211–217, 259–260 problem:magdipdecay, 509
Index procedure inlining, 1117 product ordering, 959–963 profiler, 1118 projection theorem, 336–337 projection-valued measure, 869 propagator, 655, 907–962 classical limit, 913–922 free particle, 911–913 imaginary-time, 922–933 semiclassical, 913–922 Purcell factor, 682 pure state, 107 purification, 142–143 purity, 149–150 quadratures, 857–863 quadrupole moment, 484 electric, 474 quantum beats, 291–295 steady-state, 293–297 quantum jumps, 833–840 quantum Langevin equation, 521 quantum measurement, 868 continuous, of position, 875–904 strong, 874 weak, 874 quantum regression theorem, 199–201, 209–210, 213–214, 218, 241–243, 256, 260–261 quantum Zeno effect, 505, 507–508, 678 qubit, 107 quotient of two normal deviates, 710 Röntgen interaction center-of-mass, 475–482 Rabi flopping, 157–162, 487, 657–658 damped, 180–185 resolvent approach, 657–658 Rabi frequency, 153–154, 185–186, 229, 230, 270, 487 dependence on photon number, 487 for quadrupole interaction, 484 generalized, 230 Raman, 233, 272 radiation atom near macroscopic body, 635–638 atom near mirror, 50–53 two atoms, 53–56 radiation pressure, 45–49, 236 radiation reaction, 38–39 Abraham–Lorentz force, 39 radiation, dipole, 36–39, 194–195 angular distribution, 38 Radon–Nikodym derivative, 950 Raman scattering spontaneous, 675–683 Ramsey fringes, 171–177, 255 random number generation, 1181–1190 angular distributions, 1186–1188 exponential distribution, 1188 gamma distribution, 1189–1190 power-law distribution, 1188–1189 random process complex Gaussian, 74–76 random walk, 48, 697–706 random-number generation Gaussian deviates, 1185–1186 normal deviates, 1185–1186 shuffling, 1184–1185 uniform deviates, 1181–1185 rate equations Einstein’s, 98–105, 178, 191–194 Rayleigh’s formula, 418, 419 recoil energy, 233 Raman, 273 recoil frequency, 234, 235 recoil temperature, 49 recoil velocity, 45 reduced matrix element, 334 reflectance, 441 reflected Brownian motion, 776–777, 802 reflection coefficient, 58, 602–603 Fresnel, 564, 603 Reflection Principle, 749–751 registers, cpu, 1110–1111 regularization dimensional, 1028 reparameterization independence, 992 reservoir, 143, 144 residue, 589 resolved-sideband Raman cooling, 278–281 resolvent of the Liouvillian, 181 resolvent operator, 653–681 resonance fluorescence, 194–211, 255–256 as position measurement, 885–895 coherent and incoherent, 197–199 retarded Green function, 654, 655 Richardson extrapolation, 1151–1154, 1156–1158 Riemann sheet, 666 Rihaczek distribution, 130 rotating frame, 154–157, 254, 270–271, 283 rotating-wave approximation, 34, 153, 155, 157, 515, round-trip time, 511 rubidium 87 transition wavelengths and lifetimes, 553–554
Index Runge–Kutta methods, 1125–1131 adaptive stepping, 1130 implicit, 1129–1130 second-order, 1126 stochastic, 1172–1176 Runge–Kutta–Fehlberg method, 1130 sample mean, 703 sampling theorem, 1133–1136 saturated-absorption spectroscopy, 221–228 crossover resonance, 225–228 saturation, 102–103 saturation intensity, 103, 105, 192–194 for two-level atom, 41 saturation parameter, 179 saturation spectroscopy, 221–228 crossover resonance, 225–228 scattering cross section, 195–197 scattering rate, photon, 42–45 Schrödinger equation with curvature, 935 Schrödinger picture, 112–114 Schrödinger–von Neumann equation, 108, 177 Schrödinger-cat state, 133 Schwinger proper time, 975 second definition of, 174 second quantization, 409, 972–973 second variation, 914 selection rules, 342–344 self-heterodyne spectroscopy, 81–93 separability, 131 shot noise, 729–731 sine integral Si(x), 537, 574–575 single-precision, IEEE, 1109, 1110 singular D function, 427 skew Brownian motion, 802–807 slow light, 290 slowly varying coherences, 167 Smoluchowski equation, 735 Snell’s Law, 565 sojourn time, 770–807 spectral density of frequency fluctuations, 76 of frequency fluctuations(, 94 of frequency fluctuations), 95 of phase fluctuations, 76 of phase fluctuations(, 94 of phase fluctuations), 95 spherical Bessel functions, 413 spherical cavity modes normalization, 454–456 scalar, 411–414 vector, 414–418 spherical Hankel functions, 419 spherical harmonics, 412–413 spin echo, 174–177 split-operator methods, 1149–1165 spontaneous decay rate, 186, 498–499 between planar mirrors, 509 near planar mirror, 507, 637–638 with angular-momentum degeneracy, 344–346 spontaneous emission, 98, 146, 246–247, 495–506, 508–509, 662–675 enhanced by cavity, 682 Fermi’s golden rule, 502–503 in stimulated Raman scattering, 273 in thermal field, 508–509 into 1D waveguide, 509 into 2D waveguide, 509 line shape, 663–665 magnetic-dipole, 509 master equation, 499–500 near macroscopic body, 635–638 nonexponential decay, 503–506, 666–675 resolvent operator, 662–675 Weisskopf–Wigner approach, 495–499, 507, 509 spontaneous Raman scattering, 675–683 line shape, 682–683 square-well potential, 930 squeezed state, 454 squeezing, 858 Stückelberg angle, 161, 489 standard deviation, 697, 701 of mean, 703 standard topology, 870 standard-ordered distribution, 130 Stark shift ac, 162 due to thermal field, 508 stationary noise process, 88 stationary-phase approximation, 913 steady-state quantum beats, 293–297 stiff equation, 1124 stimulated emission, 98 stimulated Raman adiabatic passage, 290–291, 298 stimulated Raman cooling, 278–281 stimulated Raman transition, 232, 269–283, 298 velocity selection, 274–280 STIRAP, 290–291, 298 stochastic derivative, 819 stochastic differential equation, 712–726, 1167–1181 implicit methods for, 1174–1176 stochastic force on damped particle, 731–737 stochastic master equation, 833–868, 875–904
Index diffusion form, 840–851 jump form, 833–840, 864 linear, 881–885 multiple observers, 853–855, 878–879 position measurement, 865–904 stochastic Schrödinger equation, 835–836, 847, 864, 875–877 stochastic Taylor expansion, 1169–1172 stochastic-dipole force, 238–246 dressed-state model, 244–245 Stokes relations, 512 Stratonovich calculus, 720–726, 948–950, 954–957 validity of chain rule in, 724–725 strong convergence SDE methods, 1176–1177 strong coupling, 491 strong measurement, 874 Struve function, 276 subradiance, 56, 298–299 superradiance, 52, 56, 298–299 survival probability, 517 susceptibility, 33, 435, 581, 1009 generalized, 592–596 magnetic, 1009 symmetrized ordering, 957–959, 959, 963 symplectic integration, 1161–1165 symplectic matrix, 136–137, 149 Tanaka formula, 802 Taylor expansion stochastic, 1169–1172 tensor product, 130–131 thermal field, 508 thermal light, 74–76 thermal state Wigner distribution, 149 Thomas–Reiche–Kuhn sum rule, 36, 539–540, 543, 549, 550 time slicing, 911 topology, 870 trace, 109 transmittance, 441 transverse decay rate, 178 truncation error, 1122 tunneling, 254 two-level atom, 151–161 two-mode squeezed state, 140 two-photon absorption, 297–298 two-photon interference, 688–694 uncertainty relation, 70, 72 generalized, 135–137 unit-stride access, 1111 unitary time-evolution operator, 111–112, 907 unitary transformation, 155–156 vacuum Rabi doublet, 489 van der Waals–London potential, 616–621 Van Vleck–Morette determinant, 917–921 harmonic oscillator, 921–922 variation process, 820 variational calculus, 402–406 vector spherical harmonics, 415 Verlet method, 1163–1165 visibility, 69–70, 94 two-detector, 94 Voigt profile, 100 VSCPT, 286–287 wave equation for string, 451–452 weak convergence SDE methods, 1176–1177 weak measurement, 874 Weber’s equation, 164 Weisskopf–Wigner approximation, 495–499, 517–518 Weyl correspondence, 122–128, 149 Weyl ordering, 125–126, 933, 954–957, 963 Weyl’s quantization rule, 123 white noise, 712 Wick rotation, 922, 970 Wiener differential, 712 Wiener increment, 712 Wiener path reflected, 776–777, 802 Wiener process, 710–712, 737–811, 925 boundary crossing, 749–754, 758–761 correlation function, 715 escape, 754–758 local time, 793–807 recurrence, 807–811 sojourn time, 770–807 Wiener–Khinchin theorem, 64–69, 715 optical, 65–69, 194 Wigner 3-j symbol, 309 Wigner distribution, 115–129 area theorem, 117–118, 134 discrete Wigner transform, 1145–1147 Gaussian, 118–119 harmonic oscillator eigentstate, 121–122 superposition state, 119–121 thermal state, 149 upper bound, 149 Wigner–Eckart theorem, 333–337, 341–342 world line, 965 world lines, 980