inst/doc/examples/Nand.R

In deSolve: Solvers for Initial Value Problems of Differential Equations ('ODE', 'DAE', 'DDE')

#-----------------------------------------------------------------------
#  Note: This file was derived from the FORTRAN code nand.f
#  The file description of the original is:
#  "
#     This file is part of the Test Set for IVP solvers
#     http://www.dm.uniba.it/~testset/
#
#        NAND gate
#        index 0 IDE of dimension 14
#
#     DISCLAIMER: see
#     http://www.dm.uniba.it/~testset/disclaimer.php
#
#     The most recent version of this source file can be found at
#     http://www.dm.uniba.it/~testset/src/problems/nand.f
#
#     This is revision
#     $Id: nand.F,v 1.2 2006/10/02 10:29:14 testset Exp $
#  "
#-----------------------------------------------------------------------

library(deSolve)
      
#-----------------------------------------------------------------------
#
# The network equation describing the nand gate
#             C[Y] * Y' - f[Y,t] = 0
# 
# ---------------------------------------------------------------------

Nand <- function(t,       # time point t
                 Y,       # node potentials at time point t
                 Yprime,
                 pars)    # rate of change of Y
{
  #-----------------------------------------------------------------------
  # Voltage-dependent capacitance matrix C(Y) for the network equation
  #             C(Y) * Y' - f(Y,t) = 0
  #-----------------------------------------------------------------------

  CAP[1, 1]   <- CGS
  CAP[1, 5]   <- -CGS
  CAP[2, 2]   <- CGD
  CAP[2, 5]   <- -CGD
  CAP[3, 3]   <- CBDBS(Y[3]-Y[5])
  CAP[3, 5]   <- -CBDBS(Y[3]-Y[5])
  CAP[4, 4]   <- CBDBS(Y[4]-VDD)
  CAP[5, 1]   <- -CGS
  CAP[5, 2]   <- -CGD
  CAP[5, 3]   <- -CBDBS(Y[3]-Y[5])
  CAP[5, 5]   <- CGS+CGD-CAP[5, 3]+  CBDBS(Y[9]-Y[5])+C9
  CAP[5, 9]   <- -CBDBS(Y[9]-Y[5])
  CAP[6, 6]   <- CGS
  CAP[7, 7]   <- CGD
  CAP[8, 8]   <- CBDBS(Y[8]-Y[10])
  CAP[8, 10]  <- -CBDBS(Y[8]-Y[10])
  CAP[9, 5]   <- -CBDBS(Y[9]-Y[5])
  CAP[9, 9]   <- CBDBS(Y[9]-Y[5])
  CAP[10, 8]  <- -CBDBS(Y[8]-Y[10])
  CAP[10, 10] <- -CAP[8, 10]+CBDBS(Y[14]-Y[10])+C9
  CAP[10, 14] <- -CBDBS(Y[14]-Y[10])
  CAP[11, 11] <- CGS
  CAP[12, 12] <- CGD
  CAP[13, 13] <- CBDBS(Y[13])
  CAP[14, 10] <- -CBDBS(Y[14]-Y[10])
  CAP[14, 14] <- CBDBS(Y[14]-Y[10])

  # ---------------------------------------------------------------------
  #          PULSE: Input signal in pulse form
  # ---------------------------------------------------------------------
  P1  <- PULSE(t, 0.0, 5.0, 5.0, 5.0, 5.0, 5.0, 20.0)
  V1  <- P1$VIN
  V1D <- P1$VIND

  P2  <- PULSE(t, 0.0, 5.0, 15.0, 5.0, 15.0, 5.0, 40.0)
  V2  <- P2$VIN
  V2D <- P2$VIND

  #-----------------------------------------------------------------------
  # Right-hand side f[X,t] for the network equation
  #             C[Y] * Y' - f[Y,t] = 0
  # External reference:
  #          IDS: Drain-source current
  #          IBS: Nonlinear current characteristic for diode between
  #               bulk and source
  #          IBD: Nonlinear current characteristic for diode between
  #               bulk and drain
  #-----------------------------------------------------------------------

  F[1] <- -(Y[1]-Y[5])/RGS-IDS(1, Y[2]-Y[1], Y[5]-Y[1], Y[3]-Y[5],
            Y[5]-Y[2], Y[4]-VDD)
  F[2] <- -(Y[2]-VDD)/RGD+IDS(1, Y[2]-Y[1], Y[5]-Y[1], Y[3]-Y[5],
            Y[5]-Y[2], Y[4]-VDD)
  F[3] <- -(Y[3]-VBB)/RBS + IBS(Y[3]-Y[5])
  F[4] <- -(Y[4]-VBB)/RBD + IBD(Y[4]-VDD)
  F[5] <- -(Y[5]-Y[1])/RGS-IBS(Y[3]-Y[5])-(Y[5]-Y[7])/RGD-
            IBD(Y[9]-Y[5])
  F[6] <- CGS*V1D-(Y[6]-Y[10])/RGS -
            IDS(2, Y[7]-Y[6], V1-Y[6], Y[8]-Y[10], V1-Y[7], Y[9]-Y[5])
  F[7] <- CGD*V1D-(Y[7]-Y[5])/RGD +
            IDS(2, Y[7]-Y[6], V1-Y[6], Y[8]-Y[10], V1-Y[7], Y[9]-Y[5])
  F[8] <- -(Y[8]-VBB)/RBS + IBS(Y[8]-Y[10])
  F[9] <- -(Y[9]-VBB)/RBD + IBD(Y[9]-Y[5])
  F[10] <- -(Y[10]-Y[6])/RGS-IBS(Y[8]-Y[10]) -
            (Y[10]-Y[12])/RGD-IBD(Y[14]-Y[10])
  F[11] <- CGS*V2D-Y[11]/RGS-IDS(2, Y[12]-Y[11], V2-Y[11], Y[13],
             V2-Y[12], Y[14]-Y[10])
  F[12] <- CGD*V2D-(Y[12]-Y[10])/RGD +
             IDS(2, Y[12]-Y[11], V2-Y[11], Y[13], V2-Y[12], Y[14]-Y[10])
  F[13] <- -(Y[13]-VBB)/RBS + IBS(Y[13])
  F[14] <- -(Y[14]-VBB)/RBD + IBD(Y[14]-Y[10])

  #             C[Y] * Y' - f[Y,t] = 0
  Delta <- colSums(t(CAP)*Yprime)-F
  return(list(c(Delta), pulse1 = P1$VIN, pulse2 = P2$VIN))
}

