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
