# R/maths_util.r In dynatopmodel: Implementation of the Dynamic TOPMODEL Hydrological Model

```#-------------------------------------------------------------------------------
#  maths utilities including
#  * numerical solution methods
#  * matrix operations
#-------------------------------------------------------------------------------
# given initial conditions x0 at t=t0, obtain approximate solution using Euler's
# method at t=t1 to ODE dx/dt=func(x,t)
solve.euler <- function(x0, t1, func, t0=0, nstep=10, return.all=FALSE, implicit=FALSE, ...)
{
if(implicit)
{
return(# backwards implicit euler scheme using fixed point iteration
solve.euler.implicit(x0, t1, func, t0, nstep, return.all=return.all, ...)
)
}
nx <- length(x0)
x <- matrix(NA,nrow=nstep, ncol=nx, byrow=TRUE)
x[1,]<- x0
dt <- (t1-t0)/nstep
for(i in 2:nstep)
{
t <- t0 + dt*(i-1)
x[i,] <- x[i-1,] + dt*func(x[i-1,], t=t, ...)
}
if(return.all)
{
return(cbind(seq(t0, t1, length.out=nstep), x))
}
return(x[nstep,])
}

# backwards implicit euler scheme using fixed point iteration
solve.euler.implicit <- function(q0, t1, dqdt, t0=0, nstep=10, niter=30, eps=1e-6, return.all=FALSE, ...)
{
nq <- length(q0)
# qk+1=qk + dt * q'(tk+1, qk+1)
# first approximation
q <- q0
dt <- (t1-t0)/nstep
t <- t0
for(istep in 0:nstep)
{
unsolved <- 1:nq
i <- 1
while(i <= niter & length(unsolved)>0)
{
prev <- q[unsolved]
# qk+1 = qk + h *f(tk+1, qk+1)
#
dq.iter <- dt*dqdt(t=t+dt*(istep+1), q=q, ...)
q[unsolved] <- q[unsolved] + dq.iter[unsolved]
dq <- ifelse(prev==q[unsolved], 0, abs(prev-q[unsolved])/prev)
unsolved <- unsolved[-which(dq < eps)]
i <- i+1
}
if(length(unsolved)>0)
{
warning(paste0("Non-convergence in implicit solution, elements ", unsolved, " of ", nq, ". Try increased number of iterations (currently ", niter, ")."))
}
}

return(q)

}

#------------------------------------------------------------------------------#
identity.matrix <- function(n)
{
mat <- matrix(ncol=n,nrow=n)
mat[] <- 0
mat[col(mat)==row(mat)]<-1
return(mat)
}

# Ron E. VanNimwegen wrote:
#   >
#   > I was looking for a similar operator, because R uses scalar products
# > when raising a matrix to a power with "^".  There might be something
# > more elegant, but this little loop function will do what you need for a
# > matrix "mat" raised to a power "pow":
#   >
#   > mp <- function(mat,pow){
#     > ans <- mat
#     > for ( i in 1:(pow-1)){
#       > ans <- mat%*%ans
#       > }
#     > return(ans)
#     > }
# >
#   This function is extremely inefficient for high powers of pow
# [an unhappy choice of variables, because pow is the power
#  function in many languages]
#
# A better method would keep track of the powers of two,
# and would optimize according to it.
#
# For example, in order to compute mat^42, we just have to
# compute mat^2, mat^4, mat^8, mat^16, mat^32, and
# then mat^40 and mat^42, with "only" 7 multiplications, instead
# of 41.
#
# See the technique in...
# http://en.wikipedia.org/wiki/Exponentiation_by_squaring
matrix.power <- function(mat, n)
{
# test if mat is a square matrix
# treat n < 0 and n = 0 -- this is left as an exercise
# trap non-integer n and return an error
if (n == 1) return(mat)
result <- diag(1, ncol(mat))
while (n > 0) {
if (n %% 2 != 0) {
result <- result %*% mat
n <- n - 1
}
mat <- mat %*% mat
n <- n / 2
}
return(result)
}

# ensure sum of all rows in matrix is unity
normalise.rows <- function(w)
{
if(!is.matrix(w)){warning("matrix input expected to normalise.rows"); return(w)}
if(any(rowSums(w)!=1))
{
cat("Normalising rows...\n")

#     w <-apply(w, MARGIN=1,
#               FUN=function(x)
#               {
#                 sum <- sum(x, na.rm=TRUE)
#                 res<- ifelse(sum==0, rep(0, length(x)), x/sum)
#                 return(res)
#               }
#     )

w <- apply(w, MARGIN=1,
FUN=function(x)
{
sum <- sum(x, na.rm=TRUE)
res<- x/sum
return(res)
})
w[!is.finite(w)]  <- 0
#  w[1:nrow(w),]<-w[1:nrow(w),]/rowSums(w)
# transpose to get in rightg row. col order
w <- t(w)
}
# need to transpose to get in right row-col order
return(w)
}

# graph utilities
# produice an adjacency matrix for the graph (directed=FALSE) or digraph (TRUE)
# i, j element is 1 if vertex i connects to vertex j
# see Foulds (1992) Chapter 6 p.76-78
n=max(conns), # max vertexes considered,defaults to max encountered in connection matrix
directed=TRUE)
{
res <- matrix(ncol=n, nrow=n)
res[]<-0

# conns a list of from -> to vertex numbers
ir <-1
while(ir < nrow(conns))
{
from <- conns[ir,1]
to <- conns[ir,2]
res[from, to]<-1
ir <- ir + 1
}
# if undirected make symmetric...
return(res)
}
```

## Try the dynatopmodel package in your browser

Any scripts or data that you put into this service are public.

dynatopmodel documentation built on Jan. 20, 2018, 9:02 a.m.