Page 356 - Numerical Methods for Chemical Engineering

Brownian dynamics and SDEs 345

2
Substituting for dX ,
2
2
2
2
2
dX = [adt + bdW t ] = a (dt) + 2ab(dt)(dW t ) + b (dW t ) 2 (7.161)
While the first two terms are of higher order than 1 in time, the last term is not,as dW t ∼
√
2
dt. Thus, we must retain it, and replace it with its “average” value (dW t ) → dt.
We provide here a heuristic argument for accepting this replacement. Consider approxi-
mating the Wiener process as a random walk in which in one unit of time we make n steps
2
each of length l. After one unit time interval, the mean-squared displacement is nl = 1.
Thus, we obtain the same mean-square displacement if we replace the Wiener process by
2
2
a random walk with l = n −1 = dt. Thus, the replacement (dW t ) → dt is valid, in the
“mean-square sense”.
2
2
2
2
Using dX ≈ b (dW t ) → b dt in (7.160), we have Itˆo’s lemma

2

∂F ∂F 1 ∂ F 2 ∂ F
dF = + a(t, X) + [b(t, X)] dt + b(t, X)dW t
∂t (t,X) ∂ X (t,X) 2 ∂ X 2 (t,X) ∂ X (t,X)
(7.162)
This demonstrates the major difference between the stochastic and deterministic forms of
2
calculus. In stochasticcalculus, weexpandfunctionsto higherorders, replace (dW t ) → dt,
and then keep all terms that contribute up to the desired order.
We now use this formalism to demonstrate the derivation of a higher order integration
method than the explicit Euler one. In the explicit Euler method we neglect the time-variation
of a and b over the time step. This is particularly bad for the second integral as dW t is of
order t 1/2 , and thus the explicit Euler method is only 1/2-order accurate for predicting the
actual trajectory. Thus, let us increase the order of accuracy of this term by using a t 1/2
accurate expansion of b in t k ≤ t ≤ t k+1 :

∂b ∂ X
) + ] (7.163)
b(t, X t ) ≈ b(t k , X t k [W t − W t k
∂ X ∂W
(t k ,X t k ) (t k ,X t k )
Using ∂ X/∂W = b, we have for the second integral of (7.155),
(
' '
t k−1 t k+1
∂b
b(t, X(t))dW t ≈ b(t k , X t k ) + b(t k , X t k )[W t − W t k ] dW t
∂ X
t k t k )
(t k ,X t k
(7.164)
This yields

X t k+1 ≈ X t k + a t k , X t k (t k+1 − t k ) + b t k , X t k W t k+1 − W t k

'
∂b t k+1
+ b t k , X t k W t − W t k dW t (7.165)
∂ X
(t k ,X t k ) t k
The last term is the leading-order correction to the explicit Euler method to raise the order
of accuracy to 1. We next evaluate the stochastic integral
'
1
2
t k+1
= dW t = − (t k+1 − t k ) (7.166)
2
I t k ,t k+1 W t − W t k W t k+1 − W t k
t k
To show that (7.166) is valid, we define
2
] − (t − t k ) (7.167)
G(t, Y) = 2I t,t k+1 = [Y − W t k dY = dW t