# ---------------------------------------------------------------------------
#
# Function evaluating the drain-current due to the model of
# Shichman and Hodges
#
# ---------------------------------------------------------------------------

IDS <- function (NED,   #   NED  Integer parameter for MOSFET-type
                 VDS,   #   VDS  Voltage between drain and source
                 VGS,   #   VGS  Voltage between gate and source
                 VBS,   #   VBS  Voltage between bulk and source
                 VGD,   #   VGD  Voltage between gate and drain
                 VBD)   #   VBD  Voltage between bulk and drain

{
  if ( VDS == 0 ) return(0)
   
  if (NED== 1) { #--- Depletion-type
    VT0     <- -2.43
    CGAMMA  <-  0.2
    PHI     <-  1.28
    BETA    <-  5.35e-4
  } else       { # --- Enhancement-type
    VT0     <-  0.2
    CGAMMA  <-  0.035
    PHI     <-  1.01
    BETA    <-  1.748e-3
  }

  if ( VDS >  0 )  # drain function for VDS>0
  {
    SQRT1<-ifelse (PHI-VBS>0, sqrt(PHI-VBS), 0)
    VTE  <-  VT0 + CGAMMA * ( SQRT1 - sqrt(PHI) )

    if ( VGS-VTE <= 0.0) IDS  <-  0.  else
    if ( 0.0 < VGS-VTE & VGS-VTE <= VDS ) 
     IDS  <-  - BETA * (VGS - VTE)^ 2.0 * (1.0 + DELTA*VDS) else
    if ( 0.0 < VDS & VDS < VGS-VTE ) 
     IDS  <-  - BETA * VDS * (2 *(VGS - VTE) - VDS) * (1.0 + DELTA*VDS)

  }  else    {

    SQRT2<-ifelse (PHI-VBD>0, sqrt(PHI-VBD), 0)
    VTE  <-  VT0 + CGAMMA * (SQRT2 - sqrt(PHI) )

    if ( VGD-VTE <= 0.0) IDS  <-  0.0  else
    if ( 0.0 < VGD-VTE & VGD-VTE <= -VDS ) 
     IDS  <-  BETA * (VGD - VTE)^2.0 * (1.0 - DELTA*VDS) else
    if ( 0.0 < -VDS & -VDS < VGD-VTE ) 
     IDS  <-  - BETA * VDS * (2 *(VGD - VTE) + VDS) *(1.0 - DELTA*VDS)
  }
  return(IDS)
}

# ---------------------------------------------------------------------------
#
# Function evaluating the current of the pn-junction between bulk and
# source due to the model of Shichman and Hodges
#
# ---------------------------------------------------------------------------

IBS <- function(VBS)     #   VBS  Voltage between bulk and source
         ifelse (VBS <= 0.0, -CURIS * (exp(VBS/VTH) - 1.0), 0.0)

# ---------------------------------------------------------------------------
#
# Function evaluating the current of the pn-junction between bulk and
# drain  due to the model of Shichman and Hodges
#
# ---------------------------------------------------------------------------
IBD <- function(VBD)     #   VBD  Voltage between bulk and drain
         ifelse(VBD <= 0.0, -CURIS * (exp(VBD/VTH) - 1.0), 0.0)

# ---------------------------------------------------------------------------
#
# Evaluating input signal at time point X
#
# ---------------------------------------------------------------------------

PULSE <- function (X,     # Time-point at which input signal is evaluated
                   LOW,   # Low-level of input signal
                   HIGH,  # High-level of input signal
                   DELAY, T1, T2, T3, PERIOD)  # Parameters to specify signal structure
