C.hr weak2.f C....*...1.........2.........3.........4.........5.........6.........7.* C C weak2.f C C Copyright (C) 1995. C C A. Ronald Gallant, P.O. Box 659, Chapel Hill NC 27514-0659, USA C C Permission to use, copy, modify, and distribute this software and C its documentation for any purpose and without fee is hereby C granted, provided that the above copyright notice appear in all C copies and that both that copyright notice and this permission C notice appear in supporting documentation. C C This software is provided "as is" without any expressed or implied C warranty. C C....*...1.........2.........3.........4.........5.........6.........7.* C.hr weak2 C@ *....*...1.........2.........3.........4.........5.........6.........7.* * * weak2 12/21/94 * * purpose * compute a realization of length n from an explicit order 2 weak * scheme for the autonomous Ito stochastic differential equation * * dX = A(X)dt + B(X)dW * * where X has dimension d and W has dimension m. * * usage * call weak2(drift,difuse,t,X,A,B,d,m,t0,X0,dt,iseed,n,Y) * * libf calls * unsk, ran * * arguments * drift - subroutine with usage call drift(t,X,d,A) that is * declared external in the calling program. A must depend * on X only, not on t, because the system is autonomous. * input, external * difuse - subroutine with usage call difuse(t,X,d,m,B) that is * declared external in the calling program. B must depend * on X only, not on t, because the system is autonomous. * input, external * t - input to drift and difuse, on return t=t0+n*dt. * workspace, real*8 * X - input to drift and difuse, a vector of length d. on * return X(i)=Y(i,n) for i=1,...,d. * workspace, real*8 * A - output of drift, a vector of length d. * workspace, real*8 * B - output of difuse, a matrix of order d by m. * workspace, real*8 * d - dimension of X and A, row dimension of B and Y. d must * be less than parameter md which is currently set to 15. * input, integer*4 * m - column dimension of B. m must be less than parameter mm * which is currently set to 5. * input, integer*4 * t0 - time of initial condition. * input, real*8 * X0 - initial condition at time t=t0, a vector of length d. * input, real*8 * dt - time increment for the simulations. * input, real*8 * iseed - seed for the simulations. * input and output, integer*4 * n - number of simulations, column dimension of Y. * input, integer*4 * Y - computed simulations, a d by n matrix. * - output, real*8 * * remarks * 1. To reduce memory requirements for a computation such as * (1/T)(integral from 0 to T of func[X(t)]) where T = n*dt, * dimension as Y(d,n0) and use * do k=1,n,n0 * call weak2(drift,difuse,t0,X0,A,B,d,m,t0,X0,dt,iseed,n0,Y) * do it=1,n0 * do i=1,d * ave(i) = ave(i) + func(Y(i,it))/n * end do * end do * end do * instead of dimensioning as Y(d,n) and using * call weak2(drift,difuse,t,X,A,B,d,m,t0,X0,dt,iseed,n,Y) * do it=1,n * do i=1,d * ave(i) = ave(i) + func(Y(i,it))/n * end do * end do * 2. Let ave[dt,n] denote an average over the time interval [0,T] * where T=n*dt computed as in remark 1. Then * (1/21)*(32*ave[dt,n] - 12*ave[2*dt,n/2] + ave[4*dt,n/4]) * is an order 4 weak estimate (Kloeden and Platen, 1992, p. 492). * * reference * Kloeden, Peter E., and Eckhard Platen (1992), Numerical Solution * of Stochastic Differential Equations, Berlin, Springer-Verlag, * pp. 486-487. * subroutine weak2(drift,difuse,t,X,A,B,d,m,t0,X0,dt,iseed,n,Y) implicit none * arguments external drift, difuse real*8 X, A, B, X0, Y real*8 t, t0, dt integer*4 d, m, iseed, n dimension X(d), A(d), B(d,m), X0(d), Y(d,n) * parameters integer*4 md, mm integer*4 nout parameter (md=15, mm=5) parameter (nout=3) * local variables intrinsic DSQRT real*4 u,ran,z,unsk real*8 dW, V real*8 Tn, ATn real*8 Rp, Rm, Up, Um real*8 BRp, BRm, BUp, BUm real*8 sum real*8 rdt, r3dt integer*4 i,j,k,it dimension dW(mm), V(mm,mm) dimension Tn(md), ATn(md) dimension Rp(md,mm), Rm(md,mm), Up(md,mm), Um(md,mm) dimension BRp(md*mm), BRm(md*mm), BUp(md*mm), BUm(md*mm) dimension sum(md) * error traps if (m.gt.mm) then write(nout,'(a,i5)') 'Error, weak2, mm < m =', m stop end if if (d.gt.md) then write(nout,'(a,i5)') 'Error, weak2, md < d =', d stop end if * explicit order 2 weak scheme do i=1,d X(i)=X0(i) end do t = t0 rdt = DSQRT(dt) r3dt = DSQRT(3.d0*dt) do it=1,n call drift(t,X,d,A) call difuse(t,X,d,m,B) do j=1,m * u = 6.e0*ran(iseed) * if (u .le. 4.e0) then * dW(j) = 0.d0 * else if (u .le. 5.e0) then * dW(j) = r3dt * else * dW(j) = -r3dt * end if z = unsk(iseed) dW(j) = rdt*z V(j,j) = -dt do k=j+1,m u = ran(iseed) if (u .le. 0.5e0) then V(k,j) = dt V(j,k) = -dt else V(k,j) = -dt V(j,k) = dt end if end do end do do i=1,d Tn(i) = X(i) + A(i)*dt end do do j=1,m do i=1,d Tn(i) = Tn(i) + B(i,j)*dW(j) Up(i,j) = X(i) + B(i,j)*rdt Rp(i,j) = Up(i,j) + A(i)*dt Um(i,j) = X(i) - B(i,j)*rdt Rm(i,j) = Um(i,j) + A(i)*dt end do end do do i=1,d sum(i) = 0.d0 end do do j=1,m call difuse(t,Rp(1,j),d,m,BRp) call difuse(t,Rm(1,j),d,m,BRm) do i=1,d sum(i) = sum(i) & + (BRp(-d+d*j+i) + BRm(-d+d*j+i) + 2.d0*B(i,j)) & * dW(j) & + (BRp(-d+d*j+i) - BRm(-d+d*j+i)) & * (dW(j)**2 - dt)/rdt end do do k=1,m if (k.ne.j) then call difuse(t,Up(1,k),d,m,BUp) call difuse(t,Um(1,k),d,m,BUm) do i=1,d sum(i) = sum(i) & + (BUp(-d+d*j+i) + BUm(-d+d*j+i) - 2.d0*B(i,j)) & * dW(j)/rdt & + (BUp(-d+d*j+i) - BUm(-d+d*j+i)) & * (dW(j)*dW(k) + V(k,j))/rdt end do end if end do end do call drift(t,Tn,d,ATn) do i=1,d Y(i,it) = X(i) + 0.5d0*(ATn(i) + A(i))*dt + 0.25d0*sum(i) end do * For a stiff system, the call to subroutine drift and the do loop * above are replaced by an equation solver to solve the equation * Y = X + 0.5d0*(A(it,Y)dt + A)*dt + 0.25d0*sum for Y. See p. 500 * of Kloeden and Platen (1992). To allow use of a solver that * uses first derivatives, subroutine drift will need to be revised * to return first derivatives with respect its argument X. do i=1,d X(i) = Y(i,it) end do t = t + dt end do return end