# ---------------------------------------------------------------------------
# Structure of input signal:
#
#                -----------------------                       HIGH
#               /                       \
#              /                         \
#             /                           \
#            /                             \
#           /                               \
#          /                                 \
#         /                                   \
#        /                                     \
#  ------                                       ---------      LOW
#
# |DELAY|   T1  |         T2           |   T3  |
# |          P     E     R     I     O     D            |
#
# ---------------------------------------------------------------------------
{
  TIME  <-  X %% PERIOD
  VIN   <-  LOW
  VIND  <-  0.0

  if (TIME > (DELAY+T1+T2)) {
    VIN   <-  ((HIGH-LOW)/T3)*(DELAY+T1+T2+T3-TIME) + LOW
    VIND  <-  -((HIGH-LOW)/T3)
  } else if (TIME > (DELAY+T1)) {
    VIN  <-  HIGH
    VIND <-  0.0
  } else if (TIME > DELAY) {
    VIN  <-  ((HIGH-LOW)/T1)*(TIME-DELAY) + LOW
    VIND <-  ((HIGH-LOW)/T1)
  }

  return  (list(VIN = VIN,  # Voltage of input signal at time point X
                VIND = VIND))  # Derivative of VIN at time point X

}

# ---------------------------------------------------------------------------
#
# Function evaluating the voltage-dependent capacitance between bulk and
# drain gevalp. source  due to the model of Shichman and Hodges
#
# ---------------------------------------------------------------------------

CBDBS <- function(V)   # Voltage between bulk and drain gevalp. source
           ifelse(V <= 0.0, CBD/sqrt(1.0-V/0.87), CBD*(1.0+V/(2.0*0.87)))


#-----------------------------------------------------------------------
# solution
# computed at Cray C90, using Cray double precision:
# Solving NAND gate using PSIDE
#
# User input:
#
# give relative error tolerance: 1d-16
# give absolute error tolerance: 1d-16
#
#
# Integration characteristics:
#
#    number of integration steps       22083
#    number of accepted steps          21506
#    number of f evaluations          308562
#    number of Jacobian evaluations      337
#    number of LU decompositions       10532
#
# CPU-time used:                         451.71 sec
#
#      y[  1] =  0.4971088699385777d+1
#      y[  2] =  0.4999752103929311d+1
#      y[  3] = -0.2499998781491227d+1
#      y[  4] = -0.2499999999999975d+1
#      y[  5] =  0.4970837023296724d+1
#      y[  6] = -0.2091214032073855d+0
#      y[  7] =  0.4970593243278363d+1
#      y[  8] = -0.2500077409198803d+1
#      y[  9] = -0.2499998781491227d+1
#      y[ 10] = -0.2090289583878100d+0
#      y[ 11] = -0.2399999999966269d-3
#      y[ 12] = -0.2091214032073855d+0
#      y[ 13] = -0.2499999999999991d+1
#      y[ 14] = -0.2500077409198803d+1
#-----------------------------------------------------------------------

RGS     <-  4
RGD     <-  4
RBS     <-  10
RBD     <-  10
CGS     <-  0.6e-4
CGD     <-  0.6e-4
CBD     <-  2.4e-5
CBS     <-  2.4e-5
C9      <-  0.5e-4
DELTA   <-  0.2e-1
CURIS   <-  1.e-14
VTH     <-  25.85
VDD     <-  5.
VBB     <-  -2.5

#-----------------------------------------------------------------------
# initialising
VBB    <-  -2.5
Y      <- c(5, 5, VBB, VBB, 5, 3.62385, 5, VBB, VBB, 3.62385, 0, 3.62385, VBB, VBB)
Yprime <- rep(0, 14)

#-----------------------------------------------------------------------
# memory allocation
CAP   <- matrix(nrow = 14, ncol = 14, 0)
F     <- vector("double", 14)
times <- seq(0, 80, by = 1)         # time: from 0 to 80 hours, steps of 1 hour

# integrate the model: low tolerances to restrict integration time
out <- daspk(y = Y, dy = NULL, times, res = Nand, parms = 0,
             rtol = 1e-6, atol = 1e-6)

# plot output
par(mfrow = c(4, 4), mar = c(4, 2, 3, 2))
for(i in 2:15) plot(out[, 1], out[, i], type = "l", ylab = "",
                    main = paste("y[", i-1, "]"), xlab = "time")

# reference solution
ref<-c(4.971088699385777,   4.999752103929311,     -2.499998781491227,
      -2.499999999999975,   4.970837023296724,     -0.2091214032073855,
       4.970593243278363,  -2.500077409198803,     -2.499998781491227,
      -0.2090289583878100, -0.2399999999966269e-3, -0.2091214032073855,
      -2.499999999999991,  -2.500077409198803)

t(rbind(daspk = out [nrow(out), 2:15] ,
    reference = ref,
         delt = out [nrow(out), 2:15] - ref)
)