Practical Numerical Math Routines

Description

This package provides R implementations of more advanced functions in numerical analysis, with a special view on on optimization and time series routines. Uses Matlab/Octave function names where appropriate to simplify porting.

Some of these implementations are the result of courses on Scientific Computing (“Wissenschaftliches Rechnen”) and are mostly intended to demonstrate how to implement certain algorithms in R/S. Others are implementations of algorithms found in textbooks.

Details

The package encompasses functions from all areas of numerical analysis, for example:

• Root finding and minimization of univariate functions,
  e.g. Newton-Raphson, Brent-Dekker, Fibonacci or ‘golden ratio’ search.

• Handling polynomials, including roots and polynomial fitting,
  e.g. Laguerre's and Muller's methods.

• Interpolation and function approximation,
  barycentric Lagrange interpolation, Pade and rational interpolation, Chebyshev or trigonometric approximation.

• Some special functions,
  e.g. Fresnel integrals, Riemann's Zeta or the complex Gamma function, and Lambert's W computed iteratively through Newton's method.

• Special matrices, e.g. Hankel, Rosser, Wilkinson

• Numerical differentiation and integration,
  Richardson approach and “complex step” derivatives, adaptive Simpson and Lobatto integration and adaptive Gauss-Kronrod quadrature.

• Solvers for ordinary differential equations and systems,
  Euler-Heun, classical Runge-Kutta, ode23, or predictor-corrector method such as the Adams-Bashford-Moulton.

• Some functions from number theory,
  such as primes and prime factorization, extended Euclidean algorithm.

• Sorting routines, e.g. recursive quickstep.

• Several functions for string manipulation and regular search, all wrapped and named similar to their Matlab analogues.

It serves three main goals:

• Collecting R scripts that can be demonstrated in courses on ‘Numerical Analysis’ or ‘Scientific Computing’ using R/S as the chosen programming language.

• Wrapping functions with appropriate Matlab names to simplify porting programs from Matlab or Octave to R.

• Providing an environment in which R can be used as a full-blown numerical computing system.

Besides that, many of these functions could be called in R applications as they do not have comparable counterparts in other R packages (at least at this moment, as far as I know).

All referenced books have been utilized in one way or another. Web links have been provided where reasonable.

Note

The following 220 functions are emulations of correspondingly named Matlab functions and bear the same signature as their Matlab cousins if possible:

accumarray, acosd, acot, acotd, acoth, acsc, acscd, acsch, and, angle, ans,
arrayfun, asec, asecd, asech, asind, atand, atan2d,
beep, bernoulli, blank, blkdiag, bsxfun,
cart2pol, cart2sph, cd, ceil, circshift, clear, compan, cond, conv,
cosd, cot, cotd, coth, cross, csc, cscd, csch, cumtrapz,
dblquad, deblank, deconv, deg2rad, detrend, deval, disp, dot,
eig, eigint, ellipj, ellipke, eps, erf, erfc, erfcinv, erfcx, erfi, erfinv,
errorbar, expint, expm, eye, ezcontour, ezmesh, ezplot, ezpolar, ezsurf,
fact, fftshift, figure, findpeaks, findstr, flipdim, fliplr, flipud,
fminbnd, fmincon, fminsearch, fminunc, fplot, fprintf, fsolve, fzero,
gammainc, gcd, geomean, gmres, gradient,
hadamard, hankel, harmmean, hilb, histc, humps, hypot,
idivide, ifft, ifftshift, inpolygon, integral, integral2, integral3,
interp1, interp2, inv, isempty, isprime,
kron,
legendre, linprog, linspace, loglog, logm, logseq, logspace, lsqcurvefit,
lsqlin, lsqnonlin, lsqnonneg, lu,
magic, meshgrid, mkpp, mldivide, mod, mrdivide,
nchoosek, ndims, nextpow2, nnz, normest, nthroot, null, num2str, numel,
ode23, ode23s, ones, or, orth,
pascal, pchip, pdist, pdist2, peaks, perms, piecewise, pinv, plotyy,
pol2cart, polar, polyfit, polyint, polylog, polyval, pow2, ppval,
primes, psi, pwd,
quad, quad2d, quadgk, quadl, quadprog, quadv, quiver,
rad2deg, randi, randn, randsample, rat, rats, regexp, regexpi,
regexpreg, rem, repmat, roots, rosser, rot90, rref, runge,
sec, secd, sech, semilogx, semilogy, sinc, sind, size, sortrows, sph2cart,
sqrtm, squareform, std, str2num, strcat, strcmp, strcmpi,
strfind, strfindi, strjust, subspace,
tand, tic, toc, trapz, tril, trimmean, triplequad, triu,
vander, vectorfield, ver,
what, who, whos, wilkinson,
zeros, zeta

The following Matlab function names have been capitalized in ‘pracma’ to avoid shadowing functions from R base or one of its recommended packages (on request of Bill Venables and because of Brian Ripley's CRAN policies):

Diag, factos, finds, Fix, Imag, Lcm, Mode, Norm, nullspace (<- null),
Poly, Rank, Real, Reshape, strRep, strTrim, Toeplitz, Trace, uniq (<- unique).

To use “ans” instead of “ans()” – as is common practice in Matlab – type (and similar for other Matlab commands):

makeActiveBinding("ans", function() .Last.value, .GlobalEnv)
makeActiveBinding("who", who(), .GlobalEnv)

Author(s)

Hans Werner Borchers

Maintainer: Hans W Borchers <hwborchers@googlemail.com>

References

Abramowitz, M., and I. A. Stegun (1972). Handbook of Mathematical Functions (with Formulas, Graphs, and Mathematical Tables). Dover, New York. URL: https://www.math.ubc.ca/~cbm/aands/notes.htm

Arndt, J. (2010). Matters Computational: Ideas, Algorithms, Source Code. Springer-Verlag, Berlin Heidelberg Dordrecht. FXT: a library of algorithms: https://www.jjj.de/fxt/.

Cormen, Th. H., Ch. E. Leiserson, and R. L. Rivest (2009). Introduction to Algorithms. Third Edition, The MIT Press, Cambridge, MA.

Encyclopedia of Mathematics (2012). Editor-in-Chief: Ulf Rehmann. https://encyclopediaofmath.org/wiki/Main_Page.

Gautschi, W. (1997). Numerical Analysis: An Introduction. Birkhaeuser, Boston.

Gentle, J. E. (2009). Computational Statistics. Springer Science+Business Media LCC, New York.

MathWorld.com (2011). Matlab Central: https://www.mathworks.com/matlabcentral/.

NIST: National Institute of Standards and Technology. Olver, F. W. J., et al. (2010). NIST Handbook of Mathematical Functions. Cambridge University Press. Internet: NIST Digital Library of Mathematical Functions, https://dlmf.nist.gov/; Guide to Available Mathematical Software, https://gams.nist.gov/.

Press, W. H., S. A. Teukolsky, W. T Vetterling, and B. P. Flannery (2007). Numerical Recipes: The Art of Numerical Computing. Third Edition, incl. Numerical Recipes Software, Cambridge University Press, New York. URL: numerical.recipes/book/book.html.

Quarteroni, A., R. Sacco, and F. Saleri (2007). Numerical Mathematics. Second Edition, Springer-Verlag, Berlin Heidelberg.

Skiena, St. S. (2008). The Algorithm Design Manual. Second Edition, Springer-Verlag, London. The Stony Brook Algorithm Repository: https://algorist.com/algorist.html.

Stoer, J., and R. Bulirsch (2002). Introduction to Numerical Analysis. Third Edition, Springer-Verlag, New York.

Strang, G. (2007). Computational Science and Engineering. Wellesley-Cambridge Press.

Weisstein, E. W. (2003). CRC Concise Encyclopedia of Mathematics. Second Edition, Chapman & Hall/CRC Press. Wolfram MathWorld: https://mathworld.wolfram.com/.

Zhang, S., and J. Jin (1996). Computation of Special Functions. John Wiley & Sons.

See Also

The R package ‘matlab’ contains some of the basic routines from Matlab, but unfortunately not any of the higher math routines.

Examples

## Not run: 
##  See examples in the help files for all functions.
    
## End(Not run)

Adams-Bashford-Moulton

Description

Third-order Adams-Bashford-Moulton predictor-corrector method.

Usage

abm3pc(f, a, b, y0, n = 50, ...)

Arguments

Details

Combined Adams-Bashford and Adams-Moulton (or: multi-step) method of third order with corrections according to the predictor-corrector approach.

Value

List with components x for grid points between a and b and y a vector y the same length as x; additionally an error estimation est.error that should be looked at with caution.

Note

This function serves demonstration purposes only.

References

Fausett, L. V. (2007). Applied Numerical Analysis Using Matlab. Second edition, Prentice Hall.

See Also

rk4, ode23

Examples

##  Attempt on a non-stiff equation
#   y' = y^2 - y^3, y(0) = d, 0 <= t <= 2/d, d = 0.01
f <- function(t, y) y^2 - y^3
d <- 1/250
abm1 <- abm3pc(f, 0, 2/d, d, n = 1/d)
abm2 <- abm3pc(f, 0, 2/d, d, n = 2/d)
## Not run: 
plot(abm1$x, abm1$y, type = "l", col = "blue")
lines(abm2$x, abm2$y, type = "l", col = "red")
grid()
## End(Not run)

Accumulate Vector Elements

Description

accumarray groups elements from a data set and applies a function to each group.

Usage

accumarray(subs, val, sz = NULL, func = sum, fillval = 0)

uniq(a, first = FALSE)

Arguments

Details

A <- accumarray(subs, val) creates an array A by accumulating elements of the vector val using the lines of subs as indices and applying func to that accumulated vector. The size of the array can be predetermined by the size vector sz.

A = uniq(a) returns a vector b identical to unique(a) and two other vectors of indices m and n such that b == a[m] and a == b[n].

Value

accumarray returns an array of size the maximum in each column of subs, or by sz.

uniq returns a list with components

Note

The Matlab function accumarray can also handle sparse matrices.

See Also

unique

Examples

##  Examples for accumarray
val = 101:105
subs = as.matrix(c(1, 2, 4, 2, 4))
accumarray(subs, val)
# [101; 206; 0; 208]

val = 101:105
subs <- matrix(c(1,2,2,2,2, 1,1,3,1,3, 1,2,2,2,2), ncol = 3)
accumarray(subs, val)
# , , 1
# [,1] [,2] [,3]
# [1,]  101    0    0
# [2,]    0    0    0
# , , 2
# [,1] [,2] [,3]
# [1,]    0    0    0
# [2,]  206    0  208

val = 101:106
subs <- matrix(c(1, 2, 1, 2, 3, 1, 4, 1, 4, 4, 4, 1), ncol = 2, byrow = TRUE)
accumarray(subs, val, func = function(x) sum(diff(x)))
# [,1] [,2] [,3] [,4]
# [1,]    0    1    0    0
# [2,]    0    0    0    0
# [3,]    0    0    0    0
# [4,]    2    0    0    0

val = 101:105
subs = matrix(c(1, 1, 2, 1, 2, 3, 2, 1, 2, 3), ncol = 2, byrow = TRUE)
accumarray(subs, val, sz = c(3, 3), func = max, fillval = NA)
# [,1] [,2] [,3]
# [1,]  101   NA   NA
# [2,]  104   NA  105
# [3,]   NA   NA   NA

##  Examples for uniq
a <- c(1, 1, 5, 6, 2, 3, 3, 9, 8, 6, 2, 4)
A <- uniq(a); A
# A$b  1  5  6  2  3  9  8  4
# A$m  2  3 10 11  7  8  9 12
# A$n  1  1  2  3  4  5  5  6  7  3  4  8
A <- uniq(a, first = TRUE); A
# A$m  1  3  4  5  6  8  9 12

##  Example: Subset sum problem
# Distribution of unique sums among all combinations of a vectors.
allsums <- function(a) {
    S <- c(); C <- c()
    for (k in 1:length(a)) {
        U <- uniq(c(S, a[k], S + a[k]))
        S <- U$b
        C <- accumarray(U$n, c(C, 1, C))
    }
    o <- order(S); S <- S[o]; C <- C[o]
    return(list(S = S, C = C))
}
A <- allsums(seq(1, 9, by=2)); A
# A$S  1  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 24 25
# A$C  1  1  1  1  1  1  2  2  2  1  2  2  1  2  2  2  1  1  1  1  1  1  1

Arithmetic-geometric Mean

Description

The arithmetic-geometric mean of real or complex numbers.

Usage

agmean(a, b)

Arguments

Details

The arithmetic-geometric mean is defined as the common limit of the two sequences a_{n+1} = (a_n + b_n)/2 and b_{n+1} = \sqrt(a_n b_n).

When used for negative or complex numbers, the complex square root function is applied.

Value

Returns a list with compoinents: agm a vector of arithmetic-geometric means, component-wise, niter the number of iterations, and prec the overall estimated precision.

Note

Gauss discovered that elliptic integrals can be effectively computed via the arithmetic-geometric mean (see example below), for example:

\int_0^{\pi/2} \frac{dt}{\sqrt{1 - m^2 sin^2(t)}} = \frac{(a+b) \pi}{4 \cdot agm(a,b)}

where m = (a-b)/(a+b)

References

https://mathworld.wolfram.com/Arithmetic-GeometricMean.html

See Also

Arithmetic, geometric, and harmonic mean.

Examples

##  Accuracy test: Gauss constant
1/agmean(1, sqrt(2))$agm - 0.834626841674073186  # 1.11e-16 < eps = 2.22e-16

## Gauss' AGM-based computation of \pi
a <- 1.0
b <- 1.0/sqrt(2)
s <- 0.5
d <- 1L
while (abs(a-b) > eps()) {
    t <- a
    a <- (a + b)*0.5
    b <- sqrt(t*b)
    c <- (a-t)*(a-t)
    d <- 2L * d
    s <- s - d*c
}
approx_pi <- (a+b)^2 / s / 2.0
abs(approx_pi - pi)             # 8.881784e-16 in 4 iterations

##  Example: Approximate elliptic integral
N <- 20
m <- seq(0, 1, len = N+1)[1:N]
E <- numeric(N)
for (i in 1:N) {
    f <- function(t) 1/sqrt(1 - m[i]^2 * sin(t)^2)
    E[i] <- quad(f, 0, pi/2)
}
A <- numeric(2*N-1)
a <- 1
b <- a * (1-m) / (m+1)

## Not run: 
plot(m, E, main = "Elliptic Integrals vs. arith.-geom. Mean")
lines(m, (a+b)*pi / 4 / agmean(a, b)$agm, col="blue")
grid()
## End(Not run)

Aitken' Method

Description

Aitken's acceleration method.

Usage

aitken(f, x0, nmax = 12, tol = 1e-8, ...)

Arguments

Details

Aitken's acceleration method, or delta-squared process, is used for accelerating the rate of convergence of a sequence (from linear to quadratic), here applied to the fixed point iteration scheme of a function.

Value

The fixpoint (as found so far).

Note

Sometimes used to accerate Newton-Raphson (Steffensen's method).

References

Quarteroni, A., and F. Saleri (2006). Scientific Computing with Matlab and Octave. Second Edition, Springer-Verlag, Berlin Heidelberg.

See Also

lambertWp

Examples

# Find a zero of    f(x) = cos(x) - x*exp(x)
# as fixpoint of  phi(x) = x + (cos(x) - x*exp(x))/2
phi <- function(x) x + (cos(x) - x*exp(x))/2
aitken(phi, 0)  #=> 0.5177574

Univariate Akima Interpolation

Description

Interpolate smooth curve through given points on a plane.

Usage

  akimaInterp(x, y, xi)

Arguments

Details

Implementation of Akima's univariate interpolation method, built from piecewise third order polynomials. There is no need to solve large systems of equations, and the method is therefore computationally very efficient.

Value

Returns the interpolated values at the points xi as a vector.

Note

There is also a 2-dimensional version in package ‘akima’.

Author(s)

Matlab code by H. Shamsundar under BSC License; re-implementation in R by Hans W Borchers.

References

Akima, H. (1970). A New Method of Interpolation and Smooth Curve Fitting Based on Local Procedures. Journal of the ACM, Vol. 17(4), pp 589-602.

Hyman, J. (1983). Accurate Monotonicity Preserving Cubic Interpolation. SIAM J. Sci. Stat. Comput., Vol. 4(4), pp. 645-654.

Akima, H. (1996). Algorithm 760: Rectangular-Grid-Data Surface Fitting that Has the Accurancy of a Bicubic Polynomial. ACM TOMS Vol. 22(3), pp. 357-361.

Akima, H. (1996). Algorithm 761: Scattered-Data Surface Fitting that Has the Accuracy of a Cubic Polynomial. ACM TOMS, Vol. 22(3), pp. 362-371.

See Also

kriging, akima::aspline, akima::interp

Examples

x <- c( 0,  2,  3,  5,  6,  8,  9,   11, 12, 14, 15)
y <- c(10, 10, 10, 10, 10, 10, 10.5, 15, 50, 60, 85)
xs <- seq(12, 14, 0.5)          # 12.0 12.5     13.0     13.5     14.0
ys <- akimaInterp(x, y, xs)     # 50.0 54.57405 54.84360 55.19135 60.0
xs; ys

## Not run: 
plot(x, y, col="blue", main = "Akima Interpolation")
xi <- linspace(0,15,51)
yi <- akimaInterp(x, y, xi)
lines(xi, yi, col = "darkred")
grid()
## End(Not run)

Logical AND, OR (Matlab Style)

Description

and(l, k) resp. or(l, k) the same as (l & k) + 0 resp. (l | k) + 0.

Usage

and(l, k)
or(l, k)

Arguments

Details

Performs a logical operation of arrays l and k and returns an array containing elements set to either 1 (TRUE) or 0 (FALSE), that is in Matlab style.

Value

Logical vector.

Examples

A <- matrix(c(0.5,  0.5,  0,    0.75, 0,
              0.5,  0,    0.75, 0.05, 0.85,
              0.35, 0,    0,    0,    0.01,
              0.5,  0.65, 0.65, 0.05, 0), 4, 5, byrow=TRUE)
B <- matrix(c( 0, 1, 0, 1, 0,
               1, 1, 1, 0, 1,
               0, 1, 1, 1, 0,
               0, 1, 0, 0, 1), 4, 5, byrow=TRUE)

and(A, B)
or(A, B)

Andrews' Curves

Description

Plots Andrews' curves in cartesian or polar coordinates.

Usage

andrewsplot(A, f, style = "pol", scaled = FALSE, npts = 101)

Arguments

Details

andrewsplot creates an Andrews plot of the multivariate data in the matrix A, assigning different colors according to the factor or integer vector f.

Andrews' plot represent each observation (row) by a periodic function over the interval [0, 2*pi]. This function for the i-th observation is defined as ...

The plot can be seen in cartesian or polar coordinates — the latter seems appropriate as all these functions are periodic.

Value

Generates a plot, no return value.

Note

Please note that a different ordering of the columns will result in quite different functions and overall picture.

There are variants utilizing principal component scores, in order of decreasing eigenvalues.

References

R. Khattree and D. N. Naik (2002). Andrews PLots for Multivariate Data: Some New Suggestions and Applications. Journal of Statistical Planning and Inference, Vol. 100, No. 2, pp. 411–425.

See Also

polar, andrews::andrews

Examples

## Not run: 
data(iris)
s <- sample(1:4, 4)
A <- as.matrix(iris[, s])
f <- as.integer(iris[, 5])
andrewsplot(A, f, style = "pol")

## End(Not run)

Basic Complex Functions

Description

Basic complex functions (Matlab style)

Usage

Real(z)
Imag(z)
angle(z)

Arguments

Details

These are just Matlab names for the corresponding functions in R. The angle function is simply defined as atan2(Im(z), Re(z)).

Value

returning real or complex values; angle returns in radians.

Note

The true Matlab names are real, imag, and conj, but as real was taken in R, all these beginnings are changed to capitals.

The function Mod has no special name in Matlab; use abs() instead.

See Also

Mod, abs

Examples

z <- c(0, 1, 1+1i, 1i)
Real(z)   # Re(z)
Imag(z)   # Im(z)
Conj(z)   # Conj(z)
abs(z)    # Mod(z)
angle(z)

Adaptive Nelder-Mead Minimization

Description

An implementation of the Nelder-Mead algorithm for derivative-free optimization / function minimization.

Usage

anms(fn, x0, ...,
     tol = 1e-10, maxfeval = NULL)

Arguments

Details

Also called a ‘simplex’ method for finding the local minimum of a function of several variables. The method is a pattern search that compares function values at the vertices of the simplex. The process generates a sequence of simplices with ever reducing sizes.

anms can be used up to 20 or 30 dimensions (then ‘tol’ and ‘maxfeval’ need to be increased). It applies adaptive parameters for simplicial search, depending on the problem dimension – see Fuchang and Lixing (2012).

With upper and/or lower bounds, anms will apply a transformation of bounded to unbounded regions before utilizing Nelder-Mead. Of course, if the optimum is near to the boundary, results will not be as accurate as when the minimum is in the interior.

Value

List with following components:

Note

Copyright (c) 2012 by F. Gao and L. Han, implemented in Matlab with a permissive license. Implemented in R by Hans W. Borchers. For another elaborate implementation of Nelder-Mead see the package ‘dfoptim’.

References

Nelder, J., and R. Mead (1965). A simplex method for function minimization. Computer Journal, Volume 7, pp. 308-313.

O'Neill, R. (1971). Algorithm AS 47: Function Minimization Using a Simplex Procedure. Applied Statistics, Volume 20(3), pp. 338-345.

J. C. Lagarias et al. (1998). Convergence properties of the Nelder-Mead simplex method in low dimensions. SIAM Journal for Optimization, Vol. 9, No. 1, pp 112-147.

Fuchang Gao and Lixing Han (2012). Implementing the Nelder-Mead simplex algorithm with adaptive parameters. Computational Optimization and Applications, Vol. 51, No. 1, pp. 259-277.

See Also

optim

Examples

##  Rosenbrock function
rosenbrock <- function(x) {
    n <- length(x)
    x1 <- x[2:n]
    x2 <- x[1:(n-1)]
    sum(100*(x1-x2^2)^2 + (1-x2)^2)
}

anms(rosenbrock, c(0,0,0,0,0))
# $xmin
# [1] 1 1 1 1 1
# $fmin
# [1] 8.268732e-21
# $nfeval
# [1] 1153

# To add constraints to the optimization problem, use a slightly 
# modified objective function. Equality constraints not possible.
# Warning: Avoid a starting value too near to the boundary !

## Not run: 
# Example: 0.0 <= x <= 0.5
fun <- function(x) {
    if (any(x < 0) || any(x > 0.5)) 100
    else rosenbrock(x)
}
x0 <- rep(0.1, 5)

anms(fun, x0)
## $xmin
## [1] 0.500000000 0.263051265 0.079972922 0.016228138 0.000267922
## End(Not run)

Approximate and Sample Entropy

Description

Calculates the approximate or sample entropy of a time series.

Usage

approx_entropy(ts, edim = 2, r = 0.2*sd(ts), elag = 1)

sample_entropy(ts, edim = 2, r = 0.2*sd(ts), tau = 1)

Arguments

Details

Approximate entropy was introduced to quantify the the amount of regularity and the unpredictability of fluctuations in a time series. A low value of the entropy indicates that the time series is deterministic; a high value indicates randomness.

Sample entropy is conceptually similar with the following differences: It does not count self-matching, and it does not depend that much on the length of the time series.

Value

The approximate, or sample, entropy, a scalar value.

Note

This code here derives from Matlab versions at Mathwork's File Exchange, “Fast Approximate Entropy” and “Sample Entropy” by Kijoon Lee under BSD license.

References

Pincus, S.M. (1991). Approximate entropy as a measure of system complexity. Proc. Natl. Acad. Sci. USA, Vol. 88, pp. 2297–2301.

Kaplan, D., M. I. Furman, S. M. Pincus, S. M. Ryan, L. A. Lipsitz, and A. L. Goldberger (1991). Aging and the complexity of cardiovascular dynamics, Biophysics Journal, Vol. 59, pp. 945–949.

Yentes, J.M., N. Hunt, K.K. Schmid, J.P. Kaipust, D. McGrath, N. Stergiou (2012). The Appropriate use of approximate entropy and sample entropy with short data sets. Ann. Biomed. Eng.

See Also

RHRV::CalculateApEn

Examples

ts <- rep(61:65, 10)
approx_entropy(ts, edim = 2)                      # -0.0004610253
sample_entropy(ts, edim = 2)                      #  0

set.seed(8237)
approx_entropy(rnorm(500), edim = 2)              # 1.351439  high, random
approx_entropy(sin(seq(1,100,by=0.2)), edim = 2)  # 0.171806  low,  deterministic
sample_entropy(sin(seq(1,100,by=0.2)), edim = 2)  # 0.2359326

## Not run: (Careful: This will take several minutes.)
# generate simulated data
N <- 1000; t <- 0.001*(1:N)
sint   <- sin(2*pi*10*t);    sd1 <- sd(sint)    # sine curve
whitet <- rnorm(N);          sd2 <- sd(whitet)  # white noise
chirpt <- sint + 0.1*whitet; sd3 <- sd(chirpt)  # chirp signal

# calculate approximate entropy
rnum <- 30; result <- zeros(3, rnum)
for (i in 1:rnum) {
    r <- 0.02 * i
    result[1, i] <- approx_entropy(sint,   2, r*sd1)
    result[2, i] <- approx_entropy(chirpt, 2, r*sd2)
    result[3, i] <- approx_entropy(whitet, 2, r*sd3)
}

# plot curves
r <- 0.02 * (1:rnum)
plot(c(0, 0.6), c(0, 2), type="n",
     xlab = "", ylab = "", main = "Approximate Entropy")
points(r, result[1, ], col="red");    lines(r, result[1, ], col="red")
points(r, result[2, ], col="green");  lines(r, result[2, ], col="green")
points(r, result[3, ], col="blue");   lines(r, result[3, ], col="blue")
grid()
## End(Not run)

Arc Length of a Curve

Description

Calculates the arc length of a parametrized curve.

Usage

arclength(f, a, b, nmax = 20, tol = 1e-05, ...)

Arguments

Details

Calculates the arc length of a parametrized curve in R^n. It applies Richardson's extrapolation by refining polygon approximations to the curve.

The parametrization of the curve must be vectorized: if t-->F(t) is the parametrization, F(c(t1,t1,...)) must return c(F(t1),F(t2),...).

Can be directly applied to determine the arc length of a one-dimensional function f:R-->R by defining F (if f is vectorized) as F:t-->c(t,f(t)).

Value

Returns a list with components length the calculated arc length, niter the number of iterations, and rel.err the relative error generated from the extrapolation.

Note

If by chance certain equidistant points of the curve lie on a straight line, the result may be wrong, then use polylength below.

Author(s)

HwB <hwborchers@googlemail.com>

See Also

poly_length

Examples

##  Example: parametrized 3D-curve with t in 0..3*pi
f <- function(t) c(sin(2*t), cos(t), t)
arclength(f, 0, 3*pi)
# $length:  17.22203            # true length 17.222032...

##  Example: length of the sine curve
f <- function(t) c(t, sin(t))
arclength(f, 0, pi)             # true length  3.82019...

## Example: Length of an ellipse with axes a = 1 and b = 0.5
# parametrization x = a*cos(t), y = b*sin(t)
a <- 1.0; b <- 0.5
f <- function(t) c(a*cos(t), b*sin(t))
L <- arclength(f, 0, 2*pi, tol = 1e-10)     #=> 4.84422411027
# compare with elliptic integral of the second kind
e <- sqrt(1 - b^2/a^2)                      # ellipticity
L <- 4 * a * ellipke(e^2)$e                 #=> 4.84422411027

## Not run: 
##  Example: oscillating 1-dimensional function (from 0 to 5)
f <- function(x) x * cos(0.1*exp(x)) * sin(0.1*pi*exp(x))
F <- function(t) c(t, f(t))
L <- arclength(F, 0, 5, tol = 1e-12, nmax = 25)
print(L$length, digits = 16)
# [1] 82.81020372882217         # true length 82.810203728822172...

# Split this computation in 10 steps (run time drops from 2 to 0.2 secs)
L <- 0
for (i in 1:10)
	L <- L + arclength(F, (i-1)*0.5, i*0.5, tol = 1e-10)$length
print(L, digits = 16)
# [1] 82.81020372882216

# Alternative calculation of arc length
f1 <- function(x) sqrt(1 + complexstep(f, x)^2)
L1 <- quadgk(f1, 0, 5, tol = 1e-14)
print(L1, digits = 16)
# [1] 82.81020372882216
  
## End(Not run)

## Not run: 
#-- --------------------------------------------------------------------
#   Arc-length parametrization of Fermat's spiral
#-- --------------------------------------------------------------------
# Fermat's spiral: r = a * sqrt(t) 
f <- function(t) 0.25 * sqrt(t) * c(cos(t), sin(t))

t1 <- 0; t2 <- 6*pi
a  <- 0; b  <- arclength(f, t1, t2)$length
fParam <- function(w) {
    fct <- function(u) arclength(f, a, u)$length - w
    urt <- uniroot(fct, c(a, 6*pi))
    urt$root
}

ts <- linspace(0, 6*pi, 250)
plot(matrix(f(ts), ncol=2), type='l', col="blue", 
     asp=1, xlab="", ylab = "",
     main = "Fermat's Spiral", sub="20 subparts of equal length")

for (i in seq(0.05, 0.95, by=0.05)) {
    v <- fParam(i*b); fv <- f(v)
    points(fv[1], f(v)[2], col="darkred", pch=20)
} 
## End(Not run)

Arnoldi Iteration

Description

Arnoldi iteration generates an orthonormal basis of the Krylov space and a Hessenberg matrix.

Usage

arnoldi(A, q, m)

Arguments

Details

arnoldi(A, q, m) carries out m iterations of the Arnoldi iteration with n-by-n matrix A and starting vector q (which need not have unit 2-norm). For m < n it produces an n-by-(m+1) matrix Q with orthonormal columns and an (m+1)-by-m upper Hessenberg matrix H such that A*Q[,1:m] = Q[,1:m]*H[1:m,1:m] + H[m+1,m]*Q[,m+1]*t(E_m), where E_m is the m-th column of the m-by-m identity matrix.

Value

Returns a list with two elements:

Q A matrix of orthonormal columns that generate the Krylov space (A, A q, A^2 q, ...).

H A Hessenberg matrix such that A = Q * H * t(Q).

References

Nicholas J. Higham (2008). Functions of Matrices: Theory and Computation, SIAM, Philadelphia.

See Also

hessenberg

Examples

A <- matrix(c(-149,   -50,  -154,
               537,   180,   546,
               -27,    -9,   -25), nrow = 3, byrow = TRUE)
a <- arnoldi(A, c(1,0,0))
a
## $Q
##      [,1]       [,2]       [,3]
## [1,]    1  0.0000000  0.0000000
## [2,]    0  0.9987384 -0.0502159
## [3,]    0 -0.0502159 -0.9987384
## 
## $H
##           [,1]         [,2]        [,3]
## [1,] -149.0000 -42.20367124  156.316506
## [2,]  537.6783 152.55114875 -554.927153
## [3,]    0.0000   0.07284727    2.448851

a$Q %*% a$H %*% t(a$Q)
##      [,1] [,2] [,3]
## [1,] -149  -50 -154
## [2,]  537  180  546
## [3,]  -27   -9  -25

Barycentric Lagrange Interpolation

Description

Barycentric Lagrange interpolation in one dimension.

Usage

barylag(xi, yi, x)

Arguments

Details

barylag interpolates the given data using the barycentric Lagrange interpolation formula (vectorized to remove all loops).

Value

Values of interpolated data at points x.

Note

Barycentric interpolation is preferred because of its numerical stability.

References

Berrut, J.-P., and L. Nick Trefethen (2004). “Barycentric Lagrange Interpolation”. SIAM Review, Vol. 46(3), pp.501–517.

See Also

Lagrange or Newton interpolation.

Examples

##  Generates an example with plot.
# Input:
#   fun  ---  function that shall be 'approximated'
#   a, b ---  interval [a, b] to be used for the example
#   n    ---  number of supporting nodes
#   m    ---  number of interpolation points
# Output
#   plot of function, interpolation, and nodes
#   return value is NULL (invisible)
## Not run: 
barycentricExample <- function(fun, a, b, n, m)
{
	xi <- seq(a, b, len=n)
	yi <- fun(xi)
	x  <- seq(a, b, len=m)

	y <- barylag(xi, yi, x)
	plot(xi, yi, col="red", xlab="x", ylab="y",
		main="Example of barycentric interpolation")

	lines(x, fun(x), col="yellow", lwd=2)
	lines(x, y, col="darkred")

	grid()
}

barycentricExample(sin, -pi, pi, 11, 101)  # good interpolation
barycentricExample(runge, -1, 1, 21, 101)  # bad interpolation

## End(Not run)

2-D Barycentric Lagrange Interpolation

Description

Two-dimensional barycentric Lagrange interpolation.

Usage

barylag2d(F, xn, yn, xf, yf)

Arguments

Details

Well-known Lagrange interpolation using barycentric coordinates, here extended to two dimensions. The function is completely vectorized.

x-coordinates run downwards in F, y-coordinates to the right. That conforms to the usage in image or contour plots, see the example below.

Value

Matrix of size length(xf)-by-length(yf) giving the interpolated values at al the grid points (xf, yf).

Note

Copyright (c) 2004 Greg von Winckel of a Matlab function under BSD license; translation to R by Hans W Borchers with permission.

References

Berrut, J.-P., and L. Nick Trefethen (2004). “Barycentric Lagrange Interpolation”. SIAM Review, Vol. 46(3), pp.501–517.

See Also

interp2, barylag

Examples

##  Example from R-help
xn <- c(4.05, 4.10, 4.15, 4.20, 4.25, 4.30, 4.35)
yn <- c(60.0, 67.5, 75.0, 82.5, 90.0)
foo <- matrix(c(
        -137.8379, -158.8240, -165.4389, -166.4026, -166.2593,
        -152.1720, -167.3145, -171.1368, -170.9200, -170.4605,
        -162.2264, -172.5862, -174.1460, -172.9923, -172.2861,
        -168.7746, -175.2218, -174.9667, -173.0803, -172.1853,
        -172.4453, -175.7163, -174.0223, -171.5739, -170.5384,
        -173.7736, -174.4891, -171.6713, -168.8025, -167.6662,
        -173.2124, -171.8940, -168.2149, -165.0431, -163.8390),
            nrow = 7, ncol = 5, byrow = TRUE)
xf <- c(4.075, 4.1)
yf <- c(63.75, 67.25)
barylag2d(foo, xn, yn, xf, yf)
#  -156.7964 -163.1753
#  -161.7495 -167.0424

# Find the minimum of the underlying function
bar <- function(xy) barylag2d(foo, xn, yn, xy[1], xy[2])
optim(c(4.25, 67.5), bar)  # "Nelder-Mead"
# $par
# 4.230547 68.522747
# $value
# -175.7959

## Not run: 
# Image and contour plots
image(xn, yn, foo)
contour(xn, yn, foo, col="white", add = TRUE)
xs <- seq(4.05, 4.35, length.out = 51)
ys <- seq(60.0, 90.0, length.out = 51)
zz <- barylag2d(foo, xn, yn, xs, ys)
contour(xs, ys, zz, nlevels = 20, add = TRUE)
contour(xs, ys, zz, levels=c(-175, -175.5), add = TRUE)
points(4.23, 68.52)
## End(Not run)

Bernoulli Numbers and Polynomials

Description

The Bernoulli numbers are a sequence of rational numbers that play an important role for the series expansion of hyperbolic functions, in the Euler-MacLaurin formula, or for certain values of Riemann's function at negative integers.

Usage

bernoulli(n, x)

Arguments

Details

The calculation of the Bernoulli numbers uses the values of the zeta function at negative integers, i.e. B_n = -n \, zeta(1-n). Bernoulli numbers B_n for odd n are 0 except B_1 which is set to -0.5 on purpose.

The Bernoulli polynomials can be directly defined as

B_n(x) = \sum_{k=0}^n {n \choose k} b_{n-k}\, x^k

and it is immediately clear that the Bernoulli numbers are then given as B_n = B_n(0).

Value

Returns the first n+1 Bernoulli numbers, if x is missing, or the value of the Bernoulli polynomial at point(s) x.

Note

The definition uses B_1 = -1/2 in accordance with the definition of the Bernoulli polynomials.

References

See the entry on Bernoulli numbers in the Wikipedia.

See Also

zeta

Examples

bernoulli(10)
# 1.00000000 -0.50000000  0.16666667  0.00000000 -0.03333333
# 0.00000000  0.02380952  0.00000000 -0.03333333  0.00000000  0.07575758
                #
## Not run: 
x1 <- linspace(0.3, 0.7, 2)
y1 <- bernoulli(1, x1)
plot(x1, y1, type='l', col='red', lwd=2,
     xlim=c(0.0, 1.0), ylim=c(-0.2, 0.2),
     xlab="", ylab="", main="Bernoulli Polynomials")
grid()
xs <- linspace(0, 1, 51)
lines(xs, bernoulli(2, xs), col="green", lwd=2)
lines(xs, bernoulli(3, xs), col="blue", lwd=2)
lines(xs, bernoulli(4, xs), col="cyan", lwd=2)
lines(xs, bernoulli(5, xs), col="brown", lwd=2)
lines(xs, bernoulli(6, xs), col="magenta", lwd=2)
legend(0.75, 0.2, c("B_1", "B_2", "B_3", "B_4", "B_5", "B_6"),
       col=c("red", "green", "blue", "cyan", "brown", "magenta"),
       lty=1, lwd=2)
  
## End(Not run)

Bernstein Polynomials

Description

Bernstein base polynomials and approximations.

Usage

bernstein(f, n, x)

bernsteinb(k, n, x)

Arguments

Details

The Bernstein basis polynomials B_{k,n}(x) are defined as

B_{k,n}(x) = {{n}\choose{k}} x^k (1-x)^{n-k}

and form a basis for the vector space of polynomials of degree n over the interval [0,1].

bernstein(f, n, x) computes the approximation of function f through Bernstein polynomials of degree n, resp. computes the value of this approximation at x. The function is vectorized and applies a brute force calculation.

But if x is a scalar, the value will be calculated using De Casteljau's algorithm for higher accuracy. For bigger n the binomial coefficients may be in for problems.

Value

Returns a scalar or vector of function values.

References

See https://en.wikipedia.org/wiki/Bernstein_polynomial

Examples

## Example
f <- function(x) sin(2*pi*x)
xs <- linspace(0, 1)
ys <- f(xs)
## Not run: 
plot(xs, ys, type='l', col="blue",
     main="Bernstein Polynomials")
grid()
b10  <- bernstein(f,  10, xs)
b100 <- bernstein(f, 100, xs)
lines(xs, b10,  col="magenta")
lines(xs, b100, col="red") 
## End(Not run)

# Bernstein basis polynomials
## Not run: 
xs <- linspace(0, 1)
plot(c(0,1), c(0,1), type='n',
     main="Bernstein Basis Polynomials")
grid()
n = 10
for (i in 0:n) {
    bs <- bernsteinb(i, n, xs)
    lines(xs, bs, col=i+1)
} 
## End(Not run)

Rootfinding Through Bisection or Secant Rule

Description

Finding roots of univariate functions in bounded intervals.

Usage

bisect(fun, a, b, maxiter = 500, tol = NA, ...)

secant(fun, a, b, maxiter = 500, tol = 1e-08, ...)

regulaFalsi(fun, a, b, maxiter = 500, tol = 1e-08, ...)

Arguments

Details

“Bisection” is a well known root finding algorithms for real, univariate, continuous functions. Bisection works in any case if the function has opposite signs at the endpoints of the interval.

bisect stops when floating point precision is reached, attaching a tolerance is no longer needed. This version is trimmed for exactness, not speed. Special care is taken when 0.0 is a root of the function. Argument 'tol' is deprecated and not used anymore.

The “Secant rule” uses a succession of roots of secant lines to better approximate a root of a function. “Regula falsi” combines bisection and secant methods. The so-called ‘Illinois’ improvement is used here.

Value

Return a list with components root, f.root, the function value at the found root, iter, the number of iterations done, and root, and the estimated accuracy estim.prec

References

Quarteroni, A., R. Sacco, and F. Saleri (2007). Numerical Mathematics. Second Edition, Springer-Verlag, Berlin Heidelberg.

See Also

ridders

Examples

bisect(sin, 3.0, 4.0)
# $root             $f.root             $iter   $estim.prec
# 3.1415926536      1.2246467991e-16    52      4.4408920985e-16

bisect(sin, -1.0, 1.0)
# $root             $f.root             $iter   $estim.prec
# 0                 0                   2       0

# Legendre polynomial of degree 5
lp5 <- c(63, 0, -70, 0, 15, 0)/8
f <- function(x) polyval(lp5, x)
bisect(f, 0.6, 1)       # 0.9061798453      correct to 15 decimals
secant(f, 0.6, 1)       # 0.5384693         different root
regulaFalsi(f, 0.6, 1)  # 0.9061798459      correct to 10 decimals

Binary Representation

Description

Literal bit representation.

Usage

bits(x, k = 54, pos_sign = FALSE, break0 = FALSE)

Arguments

Details

The literal bit/binary representation of a floating point number is computed by subtracting powers of 2.

Value

Returns a string containing the binary representation.

See Also

nextpow2

Examples

bits(2^10)        # "10000000000"
bits(1 + 2^-10)   #  "1.000000000100000000000000000000000000000000000000000000"
bits(pi)          # "11.001001000011111101101010100010001000010110100011000000"
bits(1/3.0)       #  "0.010101010101010101010101010101010101010101010101010101"
bits(1 + eps())   #  "1.000000000000000000000000000000000000000000000000000100"

String of Blank Carakters

Description

Create a string of blank characters.

Usage

blanks(n)

Arguments

Details

blanks(n) is a string of n blanks.

Value

String of n blanks.

See Also

deblank

Examples

blanks(6)

Block Diagonal Matrix

Description

Build a block diagonal matrix.

Usage

blkdiag(...)

Arguments

Details

Generate a block diagonal matrix from A, B, C, .... All the arguments must be numeric and non-empty matrices.

Value

a numeric matrix

Note

Vectors as input have to be converted to matrices before. Note that as.matrix(v) with v a vector will generate a column vector; use matrix(v, nrow=1) if a row vector is intended.

See Also

Diag

Examples

a1 <- matrix(c(1,2), 1)
a2 <- as.matrix(c(1,2))
blkdiag(a1, diag(1, 2, 2), a2)

Brent-Dekker Root Finding Algorithm

Description

Find root of continuous function of one variable.

Usage

brentDekker(fun, a, b, maxiter = 500, tol = 1e-12, ...)
brent(fun, a, b, maxiter = 500, tol = 1e-12, ...)

Arguments

Details

brentDekker implements a version of the Brent-Dekker algorithm, a well known root finding algorithms for real, univariate, continuous functions. The Brent-Dekker approach is a clever combination of secant and bisection with quadratic interpolation.

brent is simply an alias for brentDekker.

Value

brent returns a list with

References

Quarteroni, A., R. Sacco, and F. Saleri (2007). Numerical Mathematics. Second Edition, Springer-Verlag, Berlin Heidelberg.

See Also

ridders, newtonRaphson

Examples

# Legendre polynomial of degree 5
lp5 <- c(63, 0, -70, 0, 15, 0)/8
f <- function(x) polyval(lp5, x)
brent(f, 0.6, 1)                # 0.9061798459 correct to 12 places

Brownian Motion

Description

The Brown72 data set represents a fractal Brownian motion with a prescribed Hurst exponent 0f 0.72 .

Usage

data(brown72)

Format

The format is: one column.

Details

Estimating the Hurst exponent for a data set provides a measure of whether the data is a pure random walk or has underlying trends. Brownian walks can be generated from a defined Hurst exponent.

Examples

## Not run: 
data(brown72)
plot(brown72, type = "l", col = "blue")
grid()
## End(Not run)

Broyden's Method

Description

Broyden's method for the numerical solution of nonlinear systems of n equations in n variables.

Usage

broyden(Ffun, x0, J0 = NULL, ...,
        maxiter = 100, tol = .Machine$double.eps^(1/2))

Arguments

Details

F as a function must return a vector of length n, and accept an n-dim. vector or column vector as input. F must not be univariate, that is n must be greater than 1.

Broyden's method computes the Jacobian and its inverse only at the first iteration, and does a rank-one update thereafter, applying the so-called Sherman-Morrison formula that computes the inverse of the sum of an invertible matrix A and the dyadic product, uv', of a column vector u and a row vector v'.

Value

List with components: zero the best root found so far, fnorm the square root of sum of squares of the values of f, and niter the number of iterations needed.

Note

Applied to a system of n linear equations it will stop in 2n steps

References

Quarteroni, A., R. Sacco, and F. Saleri (2007). Numerical Mathematics. Second Edition, Springer-Verlag, Berlin Heidelberg.

See Also

newtonsys, fsolve

Examples

##  Example from Quarteroni & Saleri
F1 <- function(x) c(x[1]^2 + x[2]^2 - 1, sin(pi*x[1]/2) + x[2]^3)
broyden(F1, x0 = c(1, 1))
# zero: 0.4760958 -0.8793934; fnorm: 9.092626e-09; niter: 13

F <- function(x) {
    x1 <- x[1]; x2 <- x[2]; x3 <- x[3]
    as.matrix(c(x1^2 + x2^2 + x3^2 - 1,
                x1^2 + x3^2 - 0.25,
                x1^2 + x2^2 - 4*x3), ncol = 1)
}
x0 <- as.matrix(c(1, 1, 1))
broyden(F, x0)
# zero: 0.4407629 0.8660254 0.2360680; fnorm: 1.34325e-08; niter: 8

##  Find the roots of the complex function sin(z)^2 + sqrt(z) - log(z)
F2 <- function(x) {
    z  <- x[1] + x[2]*1i
    fz <- sin(z)^2 + sqrt(z) - log(z)
    c(Re(fz), Im(fz))
}
broyden(F2, c(1, 1))
# zero   0.2555197 0.8948303 , i.e.  z0 = 0.2555 + 0.8948i
# fnorm  7.284374e-10
# niter  13

##  Two more problematic examples
F3 <- function(x)
        c(2*x[1] - x[2] - exp(-x[1]), -x[1] + 2*x[2] - exp(-x[2]))
broyden(F3, c(0, 0))
# $zero   0.5671433 0.5671433   # x = exp(-x)

F4 <- function(x)   # Dennis Schnabel
        c(x[1]^2 + x[2]^2 - 2, exp(x[1] - 1) + x[2]^3 - 2)
broyden(F4, c(2.0, 0.5), maxiter = 100)

Elementwise Function Application (Matlab Style)

Description

Apply a binary function elementwise.

Usage

  bsxfun(func, x, y)

  arrayfun(func, ...)

Arguments

Details

bsxfun applies element-by-element a binary function to two vectors, matrices, or arrays of the same size. For matrices, sweep is used for reasons of speed, otherwise mapply. (For arrays of more than two dimensions this may become very slow.)

arrayfun applies func to each element of the arrays and returns an array of the same size.

Value

The result will be a vector or matrix of the same size as x, y.

Note

The underlying function mapply can be applied in a more general setting with many function parameters:

mapply(f, x, y, z, ...)

but the array structure will not be preserved in this case.

See Also

Vectorize

Examples

X <- matrix(rep(1:10, each = 10), 10, 10)
Y <- t(X)
bsxfun("*", X, Y)  # multiplication table

f <- function(x, y) x[1] * y[1]     # function not vectorized
A <- matrix(c(2, 3, 5, 7), 2, 2)
B <- matrix(c(11, 13, 17, 19), 2, 2)
arrayfun(f, A, B)

Bulirsch-Stoer Algorithm

Description

Bulirsch-Stoer algorithm for solving Ordinary Differential Equations (ODEs) very accurately.

Usage

bulirsch_stoer(f, t, y0, ..., tol = 1e-07)

midpoint(f, t0, tfinal, y0, tol = 1e-07, kmax = 25)

Arguments

Details

The Bulirsch-Stoer algorithm is a well-known method to obtain high-accuracy solutions to ordinary differential equations with reasonable computational efforts. It exploits the midpoint method to get good accuracy in each step.

The (modified) midpoint method computes the values of the dependent variable y(t) from t0 to tfinal by a sequence of substeps, applying Richardson extrapolation in each step.

Bulirsch-Stoer and midpoint shall not be used with non-smooth functions or singularities inside the interval. The best way to get intermediate points t = (t[1], ..., t[n]) may be to call ode23 or ode23s first and use the x-values returned to start bulirsch_stoer on.

Value

bulirsch_stoer returns a list with x the grid points input, and y a vector of function values at the se points.

Note

Will be extended to become a full-blown Bulirsch-Stoer implementation.

Author(s)

Copyright (c) 2014 Hans W Borchers

References

J. Stoer and R. Bulirsch (2002). Introduction to Numerical Analysis. Third Edition, Texts in Applied Mathematics 12, Springer Science + Business, LCC, New York.

See Also

ode23, ode23s

Examples

## Example: y'' = -y
f1 <- function(t, y) as.matrix(c(y[2], -y[1]))
y0 <- as.matrix(c(0.0, 1.0))
tt <- linspace(0, pi, 13)
yy <- bulirsch_stoer(f1, tt, c(0.0, 1.0))   # 13 equally-spaced grid points
yy[nrow(yy), 1]                             # 1.1e-11

## Not run: 
S  <- ode23(f1, 0, pi, c(0.0, 1.0))
yy <- bulirsch_stoer(f1, S$t, c(0.0, 1.0))  # S$x 13 irregular grid points
yy[nrow(yy), 1]                             #  2.5e-11
S$y[nrow(S$y), 1]                           # -7.1e-04

## Example: y' = -200 x y^2                 # y(x) = 1 / (1 + 100 x^2)
f2 <- function(t, y) -200 * t * y^2
y0 < 1
tic(); S <- ode23(f2, 0, 1, y0); toc()            # 0.002 sec
tic(); yy <- bulirsch_stoer(f2, S$t, y0); toc()   # 0.013 sec
## End(Not run)

Boundary Value Problems

Description

Solves boundary value problems of linear second order differential equations.

Usage

bvp(f, g, h, x, y, n = 50)

Arguments

Details

Solves the two-point boundary value problem given as a linear differential equation of second order in the form:

y'' = f(x) y' + g(x) y + h(x)

with the finite element method. The solution y(x) shall exist on the interval [a, b] with boundary conditions y(a) = y_a and y(b) = y_b.

Value

Returns a list list(xs, ys) with the grid points xs and the values ys of the solution at these points, including the boundary points.

Note

Uses a tridiagonal equation solver that may be faster then qr.solve for large values of n.

References

Kutz, J. N. (2005). Practical Scientific Computing. Lecture Notes 98195-2420, University of Washington, Seattle.

See Also

shooting

Examples

##  Solve y'' = 2*x/(1+x^2)*y' - 2/(1+x^2) * y + 1
##  with y(0) = 1.25 and y(4) = -0.95 on the interval [0, 4]:
f1 <- function(x) 2*x / (1 + x^2)
f2 <- function(x)  -2 / (1 + x^2)
f3 <- function(x) rep(1, length(x))     # vectorized constant function 1
x <- c(0.0,   4.0)
y <- c(1.25, -0.95)
sol <- bvp(f1, f2, f3, x, y)
## Not run: 
plot(sol$xs, sol$ys, ylim = c(-2, 2),
     xlab = "", ylab = "", main = "Boundary Value Problem")
# The analytic solution is
sfun <- function(x) 1.25 + 0.4860896526*x - 2.25*x^2 + 
                    2*x*atan(x) - 1/2 * log(1+x^2) + 1/2 * x^2 * log(1+x^2)
xx <- linspace(0, 4)
yy <- sfun(xx)
lines(xx, yy, col="red")
grid()
## End(Not run)

Coordinate Transformations

Description

Transforms between cartesian, spherical, polar, and cylindrical coordinate systems in two and three dimensions.

Usage

cart2sph(xyz)
sph2cart(tpr)
cart2pol(xyz)
pol2cart(prz)

Arguments

Details

The spherical coordinate system used here consists of

- theta, azimuth angle relative to the positive x-axis

- phi, elevation angle measured from the reference plane

- r, radial distance. i.e., distance to the origin

The polar angle, measured with respect from the polar axis, is then calculated as pi/2 - phi. Note that this convention differs slightly from spherical coordinates (r, theta, phi) as often used in mathematics, where phi is the polar angle.

cart2sph returns spherical coordinates as (theta, phi, r), and sph2cart expects them in this sequence.

cart2pol returns polar coordinates (phi, r) if length(xyz)==2 and cylindrical coordinates (phi, r, z) else. pol2cart needs them in this sequence and length.

To go from cylindrical to cartesian coordinates, transform to cartesian coordinates first — or write your own function, see the examples.

All transformation functions are vectorized.

Value

All functions return a (2- or 3-dimensional) vector representing a point in the requested coordinate system, or a matrix with 2 or 3 named columns where is row represents a point. The columns are named accordingly.

Note

In Matlab these functions accept two or three variables and return two or three values. In R it did not appear appropriate to return coordinates as a list.

These functions should be vectorized in the sense that they accept will accept matrices with number of rows or columns equal to 2 or 3.

Examples

x <- 0.5*cos(pi/6); y <- 0.5*sin(pi/6); z <- sqrt(1 - x^2 - y^2)
(s <-cart2sph(c(x, y, z)))      # 0.5235988 1.0471976 1.0000000
sph2cart(s)                     # 0.4330127 0.2500000 0.8660254

cart2pol(c(1,1))                # 0.7853982 1.4142136
cart2pol(c(1,1,0))              # 0.7853982 1.4142136 0.0000000
pol2cart(c(pi/2, 1))            # 6.123234e-17 1.000000e+00
pol2cart(c(pi/4, 1, 1))         # 0.7071068 0.7071068 1.0000000

##  Transform spherical to cylindrical coordinates and vice versa
#   sph2cyl <- function(th.ph.r) cart2pol(sph2cart(th.ph.r))
#   cyl2sph <- function(phi.r.z) cart2sph(pol2cart(phi.r.z))

Directory Functions (Matlab style)

Description

Displays or changes working directory, or lists files therein.

Usage

cd(dname)
pwd()

what(dname = getwd())

Arguments

Details

pwd() displays the name of the current directory, and is the same as cd(). cd(dname) changes to directory dname and if successfull displays the directory name.

what() lists all files in a directory.

Value

Name of the current working directory.

See Also

getwd, setwd, list.files

Examples

# cd()
# pwd()
# what()

Integer Functions (Matlab Style)

Description

Functions for rounding and truncating numeric values towards near integer values.

Usage

ceil(n)
Fix(n)

Arguments

Details

ceil() is an alias for ceiling() and rounds to the smallest integer equal to or above n.

Fix() truncates values towards 0 and is an alias for trunc(). Uses ml prefix to indicate Matlab style.

The corresponding functions floor() (rounding to the largest interger equal to or smaller than n) and round() (rounding to the specified number of digits after the decimal point, default being 0) are already part of R base.

Value

integer values

Examples

x <- c(-1.2, -0.8, 0, 0.5, 1.1, 2.9)
ceil(x)
Fix(x)

Characteristic Polynomial

Description

Computes the characteristic polynomial (and the inverse of the matrix, if requested) using the Faddeew-Leverrier method.

Usage

charpoly(a, info = FALSE)

Arguments

Details

Computes the characteristic polynomial recursively. In the last step the determinant and the inverse matrix can be determined without any extra cost (if the matrix is not singular).

Value

Either the characteristic polynomial as numeric vector, or a list with components cp, the characteristic polynomial, det, the determinant, and inv, the inverse matrix, will be returned.

References

Hou, S.-H. (1998). Classroom Note: A Simple Proof of the Leverrier–Faddeev Characteristic Polynomial Algorithm, SIAM Review, 40(3), pp. 706–709.

Examples

a <- magic(5)
A <- charpoly(a, info = TRUE)
A$cp
roots(A$cp)
A$det
zapsmall(A$inv %*% a)

Chebyshev Approximation

Description

Function approximation through Chebyshev polynomials (of the first kind).

Usage

chebApprox(x, fun, a, b, n)

Arguments

Details

Return approximate y-coordinates of points at x by computing the Chebyshev approximation of degree n for fun on the interval [a, b].

Value

A numeric vector of the same length as x.

Note

TODO: Evaluate the Chebyshev approximative polynomial by using the Clenshaw recurrence formula. (Not yet vectorized, that's why we still use the Horner scheme.)

References

Press, W. H., S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery (1992). Numerical Recipes in C: The Art of Scientific Computing. Second Edition, Cambridge University Press.

See Also

polyApprox

Examples

# Approximate sin(x) on [-pi, pi] with a polynomial of degree 9 !
# This polynomial has to be beaten:
# P(x) = x - 1/6*x^3 + 1/120*x^5 - 1/5040*x^7 + 1/362880*x^9

# Compare these polynomials
p1 <- rev(c(0, 1, 0, -1/6, 0, 1/120, 0, -1/5040, 0, 1/362880))
p2 <- chebCoeff(sin, -pi, pi, 9)

# Estimate the maximal distance
x  <- seq(-pi, pi, length.out = 101)
ys <- sin(x)
yp <- polyval(p1, x)
yc <- chebApprox(x, sin, -pi, pi, 9)
max(abs(ys-yp))                       # 0.006925271
max(abs(ys-yc))                       # 1.151207e-05

## Not run: 
# Plot the corresponding curves
plot(x, ys, type = "l", col = "gray", lwd = 5)
lines(x, yp, col = "navy")
lines(x, yc, col = "red")
grid()
## End(Not run)

Chebyshev Polynomials

Description

Chebyshev Coefficients for Chebyshev polynomials of the first kind.

Usage

chebCoeff(fun, a, b, n)

Arguments

Details

For a function fun on on the interval [a, b] determines the coefficients of the Chebyshev polynomials up to degree n that will approximate the function (in L2 norm).

Value

Vector of coefficients for the Chebyshev polynomials, from low to high degrees (see the example).

Note

See the “Chebfun Project” <https://www.chebfun.org/> by Nick Trefethen.

References

Weisstein, Eric W. “Chebyshev Polynomial of the First Kind." From MathWorld — A Wolfram Web Resource. https://mathworld.wolfram.com/ChebyshevPolynomialoftheFirstKind.html

See Also

chebPoly, chebApprox

Examples

##  Chebyshev coefficients for x^2 + 1
n <- 4
f2 <- function(x) x^2 + 1
cC <- chebCoeff(f2, -1, 1, n)  #  3.0   0  0.5   0   0
cC[1] <- cC[1]/2               # correcting the absolute Chebyshev term
                               # i.e.  1.5*T_0 + 0.5*T_2
cP <- chebPoly(n)              # summing up the polynomial coefficients
p <- cC %*% cP                 #  0 0 1 0 1

Chebyshev Polynomials

Description

Chebyshev polynomials and their values.

Usage

chebPoly(n, x = NULL)

Arguments

Details

Determines an (n+1)-ny-(n+1)-Matrix of Chebyshev polynomials up to degree n.

The coefficients of the first n Chebyshev polynomials are computed using the recursion formula. For computing any values at points the well known Horner schema is applied.

Value

If x is NULL, returns an (n+1)-by-(n+1) matrix with the coefficients of the first Chebyshev polynomials from 0 to n, one polynomial per row with coefficients from highest to lowest order.

If x is a numeric vector, returns the values of the n-th Chebyshev polynomial at the points of x.

Note

See the “Chebfun Project” <https://www.chebfun.org/> by Nick Trefethen.

References

Carothers, N. L. (1998). A Short Course on Approximation Theory. Bowling Green State University.

See Also

chebCoeff, chebApprox

Examples

chebPoly(6)

## Not run: 
##  Plot 6 Chebyshev Polynomials
plot(0, 0, type="n", xlim=c(-1, 1), ylim=c(-1.2, 1.2),
    main="Chebyshev Polynomials for n=1..6", xlab="x", ylab="y")
grid()
x <- seq(-1, 1, length.out = 101)
for (i in 1:6) {
    y <- chebPoly(i, x)
    lines(x, y, col=i)
}
legend(x = 0.55, y = 1.2, c("n=1", "n=2", "n=3", "n=4", "n=5", "n=6"),
    col = 1:6, lty = 1, bg="whitesmoke", cex = 0.75)

## End(Not run)

Fitting a Circle

Description

Fitting a circle from points in the plane

Usage

circlefit(xp, yp)

Arguments

Details

This routine finds an ‘algebraic’ solution based on a linear fit. The value to be minimized is the distance of the given points to the nearest point on the circle along a radius.

Value

Returns x- and y-coordinates of the center and the radius as a vector of length 3.

Writes the RMS error of the (radial) distance of the original points to the circle directly onto the console.

References

Gander, W., G. H. Golub, and R. Strebel (1994). Fitting of Circles and Ellipses — Least Squares Solutions. ETH Zürich, Technical Report 217, Institut für Wissenschaftliches Rechnen.

Examples

# set.seed(8421)
n  <- 20
w  <- 2*pi*runif(n)
xp <- cos(w) + 1 + 0.25 * (runif(n) - 0.5)
yp <- sin(w) + 1 + 0.25 * (runif(n) - 0.5)

circe <- circlefit(xp, yp)  #=> 0.9899628 1.0044920 1.0256633
                            # RMS error: 0.07631986 
## Not run: 
x0 <- circe[1]; y0 <- circe[2]; r0 <- circe[3]
plot(c(-0.2, 2.2), c(-0.2, 2.2), type="n", asp=1)
grid()
abline(h=0, col="gray"); abline(v=0, col="gray")
points(xp, yp, col="darkred")

w  <- seq(0, 2*pi, len=100)
xx <- r0 * cos(w) + x0
yy <- r0 * sin(w) + y0
lines(xx, yy, col="blue")
## End(Not run)

Clear function (Matlab style)

Description

List or remove items from workspace, or display system information.

Usage

clear(lst)
ver()

who()
whos()

Arguments

Details

Remove these or all items from the workspace, i.e. the global environment, and freeing up system memory.

who() lists all items on the workspace.
whos() lists all items and their class and size.

ver() displays version and license information for R and all the loaded packages.

Value

Invisibly NULL.

See Also

ls, rm, sessionInfo

Examples

# clear()  # DON'T
# who()
# whos()
# ver()

Clenshaw-Curtis Quadrature Formula

Description

Clenshaw-Curtis Quadrature Formula

Usage

clenshaw_curtis(f, a = -1, b = 1, n = 1024, ...)

Arguments

Details

Clenshaw-Curtis quadrature is based on sampling the integrand on Chebyshev points, an operation that can be implemented using the Fast Fourier Transform.

Value

Numerical scalar, the value of the integral.

References

Trefethen, L. N. (2008). Is Gauss Quadrature Better Than Clenshaw-Curtis? SIAM Review, Vol. 50, No. 1, pp 67–87.

See Also

gaussLegendre, gauss_kronrod

Examples

##  Quadrature with Chebyshev nodes and weights
f <- function(x) sin(x+cos(10*exp(x))/3)
## Not run: ezplot(f, -1, 1, fill = TRUE)
cc <- clenshaw_curtis(f, n = 64)  #=>  0.0325036517151 , true error > 1.3e-10

Generate Combinations

Description

Generates all combinations of length m of a vector a.

Usage

combs(a, m)

Arguments

Details

combs generates combinations of length n of the elements of the vector a.

Value

matrix representing combinations of the elements of a

See Also

perms, randcomb

Examples

combs(seq(2, 10, by=2), m = 3)

Companion Matrix

Description

Computes the companion matrix of a real or complex vector.

Usage

  compan(p)

Arguments

Details

Computes the companion matrix corresponding to the vector p with -p[2:length(p)]/p[1] as first row.

The eigenvalues of this matrix are the roots of the polynomial.

Value

A square matrix of length(p)-1 rows and columns

See Also

roots

Examples

  p <- c(1, 0, -7, 6)
  compan(p)
  # 0  7 -6
  # 1  0  0
  # 0  1  0

Complex Step Derivatives

Description

Complex step derivatives of real-valued functions, including gradients, Jacobians, and Hessians.

Usage

complexstep(f, x0, h = 1e-20, ...)

grad_csd(f, x0, h = 1e-20, ...)
jacobian_csd(f, x0, h = 1e-20, ...)
hessian_csd(f, x0, h = 1e-20, ...)
laplacian_csd(f, x0, h = 1e-20, ...)

Arguments

Details

Complex step derivation is a fast and highly exact way of numerically differentiating a function. If the following conditions are satisfied, there will be no loss of accuracy between computing a function value and computing the derivative at a certain point.

• f must have an analytical (i.e., complex differentiable) continuation into an open neighborhood of x0.

• x0 and f(x0) must be real.

• h is real and very small: 0 < h << 1.

complexstep handles differentiation of univariate functions, while grad_csd and jacobian_csd compute gradients and Jacobians by applying the complex step approach iteratively. Please understand that these functions are not vectorized, but complexstep is.

As complex step cannot be applied twice (the first derivative does not fullfil the conditions), hessian_csd works differently. For the first derivation, complex step is used, to the one time derived function Richardson's method is applied. The same applies to lapalacian_csd.

Value

complexstep(f, x0) returns the derivative f'(x_0) of f at x_0. The function is vectorized in x0.

Note

This surprising approach can be easily deduced from the complex-analytic Taylor formula.

Author(s)

HwB <hwborchers@googlemail.com>

References

Martins, J. R. R. A., P. Sturdza, and J. J. Alonso (2003). The Complex-step Derivative Approximation. ACM Transactions on Mathematical Software, Vol. 29, No. 3, pp. 245–262.

See Also

numderiv

Examples

##  Example from Martins et al.
f <- function(x) exp(x)/sqrt(sin(x)^3 + cos(x)^3)  # derivative at x0 = 1.5
# central diff formula    # 4.05342789402801, error 1e-10
# numDeriv::grad(f, 1.5)  # 4.05342789388197, error 1e-12  Richardson
# pracma::numderiv        # 4.05342789389868, error 5e-14  Richardson
complexstep(f, 1.5)       # 4.05342789389862, error 1e-15
# Symbolic calculation:   # 4.05342789389862

jacobian_csd(f, 1.5)

f1 <- function(x) sum(sin(x))
grad_csd(f1, rep(2*pi, 3))
## [1] 1 1 1

laplacian_csd(f1, rep(pi/2, 3))
## [1] -3

f2 <- function(x) c(sin(x[1]) * exp(-x[2]))
hessian_csd(f2, c(0.1, 0.5, 0.9))
##             [,1]        [,2] [,3]
## [1,] -0.06055203 -0.60350053    0
## [2,] -0.60350053  0.06055203    0
## [3,]  0.00000000  0.00000000    0

f3 <- function(u) {
    x <- u[1]; y <- u[2]; z <- u[3]
    matrix(c(exp(x^+y^2), sin(x+y), sin(x)*cos(y), x^2 - y^2), 2, 2)
  }
jacobian_csd(f3, c(1,1,1))
##            [,1]       [,2] [,3]
## [1,]  2.7182818  0.0000000    0
## [2,] -0.4161468 -0.4161468    0
## [3,]  0.2919266 -0.7080734    0
## [4,]  2.0000000 -2.0000000    0

Matrix Condition

Description

Condition number of a matrix.

Usage

cond(M, p = 2)

Arguments

Details

The condition number of a matrix measures the sensitivity of the solution of a system of linear equations to small errors in the data. Values of cond(M) and cond(M, p) near 1 are indications of a well-conditioned matrix.

Value

cond(M) returns the 2-norm condition number, the ratio of the largest singular value of M to the smallest.

c = cond(M, p) returns the matrix condition number in p-norm:

norm(X,p) * norm(inv(X),p).

(Not yet implemented.)

Note

Not feasible for large or sparse matrices as svd(M) needs to be computed. The Matlab/Octave function condest for condition estimation has not been implemented.

References

Trefethen, L. N., and D. Bau III. (1997). Numerical Linear Algebra. SIAM, Philadelphia.

See Also

normest, svd

Examples

cond(hilb(8))

Polynomial Convolution

Description

Convolution and polynomial multiplication.

Usage

conv(x, y)

Arguments

Details

r = conv(p,q) convolves vectors p and q. Algebraically, convolution is the same operation as multiplying the polynomials whose coefficients are the elements of p and q.

Value

Another vector.

Note

conv utilizes fast Fourier transformation.

See Also

deconv, polyadd

Examples

conv(c(1, 1, 1), 1)
conv(c(1, 1, 1), c(0, 0, 1))
conv(c(-0.5, 1, -1), c(0.5, 0, 1))

More Trigonometric Functions

Description

More trigonometric functions not available in R.

Usage

cot(z)
csc(z)
sec(z)
acot(z)
acsc(z)
asec(z)

Arguments

Details

The usual trigonometric cotangens, cosecans, and secans functions and their inverses, computed through the other well known – in R – sine, cosine, and tangens functions.

Value

Result vector of numeric or complex values.

Note

These function names are available in Matlab, that is the reason they have been added to the ‘pracma’ package.

See Also

Trigonometric and hyperbolic functions in R.

Examples

cot(1+1i)       # 0.2176 - 0.8680i
csc(1+1i)       # 0.6215 - 0.3039i
sec(1+1i)       # 0.4983 + 0.5911i
acot(1+1i)      # 0.5536 - 0.4024i
acsc(1+1i)      # 0.4523 - 0.5306i
asec(1+1i)      # 1.1185 + 0.5306i

Newton-Cotes Formulas

Description

Closed composite Newton-Cotes formulas of degree 2 to 8.

Usage

cotes(f, a, b, n, nodes, ...)

Arguments

Details

2 to 8 point closed and summed Newton-Cotes numerical integration formulas.

These formulas are called ‘closed’ as they include the endpoints. They are called ‘composite’ insofar as they are combined with a Lagrange interpolation over subintervals.

Value

The integral as a scalar.

Note

It is generally recommended not to apply Newton-Cotes formula of degrees higher than 6, instead increase the number n of subintervals used.

Author(s)

Standard Newton-Cotes formulas can be found in every textbook. Copyright (c) 2005 Greg von Winckel of nicely vectorized Matlab code, available from MatlabCentral, for 2 to 11 grid points. R version by Hans W Borchers, with permission.

References

Quarteroni, A., R. Sacco, and F. Saleri (2007). Numerical Mathematics. Second Edition, Springer-Verlag, Berlin Heidelberg.

See Also

simpadpt, trapz

Examples

cotes(sin, 0, pi/2, 20, 2)      # 0.999485905248533
cotes(sin, 0, pi/2, 20, 3)      # 1.000000211546591
cotes(sin, 0, pi/2, 20, 4)      # 1.000000391824184
cotes(sin, 0, pi/2, 20, 5)      # 0.999999999501637
cotes(sin, 0, pi/2, 20, 6)      # 0.999999998927507
cotes(sin, 0, pi/2, 20, 7)      # 1.000000000000363  odd degree is better
cotes(sin, 0, pi/2, 20, 8)      # 1.000000000002231

More Hyperbolic Functions

Description

More hyperbolic functions not available in R.

Usage

coth(z)
csch(z)
sech(z)
acoth(z)
acsch(z)
asech(z)

Arguments

Details

The usual hyperbolic cotangens, cosecans, and secans functions and their inverses, computed through the other well known – in R – hyperbolic sine, cosine, and tangens functions.

Value

Result vector of numeric or complex values.

Note

These function names are available in Matlab, that is the reason they have been added to the ‘pracma’ package.

See Also

Trigonometric and hyperbolic functions in R.

Examples

coth(1+1i)      # 0.8680 - 0.2176i
csch(1+1i)      # 0.3039 - 0.6215i
sech(1+1i)      # 0.4983 - 0.5911i
acoth(1+1i)     # 0.4024 - 0.5536i
acsch(1+1i)     # 0.5306 - 0.4523i
asech(1+1i)     # 0.5306 - 1.1185i

Crank-Nicolson Method

Description

The Crank-Nicolson method for solving ordinary differential equations is a combination of the generic steps of the forward and backward Euler methods.

Usage

cranknic(f, t0, t1, y0, ..., N = 100)

Arguments

Details

Adding together forward and backword Euler method in the cranknic method is by finding the root of the function merging these two formulas.

No attempt is made to catch any errors in the root finding functions.

Value

List with components t for grid (or ‘time’) points between t0 and t1, and y an n-by-m matrix with solution variables in columns, i.e. each row contains one time stamp.

Note

This is for demonstration purposes only; for real problems or applications please use ode23 or rkf54.

References

Quarteroni, A., and F. Saleri (2006). Scientific Computing With MATLAB and Octave. Second Edition, Springer-Verlag, Berlin Heidelberg.

See Also

ode23, newmark

Examples

##  Newton's example
f <- function(x, y) 1 - 3*x + y + x^2 + x*y
sol100  <- cranknic(f, 0, 1, 0, N = 100)
sol1000 <- cranknic(f, 0, 1, 0, N = 1000)

## Not run: 
# Euler's forward approach
feuler <- function(f, t0, t1, y0, n) {
    h <- (t1 - t0)/n;  x <- seq(t0, t1, by = h)
    y <- numeric(n+1); y[1] <- y0
    for (i in 1:n) y[i+1] <- y[i] + h * f(x[i], y[i])
    return(list(x = x, y = y))
}

solode <- ode23(f, 0, 1, 0)
soleul <- feuler(f, 0, 1, 0, 100)

plot(soleul$x, soleul$y, type = "l", col = "blue", 
     xlab = "", ylab = "", main = "Newton's example")
lines(solode$t, solode$y, col = "gray", lwd = 3)
lines(sol100$t, sol100$y, col = "red")
lines(sol1000$t, sol1000$y, col = "green")
grid()

##  System of differential equations
# "Herr und Hund"
fhh <- function(x, y) {
    y1 <- y[1]; y2 <- y[2]
    s <- sqrt(y1^2 + y2^2)
    dy1 <- 0.5 - 0.5*y1/s
    dy2 <- -0.5*y2/s
    return(c(dy1, dy2))
}

sol <- cranknic(fhh, 0, 60, c(0, 10))
plot(sol$y[, 1], sol$y[, 2], type = "l", col = "blue",
     xlab = "", ylab = "", main = '"Herr und Hund"')
grid()
## End(Not run)

Vector Cross Product

Description

Vector or cross product

Usage

  cross(x, y)

Arguments

Details

Computes the cross (or: vector) product of vectors in 3 dimensions. In case of matrices it takes the first dimension of length 3 and computes the cross product between corresponding columns or rows.

The more general cross product of n-1 vectors in n-dimensional space is realized as crossn.

Value

3-dim. vector if x and < are vectors, a matrix of 3-dim. vectors if x and y are matrices themselves.

See Also

dot, crossn

Examples

  cross(c(1, 2, 3), c(4, 5, 6))  # -3  6 -3

n-dimensional Vector Cross Product

Description

Vector cross product of n-1 vectors in n-dimensional space

Usage

crossn(A)

Arguments

Details

The rows of the matrix A are taken as(n-1) vectors in n-dimensional space. The cross product generates a vector in this space that is orthogonal to all these rows in A and its length is the volume of the geometric hypercube spanned by the vectors.

Value

a vector of length n

Note

The ‘scalar triple product’ in R^3 can be defined as

spatproduct <- function(a, b, c) dot(a, crossn(b, c))

It represents the volume of the parallelepiped spanned by the three vectors.

See Also

cross, dot

Examples

A <- matrix(c(1,0,0, 0,1,0), nrow=2, ncol=3, byrow=TRUE)
crossn(A)  #=> 0 0 1

x <- c(1.0, 0.0, 0.0)
y <- c(1.0, 0.5, 0.0)
z <- c(0.0, 0.0, 1.0)
identical(dot(x, crossn(rbind(y, z))), det(rbind(x, y, z)))

Interpolating Cubic Spline

Description

Computes the natural interpolation cubic spline.

Usage

cubicspline(x, y, xi = NULL, endp2nd = FALSE, der = c(0, 0))

Arguments

Details

cubicspline computes the values at xi of the natural interpolating cubic spline that interpolate the values y at the nodes x. The derivatives at the endpoints can be prescribed.

Value

Returns either the interpolated values at the points xi or, if is.null(xi), the piecewise polynomial that represents the spline.

Note

From the piecewise polynomial returned one can easily generate the spline function, see the examples.

References

Quarteroni, Q., and F. Saleri (2006). Scientific Computing with Matlab and Octave. Springer-Verlag Berlin Heidelberg.

See Also

spline

Examples

##  Example: Average temperatures at different latitudes
x <- seq(-55, 65, by = 10)
y <- c(-3.25, -3.37, -3.35, -3.20, -3.12, -3.02, -3.02,
       -3.07, -3.17, -3.32, -3.30, -3.22, -3.10)
xs <- seq(-60, 70, by = 1)

# Generate a function for this
pp <- cubicspline(x, y)
ppfun <- function(xs) ppval(pp, xs)

## Not run: 
# Plot with and without endpoint correction
plot(x, y, col = "darkblue",
           xlim = c(-60, 70), ylim = c(-3.5, -2.8),
           xlab = "Latitude", ylab = "Temp. Difference",
           main = "Earth Temperatures per Latitude")
lines(spline(x, y), col = "darkgray")
grid()

ys <- cubicspline(x, y, xs, endp2nd = TRUE)     # der = 0 at endpoints
lines(xs, ys, col = "red")
ys <- cubicspline(x, y, xs)                     # no endpoint condition
lines(xs, ys, col = "darkred")

## End(Not run)

Parametric Curve Fit

Description

Polynomial fitting of parametrized points on 2D curves, also requiring to meet some points exactly.

Usage

curvefit(u, x, y, n, U = NULL, V = NULL)

Arguments

Details

This function will attempt to fit two polynomials to parametrized curve points using the linear least squares approach with linear equality constraints in lsqlin. The requirement to meet exactly some fixed points is interpreted as a linear equality constraint.

Value

Returns a list with 4 components, xp and yp coordinates of the fitted points, and px and py the coefficients of the fitting polynomials in x- and y-direction.

Note

In the same manner, derivatives/directions could be prescribed at certain points.

See Also

circlefit, lsqlin

Examples

##  Approximating half circle arc with small perturbations
N <- 50
u <- linspace(0, pi, N)
x <- cos(u) + 0.05 * randn(1, N)
y <- sin(u) + 0.05 * randn(1, N)
n <- 8
cfit1 <- curvefit(u, x, y, n)
## Not run: 
plot(x, y, col = "darkgray", pch = 19, asp = 1)
xp <- cfit1$xp; yp <- cfit1$yp
lines(xp, yp, col="blue")
grid()
## End(Not run)

##  Fix the end points at t = 0 and t = pi
U <- c(0, pi)
V <- matrix(c(1, 0, -1, 0), 2, 2, byrow = TRUE)
cfit2 <- curvefit(u, x, y, n, U, V)
## Not run: 
xp <- cfit2$xp; yp <- cfit2$yp
lines(xp, yp, col="red")
## End(Not run)

## Not run: 
##  Archimedian spiral
n <- 8
u <- linspace(0, 3*pi, 50)
a <- 1.0
x <- as.matrix(a*u*cos(u))
y <- as.matrix(a*u*sin(u))
plot(x, y, type = "p", pch = 19, col = "darkgray", asp = 1)
lines(x, y, col = "darkgray", lwd = 3)
cfit <- curvefit(u, x, y, n)
px <- c(cfit$px); py <- c(cfit$py)
v <- linspace(0, 3*pi, 200)
xs <- polyval(px, v)
ys <- polyval(py, v)
lines(xs, ys, col = "navy")
grid()
## End(Not run)

Find Cutting Points

Description

Finds cutting points for vector s of real numbers.

Usage

cutpoints(x, nmax = 8, quant = 0.95)

Arguments

Details

Finds cutting points for vector s of real numbers, based on the gaps in the values of the vector. The number of cutting points is derived from a quantile of gaps in the values. The user can set a lower limit for this number of gaps.

Value

Returns a list with components cutp, the cutting points selected, and cutd, the gap between values of x at this cutting point.

Note

Automatically finding cutting points is often requested in Data Mining. If a target attribute is available, Quinlan's C5.0 does a very good job here. Unfortunately, the ‘C5.0’ package (of the R-Forge project “Rulebased Models”) is quite cumbersome to use.

References

Witten, I. H., and E. Frank (2005). Data Mining: Practical Machine Learning Tools and Techniques. Morgan Kaufmann Publishers, San Francisco.

See Also

cut

Examples

N <- 100; x <- sort(runif(N))
cp <- cutpoints(x, 6, 0.9)
n <- length(cp$cutp)

# Print out
nocp <- rle(findInterval(x, c(-Inf, cp$cutp, Inf)))$lengths
cbind(c(-Inf, cp$cutp), c(cp$cutp, Inf), nocp)

# Define a factor from the cutting points
fx <- cut(x, breaks = c(-Inf, cp$cutp, Inf))

## Not run: 
# Plot points and cutting points
plot(x, rep(0, N), col="gray", ann = FALSE)
points(cp$cutp, rep(0, n), pch="|", col=2)

# Compare with k-means clustering
km <- kmeans(x, n)
points(x, rep(0, N), col = km$cluster, pch = "+")

##  A 2-dimensional example
x <- y <- c()
for (i in 1:9) {
  for (j in 1:9) {
    x <- c(x, i + rnorm(20, 0, 0.2))
    y <- c(y, j + rnorm(20, 0, 0.2))
  }
}
cpx <- cutpoints(x, 8, 0)
cpy <- cutpoints(y, 8, 0)

plot(x, y, pch = 18, col=rgb(0.5,0.5,0.5), axes=FALSE, ann=FALSE)
for (xi in cpx$cutp) abline(v=xi, col=2, lty=2)
for (yi in cpy$cutp) abline(h=yi, col=2, lty=2)

km <- kmeans(cbind(x, y), 81)
points(x, y, col=km$cluster)

## End(Not run)

Double and Triple Integration

Description

Numerically evaluate double integral over rectangle.

Usage

dblquad(f, xa, xb, ya, yb, dim = 2, ..., 
        subdivs = 300, tol = .Machine$double.eps^0.5)

triplequad(f, xa, xb, ya, yb, za, zb, 
            subdivs = 300, tol = .Machine$double.eps^0.5, ...)

Arguments

Details

Function dblquad applies the internal single variable integration function integrate two times, once for each variable.

Function triplequad reduces the problem to dblquad by first integrating over the innermost variable.

Value

Numerical scalar, the value of the integral.

See Also

integrate, quad2d, simpson2d

Examples

f1 <- function(x, y) x^2 + y^2
dblquad(f1, -1, 1, -1, 1)       #   2.666666667 , i.e. 8/3 . err = 0

f2 <- function(x, y) y*sin(x)+x*cos(y)
dblquad(f2, pi, 2*pi, 0, pi)    #  -9.869604401 , i.e. -pi^2, err = 0

# f3 <- function(x, y) sqrt((1 - (x^2 + y^2)) * (x^2 + y^2 <= 1))
f3 <- function(x, y) sqrt(pmax(0, 1 - (x^2 + y^2)))
dblquad(f3, -1, 1, -1, 1)       #   2.094395124 , i.e. 2/3*pi , err = 2e-8

f4 <- function(x, y, z) y*sin(x)+z*cos(x)
triplequad(f4, 0,pi, 0,1, -1,1) # - 2.0 => -2.220446e-16

Deconvolution

Description

Deconvolution and polynomial division.

Usage

deconv(b, a)

Arguments

Details

deconv(b,a) deconvolves vector a out of vector b. The quotient is returned in vector q and the remainder in vector r such that b = conv(a,q)+r.

If b and a are vectors of polynomial coefficients, convolving them is equivalent to multiplying the two polynomials, and deconvolution is polynomial division.

Value

List with elements named q and r.

Note

TODO: Base deconv on some filter1d function.

See Also

conv, polymul

Examples

b <- c(10, 40, 100, 160, 170, 120)
a <- c(1, 2, 3, 4)

p <- deconv(b, a)
p$q                #=> 10 20 30
p$r                #=>  0  0  0

Event Detection in ODE solution

Description

Detect events in solutions of a differential equation.

Usage

deeve(x, y, yv = 0, idx = NULL)

Arguments

Details

Determines when (in x coordinates) the idx-th solution function will take on the value yv.

The interpolation is linear for the moment. For points outside the x interval NA is returned.

Value

A (time) point x0 at which the event happens.

Note

The interpolation is linear only for the moment.

See Also

deval

Examples

##  Damped pendulum:  y'' = -0.3 y' - sin(y)
#   y1 = y, y2 = y':  y1' = y2,  y2' = -0.3*y2 - sin(y1)
f <- function(t, y) {
	dy1 <- y[2]
	dy2 <- -0.3*y[2] - sin(y[1])
	return(c(dy1, dy2))
}
sol <- rk4sys(f, 0, 10, c(pi/2, 0), 100)
deeve(sol$x, sol$y[,1])                   # y1 = 0 : elongation in [sec]
# [1] 2.073507 5.414753 8.650250
# matplot(sol$x, sol$y); grid()

Degrees to Radians

Description

Transforms between angles in degrees and radians.

Usage

deg2rad(deg)
rad2deg(rad)

Arguments

Details

This is a simple calculation back and forth. Note that angles greater than 360 degrees are allowed and will be returned. This may appear incorrect but follows a corresponding discussion on Matlab Central.

Value

The angle in degrees or radians.

Examples

deg2rad(c(0, 10, 20, 30, 40, 50, 60, 70, 80, 90))
rad2deg(seq(-pi/2, pi/2, length = 19))

Remove Linear Trends

Description

Removes the mean value or (piecewise) linear trend from a vector or from each column of a matrix.

Usage

detrend(x, tt = 'linear', bp = c())

Arguments

Details

detrend computes the least-squares fit of a straight line (or composite line for piecewise linear trends) to the data and subtracts the resulting function from the data.

To obtain the equation of the straight-line fit, use polyfit.

Value

removes the mean or (piecewise) linear trend from x and returns it in y=detrend(x), that is x-y is the linear trend.

Note

Detrending is often used for FFT processing.

See Also

polyfit

Examples

t <- 1:9
x <- c(0, 2, 0, 4, 4, 4, 0, 2, 0)
x - detrend(x, 'constant')
x - detrend(x, 'linear')

y <- detrend(x, 'linear', 5)
## Not run: 
plot(t, x, col="blue")
lines(t, x - y, col="red")
grid()
## End(Not run)

Evaluate ODE Solution

Description

Evaluate solution of a differential equation solver.

Usage

deval(x, y, xp, idx = NULL)

Arguments

Details

Determines where the points xp lie within the vector x and interpolates linearly.

Value

An length(xp)-by-length(idx) matrix of values at points xp.

Note

The interpolation is linear only for the moment.

See Also

deeve

Examples

##  Free fall:  v' = -g - cw abs(v)^1.1,  cw = 1.6 drag coefficien
f <- function(t, y) -9.81 + 1.6*abs(y)^1.1
sol <- rk4(f, 0, 10, 0, 100)
# speed after 0.5, 1, 1.5, 2 seconds
cbind(c(0.5,1,1.5,2), -deval(sol$x, sol$y, c(0.5, 1, 1.5, 2)))
#  0.5  3.272267  m/s
#  1.0  4.507677
#  1.5  4.953259
#  2.0  5.112068
# plot(sol$x, -sol$y, type="l", col="blue"); grid()

Matrix Diagonal

Description

Generate diagonal matrices or return diagonal of a matrix

Usage

Diag(x, k = 0)

Arguments

Details

If x is a vector, Diag(x, k) generates a matrix with x as the (k-th secondary) diagonal.

If x is a matrix, Diag(x, k) returns the (k-th secondary) diagonal of x.

The k-th secondary diagonal is above the main diagonal for k > 0 and below the main diagonal for k < 0.

Value

matrix or vector

Note

In Matlab/Octave this function is called diag() and has a different signature than the corresponding function in R.

See Also

diag, Trace

Examples

Diag(matrix(1:12,3,4),  1)
Diag(matrix(1:12,3,4), -1)

Diag(c(1,5,9), 1)
Diag(c(1,5,9), -1)

Utility functions (Matlab style)

Description

Display text or array, or produce beep sound.

Usage

disp(...)
beep()

Arguments

Details

Display text or array, or produces the computer's default beep sound using ‘cat’ with closing newline.

Value

beep() returns NULL invisibly, disp() displays with newline.

Examples

disp("Some text, and numbers:", pi, exp(1))
# beep()

Distance Matrix

Description

Computes the Euclidean distance between rows of two matrices.

Usage

distmat(X, Y)
pdist(X)
pdist2(X, Y)

Arguments

Details

Computes Euclidean distance between two vectors A and B as:

||A-B|| = sqrt ( ||A||^2 + ||B||^2 - 2*A.B )

and vectorizes to rows of two matrices (or vectors).

pdist2 is an alias for distmat, while pdist(X) is the same as distmat(X, X).

Value

matrix of size m x n if x is of size m x k and y is of size n x k.

Note

If a is m x r and b is n x r then

apply(outer(a,t(b),"-"),c(1,4),function(x)sqrt(sum(diag(x*x))))

is the m x n matrix of distances between the m rows of a and n rows of b.

This can be modified as necessary, if one wants to apply distances other than the euclidean.

BUT: The code shown here is 10-100 times faster, utilizing the similarity between Euclidean distance and matrix operations.

References

Copyright (c) 1999 Roland Bunschoten for a Matlab version on MatlabCentral under the name distance.m. Translated to R by Hans W Borchers.

See Also

dist

Examples

A <- c(0.0, 0.0)
B <- matrix(c(
        0,0, 1,0, 0,1, 1,1), nrow=4, ncol = 2, byrow = TRUE)
distmat(A, B)  #=> 0 1 1 sqrt(2)

X <- matrix(rep(0.5, 5), nrow=1, ncol=5)
Y <- matrix(runif(50), nrow=10, ncol=5)
distmat(X, Y)

# A more vectorized form of distmat:
distmat2 <- function(x, y) {
    sqrt(outer(rowSums(x^2), rowSums(y^2), '+') - tcrossprod(x, 2 * y))
}

Scalar Product

Description

'dot' or 'scalar' product of vectors or pairwise columns of matrices.

Usage

  dot(x, y)

Arguments

Details

Returns the 'dot' or 'scalar' product of vectors or columns of matrices. Two vectors must be of same length, two matrices must be of the same size. If x and y are column or row vectors, their dot product will be computed as if they were simple vectors.

Value

A scalar or vector of length the number of columns of x and y.

See Also

cross

Examples

  dot(1:5, 1:5)  #=> 55
  # Length of space diagonal in 3-dim- cube:
  sqrt(dot(c(1,1,1), c(1,1,1)))  #=> 1.732051

Eigenvalue Function (Matlab Style)

Description

Eigenvalues of a matrix

Usage

  eig(a)

Arguments

Details

Computes the eigenvalues of a square matrix of real or complex numbers, using the R routine eigen without computing the eigenvectors.

Value

Vector of eigenvalues

See Also

compan

Examples

  eig(matrix(c(1,-1,-1,1), 2, 2))   #=> 2 0
  eig(matrix(c(1,1,-1,1), 2, 2))    # complex values
  eig(matrix(c(0,1i,-1i,0), 2, 2))  # real values

Jacobi Eigenvalue Method

Description

Jacobi's iteration method for eigenvalues and eigenvectors.

Usage

eigjacobi(A, tol = .Machine$double.eps^(2/3))

Arguments

Details

The Jacobi eigenvalue method repeatedly performs (Givens) transformations until the matrix becomes almost diagonal.

Value

Returns a list with components V, a matrix containing the eigenvectors as columns, and D a vector of the eigenvalues.

Note

This R implementation works well up to 50x50-matrices.

References

Mathews, J. H., and K. D. Fink (2004). Numerical Methods Using Matlab. Fourth edition, Pearson education, Inc., New Jersey.

See Also

eig

Examples

A <- matrix(c( 1.06, -0.73,  0.77, -0.67,
              -0.73,  2.64,  1.04,  0.72,
               0.77,  1.04,  3.93, -2.14,
              -0.67,  0.72, -2.14,  2.04), 4, 4, byrow = TRUE)
eigjacobi(A)
# $V
#            [,1]       [,2]       [,3]       [,4]
# [1,] 0.87019414 -0.3151209  0.1975473 -0.3231656
# [2,] 0.11138094  0.8661855  0.1178032 -0.4726938
# [3,] 0.07043799  0.1683401  0.8273261  0.5312548
# [4,] 0.47475776  0.3494040 -0.5124734  0.6244140
# 
# $D
# [1] 0.66335457 3.39813189 5.58753257 0.02098098

Einstein Functions

Description

Einstein functions.

Usage

einsteinF(d, x)

Arguments

Details

The Einstein functions are sometimes used for the Planck-Einstein oscillator in one degree of freedom.

The functions are defined as:

E1(x) = \frac{x^2 e^x}{(e^x - 1)^2}

E2(x) = \frac{x}{e^x - 1}

E3(x) = ln(1 - e^{-x})

E4(x) = \frac{x}{e^x - 1} - ln(1 - e^{-x})

E1 has an inflection point as x=2.34694130....

Value

Numeric/complex scalar or vector.

Examples

## Not run: 
x1 <- seq(-4, 4, length.out = 101)
y1 <- einsteinF(1, x1)
plot(x1, y1, type = "l", col = "red",
             xlab = "", ylab = "", main = "Einstein Function E1(x)")
grid()

y2 <- einsteinF(2, x1)
plot(x1, y2, type = "l", col = "red",
             xlab = "", ylab = "", main = "Einstein Function E2(x)")
grid()

x3 <- seq(0, 5, length.out = 101)
y3 <- einsteinF(3, x3)
plot(x3, y3, type = "l", col = "red",
             xlab = "", ylab = "", main = "Einstein Function E3(x)")
grid()

y4 <- einsteinF(4, x3)
plot(x3, y4, type = "l", col = "red",
             xlab = "", ylab = "", main = "Einstein Function E4(x)")
grid()
## End(Not run)

Elliptic and Jacobi Elliptic Integrals

Description

Complete elliptic integrals of the first and second kind, and Jacobi elliptic integrals.

Usage

ellipke(m, tol = .Machine$double.eps)

ellipj(u, m, tol = .Machine$double.eps)

Arguments

Details

ellipke computes the complete elliptic integrals to accuracy tol, based on the algebraic-geometric mean.

ellipj computes the Jacobi elliptic integrals sn, cn, and dn. For instance, sn is the inverse function for

u = \int_0^\phi dt / \sqrt{1 - m \sin^2 t}

with sn(u) = \sin(\phi).

Some definitions of the elliptic functions use the modules k instead of the parameter m. They are related by k^2=m=sin(a)^2 where a is the ‘modular angle’.

Value

ellipke returns list with two components, k the values for the first kind, e the values for the second kind.

ellipj returns a list with components the three Jacobi elliptic integrals sn, cn, and dn.

References

Abramowitz, M., and I. A. Stegun (1965). Handbook of Mathematical Functions. Dover Publications, New York.

See Also

elliptic::sn,cn,dn

Examples

x <- linspace(0, 1, 20)
ke <- ellipke(x)

## Not run: 
plot(x, ke$k, type = "l", col ="darkblue", ylim = c(0, 5),
     main = "Elliptic Integrals")
lines(x, ke$e, col = "darkgreen")
legend( 0.01, 4.5,
        legend = c("Elliptic integral of first kind",
                   "Elliptic integral of second kind"),
        col = c("darkblue", "darkgreen"), lty = 1)
grid()
## End(Not run)

## ellipse circumference with axes a, b
ellipse_cf <- function(a, b) {
    return(4*a*ellipke(1 - (b^2/a^2))$e)
}
print(ellipse_cf(1.0, 0.8), digits = 10)
# [1] 5.672333578

## Jacobi elliptic integrals
u <- c(0, 1, 2, 3, 4, 5)
m <- seq(0.0, 1.0, by = 0.2)
je <- ellipj(u, m)
# $sn       0.0000  0.8265  0.9851  0.7433  0.4771  0.9999
# $cn       1.0000  0.5630 -0.1720 -0.6690 -0.8789  0.0135
# $dn       1.0000  0.9292  0.7822  0.8176  0.9044  0.0135
je$sn^2 + je$cn^2       # 1 1 1 1 1 1
je$dn^2 + m * je$sn^2   # 1 1 1 1 1 1

Floating Point Relative Accuracy

Description

Distance from 1.0 to the next largest double-precision number.

Usage

eps(x = 1.0)

Arguments

Details

d=eps(x) is the positive distance from abs(x) to the next larger floating point number in double precision.

If x is an array, eps(x) will return eps(max(abs(x))).

Value

Returns a scalar.

Examples

for (i in -5:5) cat(eps(10^i), "\n")
# 1.694066e-21 
# 1.355253e-20 
# 2.168404e-19 
# 1.734723e-18 
# 1.387779e-17 
# 2.220446e-16 
# 1.776357e-15 
# 1.421085e-14 
# 1.136868e-13 
# 1.818989e-12 
# 1.455192e-11 

Error Functions and Inverses (Matlab Style)

Description

The error or Phi function is a variant of the cumulative normal (or Gaussian) distribution.

Usage

erf(x)
erfinv(y)
erfc(x)
erfcinv(y)
erfcx(x)

erfz(z)
erfi(z)

Arguments

Details

erf and erfinv are the error and inverse error functions.
erfc and erfcinv are the complementary error function and its inverse.
erfcx is the scaled complementary error function.
erfz is the complex, erfi the imaginary error function.

Value

Real or complex number(s), the value(s) of the function.

Note

For the complex error function we used Fortran code from the book S. Zhang & J. Jin “Computation of Special Functions” (Wiley, 1996).

Author(s)

First version by Hans W Borchers; vectorized version of erfz by Michael Lachmann.

See Also

pnorm

Examples

  x <- 1.0
  erf(x); 2*pnorm(sqrt(2)*x) - 1
# [1] 0.842700792949715
# [1] 0.842700792949715
  erfc(x); 1 - erf(x); 2*pnorm(-sqrt(2)*x)
# [1] 0.157299207050285
# [1] 0.157299207050285
# [1] 0.157299207050285
  erfz(x)
# [1] 0.842700792949715
  erfi(x)
# [1] 1.650425758797543

Plot Error Bars

Description

Draws symmetric error bars in x- and/or y-direction.

Usage

errorbar(x, y, xerr = NULL, yerr = NULL,
         bar.col = "red", bar.len = 0.01,
         grid = TRUE, with = TRUE, add = FALSE, ...)

Arguments

Details

errorbar plots y versus x with symmetric error bars, with a length determined by xerr resp. yerr in x- and/or y-direction. If xerr or yerr is NULL error bars in this direction will not be drawn.

A future version will allow to draw unsymmetric error bars by specifying upper and lower limits when xerr or yerr is a matrix of size (2 x length(x)).

Value

Generates a plot, no return value.

See Also

plotrix::plotCI, Hmisc::errbar

Examples

## Not run: 
x <- seq(0, 2*pi, length.out = 20)
y <- sin(x)
xe <- 0.1
ye <- 0.1 * y
errorbar(x, y, xe, ye, type = "l", with = FALSE)

cnt <- round(100*randn(20, 3))
y <- apply(cnt, 1, mean)
e <- apply(cnt, 1, sd)
errorbar(1:20, y, yerr = e, bar.col = "blue")

## End(Not run)

Dirichlet Eta Function

Description

Dirichlet's eta function valid in the entire complex plane.

Usage

eta(z)

Arguments

Details

Computes the eta function for complex arguments using a series expansion.

Accuracy is about 13 significant digits for abs(z)<100, drops off with higher absolute values.

Value

Returns a complex vector of function values.

Note

Copyright (c) 2001 Paul Godfrey for a Matlab version available on Mathwork's Matlab Central under BSD license.

References

Zhang, Sh., and J. Jin (1996). Computation of Special Functions. Wiley-Interscience, New York.

See Also

gammaz, zeta

Examples

z <- 0.5 + (1:5)*1i
eta(z)
z <- c(0, 0.5+1i, 1, 1i, 2+2i, -1, -2, -1-1i)
eta(z)

Euler-Heun ODE Solver

Description

Euler and Euler-Heun ODE solver.

Usage

euler_heun(f, a, b, y0, n = 100, improved = TRUE, ...)

Arguments

Details

euler_heun is an integration method for ordinary differential equations using the simple Euler resp. the improved Euler-Heun Method.

Value

List with components t for grid (or ‘time’) points, and y the vector of predicted values at those grid points.

References

Quarteroni, A., and F. Saleri (). Scientific Computing with MATLAB and Octave. Second Edition, Springer-Verlag, Berlin Heidelberg, 2006.

See Also

cranknic

Examples

##  Flame-up process
f <- function(x, y) y^2 - y^3
s1 <- cranknic(f, 0, 200, 0.01)
s2 <- euler_heun(f, 0, 200, 0.01)
## Not run: 
plot(s1$t, s1$y, type="l", col="blue")
lines(s2$t, s2$y, col="red")
grid()
## End(Not run)

Exponential and Logarithmic Integral

Description

The exponential integral functions E1 and Ei and the logarithmic integral Li.

The exponential integral is defined for x > 0 as

\int_x^\infty \frac{e^{-t}}{t} dt

and by analytic continuation in the complex plane. It can also be defined as the Cauchy principal value of the integral

\int_{-\infty}^x \frac{e^t}{t} dt

This is denoted as Ei(x) and the relationship between Ei and expint(x) for x real, x > 0 is as follows:

Ei(x) = - E1(-x) -i \pi

The logarithmic integral li(x) for real x, x > 0, is defined as

li(x) = \int_0^x \frac{dt}{log(t)}

and the Eulerian logarithmic integral as Li(x) = li(x) - li(2).

The integral Li approximates the prime number function \pi(n), i.e., the number of primes below or equal to n (see the examples).

Usage

expint(x)
expint_E1(x)

expint_Ei(x)
li(x)

Arguments

Details

For x in [-38, 2] we use a series expansion, otherwise a continued fraction, see the references below, chapter 5.

Value

Returns a vector of real or complex numbers, the vectorized exponential integral, resp. the logarithmic integral.

Note

The logarithmic integral li(10^i)-li(2) is an approximation of the number of primes below 10^i, i.e., Pi(10^i), see “?primes”.

References

Abramowitz, M., and I.A. Stegun (1965). Handbook of Mathematical Functions. Dover Publications, New York.

See Also

gsl::expint_E1,expint_Ei, primes

Examples

expint_E1(1:10)
#   0.2193839  0.0489005  0.0130484  0.0037794  0.0011483
#   0.0003601  0.0001155  0.0000377  0.0000124  0.0000042
expint_Ei(1:10)

## Not run: 
estimPi <- function(n) round(Re(li(n) - li(2))) # estimated number of primes
primesPi <- function(n) length(primes(n))       # true number of primes <= n
N <- 1e6
(estimPi(N) - primesPi(N)) / estimPi(N)         # deviation is 0.16 percent!
## End(Not run)

Matrix Exponential

Description

Computes the exponential of a matrix.

Usage

expm(A, np = 128)

logm(A)

Arguments

Details

For an analytic function f and a matrix A the expression f(A) can be computed by the Cauchy integral

f(A) = (2 \pi i)^{-1} \int_G (zI-A)^{-1} f(z) dz

where G is a closed contour around the eigenvalues of A.

Here this is achieved by taking G to be a circle and approximating the integral by the trapezoid rule.

logm is a fake at the moment as it computes the matrix logarithm through taking the logarithm of its eigenvalues; will be replaced by an approach using Pade interpolation.

Another more accurate and more reliable approach for computing these functions can be found in the R package ‘expm’.

Value

Matrix of the same size as A.

Note

This approach could be used for other analytic functions, but a point to consider is which branch to take (e.g., for the logm function).

Author(s)

Idea and Matlab code for a cubic root by Nick Trefethen in his “10 digits 1 page” project.

References

Moler, C., and Ch. Van Loan (2003). Nineteen Dubious Ways to Compute the Exponential of a Matrix, Twenty-Five Years Later. SIAM Review, Vol. 1, No. 1, pp. 1–46.

N. J. Higham (2008). Matrix Functions: Theory and Computation. SIAM Society for Industrial and Applied Mathematics.

See Also

expm::expm

Examples

##  The Ward test cases described in the help for expm::expm agree up to
##  10 digits with the values here and with results from Matlab's expm !
A <- matrix(c(-49, -64, 24, 31), 2, 2)
expm(A)
# -0.7357588 0.5518191
# -1.4715176 1.1036382

A1 <- matrix(c(10,  7,  8,  7,
                7,  5,  6,  5,
                8,  6, 10,  9,
                7,  5,  9, 10), nrow = 4, ncol = 4, byrow = TRUE)
expm(logm(A1))
logm(expm(A1))

##  System of linear differential equations: y' = M y  (y = c(y1, y2, y3))
M <- matrix(c(2,-1,1, 0,3,-1, 2,1,3), 3, 3, byrow=TRUE)
M
C1 <- 0.5; C2 <- 1.0; C3 <- 1.5
t  <- 2.0; Mt <- expm(t * M)
yt <- Mt 

Some Basic Matrices

Description

Create basic matrices.

Usage

eye(n, m = n)
ones(n, m = n)
zeros(n, m = n)

Arguments

Value

Matrix of size n x m. Defaults to a square matrix if m is missing.

No dropping of dimensions; if n = 1, still returns a matrix and not a vector.

See Also

Diag,

Examples

eye(3)
ones(3, 1)
zeros(1, 3)

Contour, Surface, and Mesh Plotter

Description

Easy-to-use contour and 3-D surface resp mesh plotter.

Usage

ezcontour(f, xlim = c(-pi,pi), ylim = c(-pi,pi), 
          n = 60, filled = FALSE, col = NULL)

ezsurf(f, xlim = c(-pi, pi), ylim = c(-pi, pi),
       n = 60, ...)

ezmesh(f, xlim = c(-pi,pi), ylim = c(-pi,pi), 
       n = 60, ...)

Arguments

Details

ezcontour generates a contour plot of the function f using contour (and image if filled=TRUE is chosen). If filled=TRUE is chosen, col should be a color scheme, the default is heat.colors(12).

ezsurf resp. ezmesh generates a surface/mesh plot of the function f using persp.

The function f needs not be vectorized in any form.

Value

Plots the function graph and invisibly returns NULL.

Note

Mimicks Matlab functions of the same names; Matlab's ezcontourf can be generated with filled=TRUE.

See Also

contour, image, persp

Examples

## Not run: 
f <- function(xy) {
    x <- xy[1]; y <- xy[2]
    3*(1-x)^2 * exp(-(x^2) - (y+1)^2) -
        10*(x/5 - x^3 - y^5) * exp(-x^2 - y^2) -
        1/3 * exp(-(x+1)^2 - y^2)
    }
ezcontour(f, col = "navy")
ezcontour(f, filled = TRUE)
ezmesh(f)
ezmesh(f, col="lightblue", theta = -15, phi = 30)
  
## End(Not run)

Easy Function Plot

Description

Easy function plot w/o the need to define x, y coordinates.

Usage

fplot(f, interval, ...)

ezplot( f, a, b, n = 101, col = "blue", add = FALSE,
        lty = 1, lwd = 1, marker = 0, pch = 1,
        grid = TRUE, gridcol = "gray", 
        fill = FALSE, fillcol = "lightgray",
        xlab = "x", ylab = "f (x)", main = "Function Plot", ...)

Arguments

Details

Calculates the x, y coordinates of points to be plotted and calls the plot function.

If fill is TRUE, also calls the polygon function with the x, y coordinates in appropriate order.

If the no. of markers is greater than 2, this number of markers will be added to the curve, with equal distances measured along the curve.

Value

Plots the function graph and invisibly returns NULL.

Note

fplot is almost an alias for ezplot as all ez... will be replaced by MATLAB with function names f... in 2017.

ezplot should mimick the Matlab function of the same name, has more functionality, misses the possibility to plot several functions.

See Also

curve

Examples

## Not run: 
fun <- function(x) x * cos(0.1*exp(x)) * sin(0.1*pi*exp(x))
ezplot(fun, 0, 5, n = 1001, fill = TRUE)
  
## End(Not run)

Easy Polar Plot

Description

Easy function plot w/o the need to define x, y coordinates.

Usage

ezpolar(fun, interv = c(0, 2*pi))

Arguments

Details

Calculates the x, y coordinates of points to be plotted and calls the polar function.

Value

Plots the function graph and invisibly returns NULL.

Note

Mimick the Matlab function of the same name.

See Also

ezplot

Examples

## Not run: 
fun <- function(x) 1 + cos(x)
ezpolar(fun)
  
## End(Not run)

Factorial Function

Description

Factorial for non-negative integers n <= 170.

Usage

fact(n)

factorial2(n)

Arguments

Details

The factorial is computed by brute force; factorials for n >= 171 are not representable as ‘double’ anymore.

Value

fact returns the factorial of each element in n. If n < 0 the value is NaN, and for n > 170 it is Inf. Non-integers will be reduced to integers through floor(n).

factorial2 returns the product of all even resp. odd integers, depending on whether n is even or odd.

Note

The R core function factorial uses the gamma function, whose implementation is not accurate enough for larger input values.

See Also

factorial

Examples

fact(c(-1, 0, 1, NA, 171))  #=> NaN   1   1  NA Inf
fact(100)                   #=> 9.332621544394410e+157
factorial(100)              #=> 9.332621544394225e+157
# correct value:                9.332621544394415e+157
# Stirling's approximation:     9.324847625269420e+157
# n! ~ sqrt(2*pi*n) * (n/e)^n

factorial2(8);  factorial2(9);  factorial2(10)  # 384   945  3840
factorial(10) / factorial2(10)                  # => factorial2(9)

Prime Factors

Description

Returns a vector containing the prime factors of n.

Usage

  factors(n)

Arguments

Details

Computes the prime factors of n in ascending order, each one as often as its multiplicity requires, such that n == prod(factors(n)).

The corresponding Matlab function is called ‘factor’, but because factors have a special meaning in R and the factor() function in R could not (or should not) be shadowed, the number theoretic function has been renamed here.

Value

Vector containing the prime factors of n.

See Also

isprime, primes

Examples

## Not run: 
  factors(1002001)       # 7  7  11  11  13  13
  factors(65537)         # is prime
  # Euler's calculation
  factors(2^32 + 1)      # 641  6700417
## End(Not run)

Numerical Differentiation

Description

Numerical function differentiation for orders n=1..4 using finite difference approximations.

Usage

fderiv(f, x, n = 1, h = 0,
        method = c("central", "forward", "backward"), ...)

Arguments

Details

Derivatives are computed applying central difference formulas that stem from the Taylor series approximation. These formulas have a convergence rate of O(h^2).

Use the ‘forward’ (right side) or ‘backward’ (left side) method if the function can only be computed or is only defined on one side. Otherwise, always use the central difference formulas.

Optimal step sizes depend on the accuracy the function can be computed with. Assuming internal functions with an accuracy 2.2e-16, appropriate step sizes might be 5e-6, 1e-4, 5e-4, 2.5e-3 for n=1,...,4 and precisions of about 10^-10, 10^-8, 5*10^-7, 5*10^-6 (at best).

For n>4 a recursion (or finite difference) formula will be applied, cd. the Wikipedia article on “finite difference”.

Value

Vector of the same length as x.

Note

Numerical differentiation suffers from the conflict between round-off and truncation errors.

References

Kiusalaas, J. (2005). Numerical Methods in Engineering with Matlab. Cambridge University Press.

See Also

numderiv, taylor

Examples

## Not run: 
f <- sin
xs <- seq(-pi, pi, length.out = 100)
ys <- f(xs)
y1 <- fderiv(f, xs, n = 1, method = "backward")
y2 <- fderiv(f, xs, n = 2, method = "backward")
y3 <- fderiv(f, xs, n = 3, method = "backward")
y4 <- fderiv(f, xs, n = 4, method = "backward")
plot(xs, ys, type = "l", col = "gray", lwd = 2,
     xlab = "", ylab = "", main = "Sinus and its Derivatives")
lines(xs, y1, col=1, lty=2)
lines(xs, y2, col=2, lty=3)
lines(xs, y3, col=3, lty=4)
lines(xs, y4, col=4, lty=5)
grid()
## End(Not run)

Fibonacci Search

Description

Fibonacci search for function minimum.

Usage

fibsearch(f, a, b, ..., endp = FALSE, tol = .Machine$double.eps^(1/2))

Arguments

Details

Fibonacci search for a univariate function minimum in a bounded interval.

Value

Return a list with components xmin, fmin, the function value at the minimum, niter, the number of iterations done, and the estimated precision estim.prec

See Also

uniroot

Examples

f <- function(x) x * cos(0.1*exp(x)) * sin(0.1*pi*exp(x))
fibsearch(f, 0, 4, tol=10^-10)   # $xmin    = 3.24848329403424
optimize(f, c(0,4), tol=10^-10)  # $minimum = 3.24848328971188

Control Plot Devices (Matlab Style)

Description

Open, activate, and close grahics devices.

Usage

figure(figno, title = "")

Arguments

Details

The number of a graphics device cannot be 0 or 1. The function will work for the operating systems Mac OS, MS Windows, and most Linux systems.

If figno is negative and a graphics device with that number does exist, it will be closed.

Value

No return value, except when a device of that number does not exist, in which case it returns a list of numbers of open graphics devices.

Note

Does not bring the activated graphics device in front.

See Also

dev.set, dev.off, dev.list

Examples

## Not run: 
figure()
figure(-2)

## End(Not run)

Find Interval Indices

Description

Find indices i in vector xs such that either x=xs[i] or such that xs[i]<x<xs[i+1] or xs[i]>x>xs[i+1].

Usage

findintervals(x, xs)

Arguments

Details

Contrary to findInterval, the vector xs in findintervals need not be sorted.

Value

Vector of indices in 1..length(xs). If none is found, returns integer(0).

Note

If x is equal to the last element in xs, the index length(xs) will also be returned.

Examples

xs <- zapsmall(sin(seq(0, 10*pi, len=100)))
findintervals(0, xs)
#   1  10  20  30  40  50  60  70  80  90 100

Find All Minima

Description

Finding all local(!) minima of a unvariate function in an interval by splitting the interval in many small subintervals.

Usage

findmins(f, a, b, n = 100, tol = .Machine$double.eps^(2/3), ...)

Arguments

Details

Local minima are found by looking for one minimum in each subinterval. It will be found by applying optimize to any two adjacent subinterval where the first slope is negative and the second one positive.

If the function is minimal on a whole subinterval, this will cause problems. If some minima are apparently not found, increase the number of subintervals.

Note that the endpoints of the interval will never be considered to be local minima. The function need not be vectorized.

Value

Numeric vector with the x-positions of all minima found in the interval.

See Also

optimize

Examples

fun <- function(x) x * cos(0.1*exp(x)) * sin(0.1*pi*exp(x))
## Not run: ezplot(fun, 0, 5, n = 1001)

# If n is smaller, the rightmost minimum will not be found.
findmins(fun, 0, 5, n= 1000)
#  2.537727 3.248481 3.761840 4.023021 4.295831
#  4.455115 4.641481 4.756263 4.897461 4.987802

Find Peaks

Description

Find peaks (maxima) in a time series.

Usage

findpeaks(x, nups = 1, ndowns = nups, zero = "0", peakpat = NULL,
          minpeakheight = -Inf, minpeakdistance = 1,
          threshold = 0, npeaks = 0, sortstr = FALSE)

Arguments

Details

This function is quite general as it relies on regular patterns to determine where a peak is located, from beginning to end.

Value

Returns a matrix where each row represents one peak found. The first column gives the height, the second the position/index where the maximum is reached, the third and forth the indices of where the peak begins and ends — in the sense of where the pattern starts and ends.

Note

On Matlab Central there are several realizations for finding peaks, for example “peakfinder”, “peakseek”, or “peakdetect”. And “findpeaks” is also the name of a function in the Matlab ‘signal’ toolbox.

The parameter names are taken from the “findpeaks” function in ‘signal’, but the implementation utilizing regular expressions is unique and fast.

See Also

hampel

Examples

x <- seq(0, 1, len = 1024)
pos <- c(0.1, 0.13, 0.15, 0.23, 0.25, 0.40, 0.44, 0.65, 0.76, 0.78, 0.81)
hgt <- c(4, 5, 3, 4, 5, 4.2, 2.1, 4.3, 3.1, 5.1, 4.2)
wdt <- c(0.005, 0.005, 0.006, 0.01, 0.01, 0.03, 0.01, 0.01, 0.005, 0.008, 0.005)

pSignal <- numeric(length(x))
for (i in seq(along=pos)) {
	pSignal <- pSignal + hgt[i]/(1 + abs((x - pos[i])/wdt[i]))^4
}
findpeaks(pSignal, npeaks=3, threshold=4, sortstr=TRUE)

## Not run: 
plot(pSignal, type="l", col="navy")
grid()
x <- findpeaks(pSignal, npeaks=3, threshold=4, sortstr=TRUE)
points(x[, 2], x[, 1], pch=20, col="maroon")
## End(Not run)

find function (Matlab Style)

Description

Finds indices of nonzero elements.

Usage

finds(v)

Arguments

Details

Finds indices of true or nonzero elements of argument v; can be used with a logical expression.

Value

Indices of elements matching the expression x.

Examples

finds(-3:3 >= 0)
finds(c(0, 1, 0, 2, 3))

Find All Roots

Description

Finding all roots of a unvariate function in an interval by splitting the interval in many small subintervals.

Usage

findzeros(f, a, b, n = 100, tol = .Machine$double.eps^(2/3), ...)

Arguments

Details

Roots, i.e. zeros in a subinterval will be found by applying uniroot to any subinterval where the sign of the function changes. The endpoints of the interval will be tested separately.

If the function points are both positive or negative and the slope in this interval is high enough, the minimum or maximum will be determined with optimize and checked for a possible zero.

The function need not be vectorized.

Value

Numeric vector with the x-positions of all roots found in the interval.

See Also

findmins

Examples

f1 <- function(x) sin(pi/x)
findzeros(f1, 1/10, 1)
#  0.1000000  0.1111028  0.1250183  0.1428641  0.1666655
#  0.2000004  0.2499867  0.3333441  0.4999794  1.0000000

f2 <- function(x) 0.5*(1 + sin(10*pi*x))
findzeros(f2, 0, 1)
#  0.15  0.35  0.55  0.75  0.95

f3 <- function(x) sin(pi/x) + 1
findzeros(f3, 0.1, 0.5)
# 0.1052632 0.1333333 0.1818182 0.2857143

f4 <- function(x) sin(pi/x) - 1
findzeros(f4, 0.1, 0.5)
# 0.1176471 0.1538462 0.2222222 0.4000000

## Not run: 
# Dini function
Dini <- function(x) x * besselJ(x, 1) + 3 * besselJ(x, 0)
findzeros(Dini, 0, 100, n = 128)
ezplot(Dini, 0, 100, n = 512)

## End(Not run)

Fletcher-Powell Conjugate Gradient Minimization

Description

Conjugate Gradient (CG) minimization through the Davidon-Fletcher-Powell approach for function minimization.

The Davidon-Fletcher-Powell (DFP) and the Broyden-Fletcher-Goldfarb-Shanno (BFGS) methods are the first quasi-Newton minimization methods developed. These methods differ only in some details; in general, the BFGS approach is more robust.

Usage

fletcher_powell(x0, f, g = NULL,
                maxiter = 1000, tol = .Machine$double.eps^(2/3))

Arguments

Details

The starting point is Newton's method in the multivariate case, when the estimate of the minimum is updated by the following equation

x_{new} = x - H^{-1}(x) grad(g)(x)

where H is the Hessian and grad the gradient.

The basic idea is to generate a sequence of good approximations to the inverse Hessian matrix, in such a way that the approximations are again positive definite.

Value

List with following components:

Note

Used some Matlab code as described in the book “Applied Numerical Analysis Using Matlab” by L. V.Fausett.

References

J. F. Bonnans, J. C. Gilbert, C. Lemarechal, and C. A. Sagastizabal. Numerical Optimization: Theoretical and Practical Aspects. Second Edition, Springer-Verlag, Berlin Heidelberg, 2006.

See Also

steep_descent

Examples

##  Rosenbrock function
rosenbrock <- function(x) {
    n <- length(x)
    x1 <- x[2:n]
    x2 <- x[1:(n-1)]
    sum(100*(x1-x2^2)^2 + (1-x2)^2)
}
fletcher_powell(c(0, 0), rosenbrock)
# $xmin
# [1] 1 1
# $fmin
# [1] 1.774148e-27
# $niter
# [1] 14

Matrix Flipping (Matlab Style)

Description

Flip matrices up and down or left and right; or circulating indices per dimension.

Usage

flipdim(a, dim)
flipud(a)
fliplr(a)
circshift(a, sz)

Arguments

Details

flipdim will flip a matrix along the dim dimension, where dim=1 means flipping rows, and dim=2 flipping the columns.

flipud and fliplr are simply shortcuts for flipdim(a, 1) resp. flipdim(a, 2).

circshift(a, sz) circulates each dimension (should be applicable to arrays).

Value

the original matrix somehow flipped or circularly shifted.

Examples

a <- matrix(1:12, nrow=3, ncol=4, byrow=TRUE)
flipud(a)
fliplr(a)

circshift(a, c(1, -1))
v <- 1:10
circshift(v, 5)

Finding Function Minimum

Description

Find minimum of single-variable function on fixed interval.

Usage

fminbnd(f, a, b, maxiter = 1000, maximum = FALSE,
        tol = 1e-07, rel.tol = tol, abs.tol = 1e-15, ...)

Arguments

Details

fminbnd finds the minimum of a function of one variable within a fixed interval. It applies Brent's algorithm, based on golden section search and parabolic interpolation.

fminbnd may only give local solutions. fminbnd never evaluates f at the endpoints.

Value

List with

Note

fminbnd mimics the Matlab function of the same name.

References

R. P. Brent (1973). Algorithms for Minimization Without Derivatives. Dover Publications, reprinted 2002.

See Also

fibsearch, golden_ratio

Examples

##  CHEBFUN example by Trefethen
f <- function(x) exp(x)*sin(3*x)*tanh(5*cos(30*x))
fminbnd(f, -1, 1)                   # fourth local minimum (from left)
g <- function(x) complexstep(f, x)  # complex-step derivative
xs <- findzeros(g, -1, 1)           # local minima and maxima
ys <- f(xs); n0 <- which.min(ys)    # index of global minimum
fminbnd(f, xs[n0-1], xs[n0+1])      # xmin:0.7036632, fmin: -1.727377

## Not run: 
ezplot(f, -1, 1, n = 1000, col = "darkblue", lwd = 2)
ezplot(function(x) g(x)/150, -1, 1, n = 1000, col = "darkred", add = TRUE)
grid()
## End(Not run)

Minimize Nonlinear Constrained Multivariable Function.

Description

Find minimum of multivariable functions with nonlinear constraints.

Usage

  fmincon(x0, fn, gr = NULL, ..., method = "SQP",
          A = NULL, b = NULL, Aeq = NULL, beq = NULL,
          lb = NULL, ub = NULL, hin = NULL, heq = NULL,
          tol = 1e-06, maxfeval = 10000, maxiter = 5000)

Arguments

Details

Wraps the function solnl in the 'NlcOptim' package. The underlying method is a Squential Quadratic Programming (SQP) approach.

Constraints can be defined in different ways, as linear constraints in matrix form, as nonlinear functions, or as bounds constraints.

Value

List with the following components:

Note

fmincon mimics the Matlab function of the same name.

Author(s)

Xianyan Chen for the package NlcOptim.

References

J. Nocedal and S. J. Wright (2006). Numerical Optimization. Second Edition, Springer Science+Business Media, New York.

See Also

fminsearch, fminunc,

Examples

# Classical Rosenbrock function
n <- 10; x0 <- rep(1/n, n)
fn <- function(x) {n <- length(x)
    x1 <- x[2:n]; x2 <- x[1:(n - 1)]
    sum(100 * (x1 - x2^2)^2 + (1 - x2)^2)
}
# Equality and inequality constraints
heq1 <- function(x) sum(x)-1.0
hin1 <- function(x) -1 * x
hin2 <- function(x) x - 0.5
ub <- rep(0.5, n)

# Apply constraint minimization
res <- fmincon(x0, fn, hin = hin1, heq = heq1)
res$par; res$value

Derivative-free Nonlinear Function Minimization

Description

Find minimum of multivariable functions using derivative-free methods.

Usage

fminsearch(fn, x0, ..., lower = NULL, upper = NULL,
           method = c("Nelder-Mead", "Hooke-Jeeves"),
           minimize = TRUE, maxiter = 1000, tol = 1e-08)

Arguments

Details

fminsearch finds the minimum of a nonlinear scalar multivariable function, starting at an initial estimate and returning a value x that is a local minimizer of the function. With minimize=FALSE it searches for a maximum, by default for a (local) minimum.

As methods/solvers "Nelder-Mead" and "Hooke-Jeeves" are available. Only Hooke-Jeeves can handle bounds constraints. For nonlinear constraints see fmincon, and for methods using gradients see fminunc.

Important: fminsearch may only give local solutions.

Value

List with

Note

fminsearch mimics the Matlab function of the same name.

References

Nocedal, J., and S. Wright (2006). Numerical Optimization. Second Edition, Springer-Verlag, New York.

See Also

nelder_mead, hooke_jeeves

Examples

# Rosenbrock function
rosena <- function(x, a) 100*(x[2]-x[1]^2)^2 + (a-x[1])^2  # min: (a, a^2)

fminsearch(rosena, c(-1.2, 1), a = sqrt(2), method="Nelder-Mead")
## $xmin                   $fmin
## [1] 1.414292 2.000231   [1] 1.478036e-08

fminsearch(rosena, c(-1.2, 1), a = sqrt(2), method="Hooke-Jeeves")
## $xmin                   $fmin
## [1] 1.414215 2.000004   [1] 1.79078e-12

Minimize Unconstrained Multivariable Function

Description

Find minimum of unconstrained multivariable functions.

Usage

  fminunc(x0, fn, gr = NULL, ...,
          tol = 1e-08, maxiter = 0, maxfeval = 0)

Arguments

Details

The method used here for unconstrained minimization is a variant of a "variable metric" resp. quasi-Newton approach.

Value

List with the following components:

Note

fminunc mimics the Matlab function of the same name.

Author(s)

The "variable metric" code provided by John Nash (package Rvmmin), stripped-down version by Hans W. Borchers.

References

J. Nocedal and S. J. Wright (2006). Numerical Optimization. Second Edition, Springer Science+Business Media, New York.

See Also

fminsearch, fmincon,

Examples

  fun = function(x) 
          x[1]*exp(-(x[1]^2 + x[2]^2)) + (x[1]^2 + x[2]^2)/20
  fminunc(x0 = c(1, 2), fun)
  ## xmin: c(-0.6691, 0.0000); fmin: -0.4052

Function Norm

Description

The fnorm function calculates several different types of function norms for depending on the argument p.

Usage

fnorm(f, g, x1, x2, p = 2, npoints = 100)

Arguments

Details

fnorm returns a scalar that gives some measure of the distance of two functions f and g on the interval [x1, x2].

It takes npoints equidistant points in the interval, computes the function values for f and g and applies Norm to their difference.

Especially p=Inf returns the maximum norm, while fnorm(f, g, x1, x2, p = 1, npoints) / npoints would return some estimate of the mean distance.

Value

Numeric scalar (or Inf), or NA if one of these functions returns NA.

Note

Another kind of ‘mean’ distance could be calculated by integrating the difference f-g and dividing through the length of the interval.

See Also

Norm

Examples

xp <- seq(-1, 1, length.out = 6)
yp <- runge(xp)
p5 <- polyfit(xp, yp, 5)
f5 <- function(x) polyval(p5, x)
fnorm(runge, f5, -1, 1, p = Inf)                  #=> 0.4303246
fnorm(runge, f5, -1, 1, p = Inf, npoints = 1000)  #=> 0.4326690

# Compute mean distance using fnorm:
fnorm(runge, f5, -1, 1, p = 1, 1000) / 1000       #=> 0.1094193

# Compute mean distance by integration:
fn <- function(x) abs(runge(x) - f5(x))
integrate(fn, -1, 1)$value / 2                    #=> 0.1095285

Fornberg's Finite Difference Approximation

Description

Finite difference approximation using Fornberg's method for the derivatives of order 1 to k based on irregulat grid values.

Usage

fornberg(x, y, xs, k = 1)

Arguments

Details

Compute coefficients for finite difference approximation for the derivative of order k at xs based on grid values at points in x. For k=0 this will evaluate the interpolating polynomial itself, but call it with k=1.

Value

Returns a matrix of size (length(xs)), where the (k+1)-th column gives the value of the k-th derivative. Especially the first column returns the polynomial interpolation of the function.

Note

Fornberg's method is considered to be numerically more stable than applying Vandermonde's matrix.

References

LeVeque, R. J. (2007). Finite Difference Methods for Ordinary and Partial Differential Equations. Society for Industrial and Applied Mathematics (SIAM), Philadelphia.

See Also

neville, newtonInterp

Examples

x <- 2 * pi * c(0.0, 0.07, 0.13, 0.2, 0.28, 0.34, 0.47, 0.5, 0.71, 0.95, 1.0)
y <- sin(0.9*x)
xs <- linspace(0, 2*pi, 51)
fornb <- fornberg(x, y, xs, 10)
## Not run: 
matplot(xs, fornb, type="l")
grid()
## End(Not run)

Formatted Printing (Matlab style)

Description

Formatted printing to stdout or a file.

Usage

fprintf(fmt, ..., file = "", append = FALSE)

Arguments

Details

fprintf applies the format string fmt to all input data ... and writes the result to standard output or a file. The usual C-style string formatting commands are used-

Value

Returns invisibly the number of bytes printed (using nchar).

See Also

sprintf

Examples

##  Examples:
nbytes <- fprintf("Results are:\n", file = "")
for (i in 1:10) {
    fprintf("%4d  %15.7f\n", i, exp(i), file = "")
}

Fractal Curves

Description

Generates the following fractal curves: Dragon Curve, Gosper Flowsnake Curve, Hexagon Molecule Curve, Hilbert Curve, Koch Snowflake Curve, Sierpinski Arrowhead Curve, Sierpinski (Cross) Curve, Sierpinski Triangle Curve.

Usage

fractalcurve(n, which = c("hilbert", "sierpinski", "snowflake",
    "dragon", "triangle", "arrowhead", "flowsnake", "molecule"))

Arguments

Details

The Hilbert curve is a continuous curve in the plane with 4^N points.

The Sierpinski (cross) curve is a closed curve in the plane with 4^(N+1)+1 points.

His arrowhead curve is a continuous curve in the plane with 3^N+1 points, and his triangle curve is a closed curve in the plane with 2*3^N+2 points.

The Koch snowflake curve is a closed curve in the plane with 3*2^N+1 points.

The dragon curve is a continuous curve in the plane with 2^(N+1) points.

The flowsnake curve is a continuous curve in the plane with 7^N+1 points.

The hexagon molecule curve is a closed curve in the plane with 6*3^N+1 points.

Value

Returns a list with x, y the x- resp. y-coordinates of the generated points describing the fractal curve.

Author(s)

Copyright (c) 2011 Jonas Lundgren for the Matlab toolbox fractal curves available on MatlabCentral under BSD license; here re-implemented in R with explicit allowance from the author.

References

Peitgen, H.O., H. Juergens, and D. Saupe (1993). Fractals for the Classroom. Springer-Verlag Berlin Heidelberg.

Examples

## The Hilbert curve transforms a 2-dim. function into a time series.
z <- fractalcurve(4, which = "hilbert")

## Not run: 
f1 <- function(x, y) x^2 + y^2
plot(f1(z$x, z$y), type = 'l', col = "darkblue", lwd = 2,
     ylim = c(-1, 2), main = "Functions transformed by Hilbert curves")

f2 <- function(x, y) x^2 - y^2
lines(f2(z$x, z$y), col = "darkgreen", lwd = 2)

f3 <- function(x, y) x^2 * y^2
lines(f3(z$x, z$y), col = "darkred", lwd = 2)
grid()
## End(Not run)

## Not run: 
## Show some more fractal surves
n <- 8
opar <- par(mfrow=c(2,2), mar=c(2,2,1,1))

z <- fractalcurve(n, which="dragon")
x <- z$x; y <- z$y
plot(x, y, type='l', col="darkgrey", lwd=2)
title("Dragon Curve")

z <- fractalcurve(n, which="molecule")
x <- z$x; y <- z$y
plot(x, y, type='l', col="darkblue")
title("Molecule Curve")

z <- fractalcurve(n, which="arrowhead")
x <- z$x; y <- z$y
plot(x, y, type='l', col="darkgreen")
title("Arrowhead Curve")

z <- fractalcurve(n, which="snowflake")
x <- z$x; y <- z$y
plot(x, y, type='l', col="darkred", lwd=2)
title("Snowflake Curve")

par(opar)
## End(Not run)

Fresnel Integrals

Description

(Normalized) Fresnel integrals S(x) and C(x)

Usage

fresnelS(x)
fresnelC(x)

Arguments

Details

The normalized Fresnel integrals are defined as

S(x) = \int_0^x \sin(\pi/2 \, t^2) dt

C(x) = \int_0^x \cos(\pi/2 \, t^2) dt

This program computes the Fresnel integrals S(x) and C(x) using Fortran code by Zhang and Jin. The accuracy is almost up to Machine precision.

The functions are not (yet) truly vectorized, but use a call to ‘apply’. The underlying function .fresnel (not exported) computes single values of S(x) and C(x) at the same time.

Value

Numeric vector of function values.

Note

Copyright (c) 1996 Zhang and Jin for the Fortran routines, converted to Matlab using the open source project ‘f2matlab’ by Ben Barrowes, posted to MatlabCentral in 2004, and then translated to R by Hans W. Borchers.

References

Zhang, S., and J. Jin (1996). Computation of Special Functions. Wiley-Interscience.

See Also

gaussLegendre

Examples

##  Compute Fresnel integrals through Gauss-Legendre quadrature
f1 <- function(t) sin(0.5 * pi * t^2)
f2 <- function(t) cos(0.5 * pi * t^2)
for (x in seq(0.5, 2.5, by = 0.5)) {
    cgl <- gaussLegendre(51, 0, x)
    fs <- sum(cgl$w * f1(cgl$x))
    fc <- sum(cgl$w * f2(cgl$x))
    cat(formatC(c(x, fresnelS(x), fs, fresnelC(x), fc),
        digits = 8, width = 12, flag = " ----"), "\n")
}

## Not run: 
xs <- seq(0, 7.5, by = 0.025)
ys <- fresnelS(xs)
yc <- fresnelC(xs)

##  Function plot of the Fresnel integrals
plot(xs, ys, type = "l", col = "darkgreen",
    xlim = c(0, 8), ylim = c(0, 1),
    xlab = "", ylab = "", main = "Fresnel Integrals")
lines(xs, yc, col = "blue")
legend(6.25, 0.95, c("S(x)", "C(x)"), col = c("darkgreen", "blue"), lty = 1)
grid()

##  The Cornu (or Euler) spiral
plot(c(-1, 1), c(-1, 1), type = "n",
    xlab = "", ylab = "", main = "Cornu Spiral")
lines(ys, yc, col = "red")
lines(-ys, -yc, col = "red")
grid()
## End(Not run)

Solve System of Nonlinear Equations

Description

Solve a system of m nonlinear equations of n variables.

Usage

fsolve(f, x0, J = NULL,
       maxiter = 100, tol = .Machine$double.eps^(0.5), ...)

Arguments

Details

fsolve tries to solve the components of function f simultaneously and uses the Gauss-Newton method with numerical gradient and Jacobian. If m = n, it uses broyden. Not applicable for univariate root finding.

Value

List with

Note

fsolve mimics the Matlab function of the same name.

References

Antoniou, A., and W.-S. Lu (2007). Practical Optimization: Algorithms and Engineering Applications. Springer Science+Business Media, New York.

See Also

broyden, gaussNewton

Examples

## Not run: 
# Find a matrix X such that X * X * X = [1, 2; 3, 4]
  F <- function(x) {
    a <- matrix(c(1, 3, 2, 4), nrow = 2, ncol = 2, byrow = TRUE)
    X <- matrix(x,             nrow = 2, ncol = 2, byrow = TRUE)
    return(c(X %*% X %*% X - a))
  }
  x0 <- matrix(1, 2, 2)
  X  <- matrix(fsolve(F, x0)$x, 2, 2)
  X
  # -0.1291489  0.8602157
  #  1.2903236  1.1611747

## End(Not run)

Root Finding Algorithm

Description

Find root of continuous function of one variable.

Usage

fzero(fun, x, maxiter = 500, tol = 1e-12, ...)

Arguments

Details

fzero tries to find a zero of f near x, if x is a scalar. Expands the interval until different signs are found at the endpoints or the maximum number of iterations is exceeded. If x is a vector of length two, fzero assumes x is an interval where the sign of x[1] differs from the sign of x[1]. An error occurs if this is not the case.

“This is essentially the ACM algorithm 748. The structure of the algorithm has been transformed non-trivially: it implement here a FSM version using one interior point determination and one bracketing per iteration, thus reducing the number of temporary variables and simplifying the structure.”

This approach will not find zeroes of quadratic order.

Value

fzero returns a list with

Note

fzero mimics the Matlab function of the same name, but is translated from Octave's fzero function, copyrighted (c) 2009 by Jaroslav Hajek.

References

Alefeld, Potra and Shi (1995). Enclosing Zeros of Continuous Functions. ACM Transactions on Mathematical Software, Vol. 21, No. 3.

See Also

uniroot, brent

Examples

fzero(sin, 3)                    # 3.141593
fzero(cos,c(1, 2))               # 1.570796
fzero(function(x) x^3-2*x-5, 2)  # 2.094551

Complex Root Finding

Description

Find the root of a complex function

Usage

fzsolve(fz, z0)

Arguments

Details

fzsolve tries to find the root of the complex and relatively smooth (i.e., analytic) function near a starting point.

The function is considered as real function R^2 --> R^2 and the newtonsys function is applied.

Value

Complex point with sufficiently small function value.

See Also

newtonsys

Examples

fz <- function(z) sin(z)^2 + sqrt(z) - log(z)
fzsolve(fz, 1+1i)
# 0.2555197+0.8948303i

Incomplete Gamma Function

Description

Lower and upper incomplete gamma function.

Usage

gammainc(x, a)

incgam(x, a)

Arguments

Details

gammainc computes the lower and upper incomplete gamma function, including the regularized gamma function. The lower and upper incomplete gamma functions are defined as

\gamma(x, a) = \int_0^x e^{-t} \, t^{a-1} \, dt

and

\Gamma(x, a) = \int_x^{\infty} e^{-t} \, t^{a-1} \, dt

while the regularized incomplete gamma function is \gamma(x, a)/\Gamma(a).

incgam (a name used in Pari/GP) computes the upper incomplete gamma function alone, applying the R function pgamma. The accuracy is thus much higher. It works for a >= -1, for even smaller values a recursion will give the result.

Value

gammainc returns a list with the values of the lower, the upper, and regularized lower incomplete gamma function. incgam only returns the value of the incomplete upper gamma function.

Note

Directly converting Fortran code is often easier than translating Matlab code generated with f2matlab.

References

Zhang, Sh., and J. Jin (1996). Computation of Special Functions. Wiley-Interscience, New York.

See Also

gamma, pgamma

Examples

gammainc( 1.5, 2)
gammainc(-1.5, 2)

incgam(3, 1.2)
incgam(3, 0.5); incgam(3, -0.5)

Complex Gamma Function

Description

Gamma function valid in the entire complex plane.

Usage

gammaz(z)

Arguments

Details

Computes the Gamma function for complex arguments using the Lanczos series approximation.

Accuracy is 15 significant digits along the real axis and 13 significant digits elsewhere.

To compute the logarithmic Gamma function use log(gammaz(z)).

Value

Returns a complex vector of function values.

Note

Copyright (c) 2001 Paul Godfrey for a Matlab version available on Mathwork's Matlab Central under BSD license.

Numerical Recipes used a 7 terms formula for a less effective approximation.

References

Zhang, Sh., and J. Jin (1996). Computation of Special Functions. Wiley-Interscience, New York.

See Also

gamma, gsl::lngamma_complex

Examples

max(gamma(1:10) - gammaz(1:10))
gammaz(-1)
gammaz(c(-2-2i, -1-1i, 0, 1+1i, 2+2i))

# Euler's reflection formula
z <- 1+1i
gammaz(1-z) * gammaz(z)  # == pi/sin(pi*z)

Gauss-Kronrod Quadrature

Description

Simple Gaussian-Kronrod quadrature formula.

Usage

gauss_kronrod(f, a, b, ...)

Arguments

Details

Gaussian quadrature of degree 7 with Gauss-Kronrod of degree 15 for error estimation, the quadQK15 procedure in the QUADPACK library.

Value

List of value and relative error.

Note

The function needs to be vectorized (though this could easily be changed), but the function does not need to be defined at the end points.

References

Fausett, L. V. (2007). Applied Numerical Analysis Using Matlab. Second edition, Prentice Hall.

See Also

quadgk, romberg

Examples

gauss_kronrod(sin, 0, pi)  #  2.000000000000000 , rel.error: 1.14e-12
gauss_kronrod(exp, 0, 1)   #  1.718281828459045 , rel.error: 0
                           #  1.718281828459045 , i.e. exp(1) - 1

Gauss-Hermite Quadrature Formula

Description

Nodes and weights for the n-point Gauss-Hermite quadrature formula.

Usage

gaussHermite(n)

Arguments

Details

Gauss-Hermite quadrature is used for integrating functions of the form

\int_{-\infty}^{\infty} f(x) e^{-x^2} dx

over the infinite interval ]-\infty, \infty[.

x and w are obtained from a tridiagonal eigenvalue problem. The value of such an integral is then sum(w*f(x)).

Value

List with components x, the nodes or points in]-Inf, Inf[, and w, the weights applied at these nodes.

Note

The basic quadrature rules are well known and can, e. g., be found in Gautschi (2004) — and explicit Matlab realizations in Trefethen (2000). These procedures have also been implemented in Matlab by Geert Van Damme, see his entries at MatlabCentral since 2010.

References

Gautschi, W. (2004). Orthogonal Polynomials: Computation and Approximation. Oxford University Press.

Trefethen, L. N. (2000). Spectral Methods in Matlab. SIAM, Society for Industrial and Applied Mathematics.

See Also

gaussLegendre, gaussLaguerre

Examples

cc <- gaussHermite(17)
# Integrate  exp(-x^2)  from -Inf to Inf
sum(cc$w)                        #=> 1.77245385090552 == sqrt(pi)
# Integrate  x^2 exp(-x^2)
sum(cc$w * cc$x^2)               #=> 0.88622692545276 == sqrt(pi) /2
# Integrate  cos(x) * exp(-x^2)
sum(cc$w * cos(cc$x))            #=> 1.38038844704314 == sqrt(pi)/exp(1)^0.25

Gauss-Laguerre Quadrature Formula

Description

Nodes and weights for the n-point Gauss-Laguerre quadrature formula.

Usage

gaussLaguerre(n, a = 0)

Arguments

Details

Gauss-Laguerre quadrature is used for integrating functions of the form

\int_0^{\infty} f(x) x^a e^{-x} dx

over the infinite interval ]0, \infty[.

x and w are obtained from a tridiagonal eigenvalue problem. The value of such an integral is then sum(w*f(x)).

Value

List with components x, the nodes or points in[0, Inf[, and w, the weights applied at these nodes.

Note

The basic quadrature rules are well known and can, e. g., be found in Gautschi (2004) — and explicit Matlab realizations in Trefethen (2000). These procedures have also been implemented in Matlab by Geert Van Damme, see his entries at MatlabCentral since 2010.

References

Gautschi, W. (2004). Orthogonal Polynomials: Computation and Approximation. Oxford University Press.

Trefethen, L. N. (2000). Spectral Methods in Matlab. SIAM, Society for Industrial and Applied Mathematics.

See Also

gaussLegendre, gaussHermite

Examples

cc <- gaussLaguerre(7)
# integrate exp(-x) from 0 to Inf
sum(cc$w)                     # 1
# integrate x^2 * exp(-x)     # integral x^n * exp(-x) is n!
sum(cc$w * cc$x^2)            # 2
# integrate sin(x) * exp(-x)
cc <- gaussLaguerre(17, 0)    # we need more nodes
sum(cc$w * sin(cc$x))         #=> 0.499999999994907 , should be 0.5

Gauss-Legendre Quadrature Formula

Description

Nodes and weights for the n-point Gauss-Legendre quadrature formula.

Usage

gaussLegendre(n, a, b)

Arguments

Details

x and w are obtained from a tridiagonal eigenvalue problem.

Value

List with components x, the nodes or points in[a,b], and w, the weights applied at these nodes.

Note

Gauss quadrature is not suitable for functions with singularities.

References

Gautschi, W. (2004). Orthogonal Polynomials: Computation and Approximation. Oxford University Press.

Trefethen, L. N. (2000). Spectral Methods in Matlab. SIAM, Society for Industrial and Applied Mathematics.

See Also

gaussHermite, gaussLaguerre

Examples

##  Quadrature with Gauss-Legendre nodes and weights
f <- function(x) sin(x+cos(10*exp(x))/3)
#\dontrun{ezplot(f, -1, 1, fill = TRUE)}
cc <- gaussLegendre(51, -1, 1)
Q <- sum(cc$w * f(cc$x))  #=> 0.0325036515865218 , true error: < 1e-15

# If f is not vectorized, do an explicit summation:
Q <- 0; x <- cc$x; w <- cc$w
for (i in 1:51) Q <- Q + w[i] * f(x[i])

# If f is infinite at b = 1, set  b <- b - eps  (with, e.g., eps = 1e-15)

# Use Gauss-Kronrod approach for error estimation
cc <- gaussLegendre(103, -1, 1)
abs(Q - sum(cc$w * f(cc$x)))     # rel.error < 1e-10

# Use Gauss-Hermite for vector-valued functions
f <- function(x) c(sin(pi*x), exp(x), log(1+x))
cc <- gaussLegendre(32, 0, 1)
drop(cc$w %*% matrix(f(cc$x), ncol = 3))  # c(2/pi, exp(1) - 1, 2*log(2) - 1)
# absolute error < 1e-15

Gauss-Newton Function Minimization

Description

Gauss-Newton method of minimizing a term f_1(x)^2 + \ldots + f_m(x)^2 or F' F where F = (f_1, \ldots, f_m) is a multivariate function of n variables, not necessarily n = m.

Usage

gaussNewton(x0, Ffun, Jfun = NULL,
                        maxiter =100, tol = .Machine$double.eps^(1/2), ...)

Arguments

Details

Solves the system of equations applying the Gauss-Newton's method. It is especially designed for minimizing a sum-of-squares of functions and can be used to find a common zero of several function.

This algorithm is described in detail in the textbook by Antoniou and Lu, incl. different ways to modify and remedy the Hessian if not being positive definite. Here, the approach by Goldfeld, Quandt and Trotter is used, and the hessian modified by the Matthews and Davies algorithm if still not invertible.

To accelerate the iteration, an inexact linesearch is applied.

Value

List with components:
xs the minimum or root found so far,
fs the square root of sum of squares of the values of f,
iter the number of iterations needed, and
relerr the absoulte distance between the last two solutions.

Note

If n=m then directly applying the newtonsys function might be a better alternative.

References

Antoniou, A., and W.-S. Lu (2007). Practical Optimization: Algorithms and Engineering Applications. Springer Business+Science, New York.

See Also

newtonsys, softline

Examples

f1 <- function(x) c(x[1]^2 + x[2]^2 - 1, x[1] + x[2] - 1)
gaussNewton(c(4, 4), f1)

f2 <- function(x) c( x[1] + 10*x[2], sqrt(5)*(x[] - x[4]),
                    (x[2] - 2*x[3])^2, 10*(x[1] - x[4])^2)
gaussNewton(c(-2, -1, 1, 2), f2)

f3 <- function(x)
        c(2*x[1] - x[2] - exp(-x[1]), -x[1] + 2*x[2] - exp(-x[2]))
gaussNewton(c(0, 0), f3)
# $xs   0.5671433 0.5671433

f4 <- function(x)  # Dennis Schnabel
        c(x[1]^2 + x[2]^2 - 2, exp(x[1] - 1) + x[2]^3 - 2)
gaussNewton(c(2.0, 0.5), f4)
# $xs    1 1

##  Examples (from Matlab)
F1 <- function(x) c(2*x[1]-x[2]-exp(-x[1]), -x[1]+2*x[2]-exp(-x[2]))
gaussNewton(c(-5, -5), F1)

# Find a matrix X such that X %*% X %*% X = [1 2; 3 4]
F2 <- function(x) {
    X <- matrix(x, 2, 2)
    D <- X %*% X %*% X - matrix(c(1,3,2,4), 2, 2)
    return(c(D))
}
sol <- gaussNewton(ones(2,2), F2)
(X  <- matrix(sol$xs, 2, 2))
#            [,1]      [,2]
# [1,] -0.1291489 0.8602157
# [2,]  1.2903236 1.1611747
X %*% X %*% X

GCD and LCM Integer Functions

Description

Greatest common divisor and least common multiple

Usage

gcd(a, b, extended = FALSE)
Lcm(a, b)

Arguments

Details

Computation based on the extended Euclidean algorithm.

If both a and b are vectors of the same length, the greatest common divisor/lowest common multiple will be computed elementwise. If one is a vektor, the other a scalar, the scalar will be replicated to the same length.

Value

A numeric (integer) value or vector of integers. Or a list of three vectors named c, d, g, g containing the greatest common divisors, such that

g = c * a + d * b.

Note

The following relation is always true:

n * m = gcd(n, m) * lcm(n, m)

See Also

numbers::extGCD

Examples

gcd(12, 1:24)
gcd(46368, 75025)  # Fibonacci numbers are relatively prime to each other
Lcm(12, 1:24)
Lcm(46368, 75025)  # = 46368 * 75025

Geometric Median

Description

Compute the “geometric median” of points in n-dimensional space, that is the point with the least sum of (Euclidean) distances to all these points.

Usage

geo_median(P, tol = 1e-07, maxiter = 200)

Arguments

Details

The task is solved applying an iterative process, known as Weiszfeld's algorithm. The solution is unique whenever the points are not collinear.

If the dimension is 1 (one column), the median will be returned.

Value

Returns a list with components p the coordinates of the solution point, d the sum of distances to all the sample points, reltol the relative tolerance of the iterative process, and niter the number of iterations.

Note

This is also known as the “1-median problem” and can be generalized to the “k-median problem” for k cluster centers; see kcca in the ‘flexclust’ package.

References

See Wikipedia's entry on “Geometric median”.

See Also

L1linreg

Examples

# Generate 100 points on the unit sphere in the 10-dim. space
set.seed(1001)
P <- rands(n=100, N=9)
( sol <- geo_median(P) )
# $p
#  [1] -0.009481361 -0.007643410 -0.001252910  0.006437703 -0.019982885 -0.045337987
#  [7]  0.036249563  0.003232175  0.035040592  0.046713023
# $d
# [1] 99.6638
# $reltol
# [1] 3.069063e-08
# $niter
# [1] 10

Geometric and Harmonic Mean (Matlab Style)

Description

Geometric and harmonic mean along a dimension of a vector, matrix, or array.
trimmean is almost the same as mean in R.

Usage

geomean(x, dim = 1)
harmmean(x, dim = 1)

trimmean(x, percent = 0)

Arguments

Details

trimmean does not call mean with the trim option, but rather calculates k<-round(n*percent/100/2) and leaves out k values at the beginning and end of the sorted x vector (or row or column of a matrix).

Value

Returns a scalar or vector (or array) of geometric or harmonic means: For dim=1 the mean of columns, dim=2 the mean of rows, etc.

Note

To have an exact analogue of mean(x) in Matlab, apply trimmean(x).

See Also

mean

Examples

A <- matrix(1:12, 3, 4)
geomean(A, dim = 1)
## [1]  1.817121  4.932424  7.958114 10.969613
harmmean(A, dim = 2)
## [1] 2.679426 4.367246 5.760000

x <- c(-0.98, -0.90, -0.68, -0.61, -0.61, -0.38, -0.37, -0.32, -0.20, -0.16,
        0.00,  0.05,  0.12,  0.30,  0.44,  0.77,  1.37,  1.64,  1.72,  2.80)
trimmean(x); trimmean(x, 20)    # 0.2  0.085
mean(x); mean(x, 0.10)          # 0.2  0.085

Givens Rotation

Description

Givens Rotations and QR decomposition

Usage

givens(A)

Arguments

Details

givens(A) returns a QR decomposition (or factorization) of the square matrix A by applying unitary 2-by-2 matrices U such that U * [xk;xl] = [x,0] where x=sqrt(xk^2+xl^2)

Value

List with two matrices Q and R, Q orthonormal and R upper triangular, such that A=Q%*%R.

References

Golub, G. H., and Ch. F. van Loan (1996). Matrix Computations. Third edition, John Hopkins University Press, Baltimore.

See Also

householder

Examples

##  QR decomposition
A <- matrix(c(0,-4,2, 6,-3,-2, 8,1,-1), 3, 3, byrow=TRUE)
gv <- givens(A)
(Q <- gv$Q); (R <- gv$R)
zapsmall(Q %*% R)

givens(magic(5))

Generalized Minimal Residual Method

Description

gmres(A,b) attempts to solve the system of linear equations A*x=b for x.

Usage

    gmres(A, b, x0 = rep(0, length(b)), 
          errtol = 1e-6, kmax = length(b)+1, reorth = 1)

Arguments

Details

Iterative method for the numerical solution of a system of linear equations. The method approximates the solution by the vector in a Krylov subspace with minimal residual. The Arnoldi iteration is used to find this vector.

Reorthogonalization method:
1 – Brown/Hindmarsh condition (default)
2 – Never reorthogonalize (not recommended)
3 – Always reorthogonalize (not cheap!)

Value

Returns a list with components x the solution, error the vector of residual norms, and niter the number of iterations.

Author(s)

Based on Matlab code from C. T. Kelley's book, see references.

References

C. T. Kelley (1995). Iterative Methods for Linear and Nonlinear Equations. SIAM, Society for Industrial and Applied Mathematics, Philadelphia, USA.

See Also

solve

Examples

A <- matrix(c(0.46, 0.60, 0.74, 0.61, 0.85,
              0.56, 0.31, 0.80, 0.94, 0.76,
              0.41, 0.19, 0.15, 0.33, 0.06,
              0.03, 0.92, 0.15, 0.56, 0.08,
              0.09, 0.06, 0.69, 0.42, 0.96), 5, 5)
x <- c(0.1, 0.3, 0.5, 0.7, 0.9)
b <- A %*% x
gmres(A, b)
# $x
#      [,1]
# [1,]  0.1
# [2,]  0.3
# [3,]  0.5
# [4,]  0.7
# [5,]  0.9
# 
# $error
# [1] 2.37446e+00 1.49173e-01 1.22147e-01 1.39901e-02 1.37817e-02 2.81713e-31
# 
# $niter
# [1] 5

Golden Ratio Search

Description

Golden Ratio search for a univariate function minimum in a bounded interval.

Usage

golden_ratio(f, a, b, ..., maxiter = 100, tol = .Machine$double.eps^0.5)

Arguments

Details

‘Golden ratio’ search for a univariate function minimum in a bounded interval.

Value

Return a list with components xmin, fmin, the function value at the minimum, niter, the number of iterations done, and the estimated precision estim.prec

See Also

uniroot

Examples

f <- function(x) x * cos(0.1*exp(x)) * sin(0.1*pi*exp(x))
golden_ratio(f, 0, 4, tol=10^-10)  # $xmin    = 3.24848329206212
optimize(f, c(0,4), tol=10^-10)    # $minimum = 3.24848328971188

Numerical Gradient

Description

Numerical function gradient.

Usage

grad(f, x0, heps = .Machine$double.eps^(1/3), ...)

Arguments

Details

Computes the gradient

(\frac{\partial f}{\partial x_1}, \ldots, \frac{\partial f}{\partial x_n})

numerically using the “central difference formula”.

Value

Vector of the same length as x0.

References

Mathews, J. H., and K. D. Fink (1999). Numerical Methods Using Matlab. Third Edition, Prentice Hall.

See Also

fderiv

Examples

f <- function(u) {
    x <- u[1]; y <- u[2]; z <- u[3]
    return(x^3 + y^2 + z^2 +12*x*y + 2*z)
 }
x0 <- c(1,1,1)
grad(f, x0)     # 15 14  4        # direction of steepest descent

sum(grad(f, x0) * c(1, -1, 0))    # 1 , directional derivative

f <- function(x) x[1]^2 + x[2]^2
grad(f, c(0,0))                   # 0 0 , i.e. a local optimum

Discrete Gradient (Matlab Style)

Description

Discrete numerical gradient.

Usage

gradient(F, h1 = 1, h2 = 1)

Arguments

Details

Returns the numerical gradient of a vector or matrix as a vector or matrix of discrete slopes in x- (i.e., the differences in horizontal direction) and slopes in y-direction (the differences in vertical direction).

A single spacing value, h, specifies the spacing between points in every direction, where the points are assumed equally spaced.

Value

If F is a vector, one gradient vector will be returned.

If F is a matrix, a list with two components will be returned:

where each matrix is of the same size as F.

Note

TODO: If h2 is missing, it will not automatically be adapted.

See Also

fderiv

Examples

x <- seq(0, 1, by=0.2)
y <- c(1, 2, 3)
(M <- meshgrid(x, y))
gradient(M$X^2 + M$Y^2)
gradient(M$X^2 + M$Y^2, x, y)

## Not run: 
# One-dimensional example
x <- seq(0, 2*pi, length.out = 100)
y <- sin(x)
f <- gradient(y, x)
max(f - cos(x))      #=> 0.00067086
plot(x, y, type = "l", col = "blue")
lines(x, cos(x), col = "gray", lwd = 3)
lines(x, f, col = "red")
grid()

# Two-dimensional example
v <- seq(-2, 2, by=0.2)
X <- meshgrid(v, v)$X
Y <- meshgrid(v, v)$Y

Z <- X * exp(-X^2 - Y^2)
image(v, v, t(Z))
contour(v, v, t(Z), col="black", add = TRUE)
grid(col="white")

grX <- gradient(Z, v, v)$X
grY <- gradient(Z, v, v)$Y

quiver(X, Y, grX, grY, scale = 0.2, col="blue")

## End(Not run)

Gram-Schmidt

Description

Modified Gram-Schmidt Process

Usage

  gramSchmidt(A, tol = .Machine$double.eps^0.5)

Arguments

Details

The modified Gram-Schmidt process uses the classical orthogonalization process to generate step by step an orthonoral basis of a vector space. The modified Gram-Schmidt iteration uses orthogonal projectors in order ro make the process numerically more stable.

Value

List with two matrices Q and R, Q orthonormal and R upper triangular, such that A=Q%*%R.

References

Trefethen, L. N., and D. Bau III. (1997). Numerical Linear Algebra. SIAM, Society for Industrial and Applied Mathematics, Philadelphia.

See Also

householder, givens

Examples

##  QR decomposition
A <- matrix(c(0,-4,2, 6,-3,-2, 8,1,-1), 3, 3, byrow=TRUE)
gs <- gramSchmidt(A)
(Q <- gs$Q); (R <- gs$R)
Q %*% R  # = A

Hadamard Matrix

Description

Generate Hadamard matrix of a certain size.

Usage

hadamard(n)

Arguments

Details

An n-by-n Hadamard matrix with n>2 exists only if rem(n,4)=0. This function handles only the cases where n, n/12, or n/20 is a power of 2.

Value

Matrix of size n-by-n of orthogonal columns consisting of 1 and -1 only.

Note

Hadamard matrices have applications in combinatorics, signal processing, and numerical analysis.

See Also

hankel, Toeplitz

Examples

hadamard(4)
H <- hadamard(8)
t(H) 

Halley's Root Finding Mathod

Description

Finding roots of univariate functions using the Halley method.

Usage

halley(fun, x0, maxiter = 500, tol = 1e-08, ...)

Arguments

Details

Well known root finding algorithms for real, univariate, continuous functions; the second derivative must be smooth, i.e. continuous. The first and second derivative are computed numerically.

Value

Return a list with components root, f.root, the function value at the found root, iter, the number of iterations done, and the estimated precision estim.prec

References

https://mathworld.wolfram.com/HalleysMethod.html

See Also

newtonRaphson

Examples

halley(sin, 3.0)        # 3.14159265358979 in 3 iterations
halley(function(x) x*exp(x) - 1, 1.0)
                        # 0.567143290409784 Gauss' omega constant

# Legendre polynomial of degree 5
lp5 <- c(63, 0, -70, 0, 15, 0)/8
f <- function(x) polyval(lp5, x)
halley(f, 1.0)          # 0.906179845938664

Hampel Filter

Description

Median absolute deviation (MAD) outlier in Time Series

Usage

hampel(x, k, t0 = 3)

Arguments

Details

The ‘median absolute deviation’ computation is done in the [-k...k] vicinity of each point at least k steps away from the end points of the interval. At the lower and upper end the time series values are preserved.

A high threshold makes the filter more forgiving, a low one will declare more points to be outliers. t0<-3 (the default) corresponds to Ron Pearson's 3 sigma edit rule, t0<-0 to John Tukey's median filter.

Value

Returning a list L with L$y the corrected time series and L$ind the indices of outliers in the ‘median absolut deviation’ sense.

Note

Don't take the expression outlier too serious. It's just a hint to values in the time series that appear to be unusual in the vicinity of their neighbors under a normal distribution assumption.

References

Pearson, R. K. (1999). “Data cleaning for dynamic modeling and control”. European Control Conference, ETH Zurich, Switzerland.

See Also

findpeaks

Examples

set.seed(8421)
x <- numeric(1024)
z <- rnorm(1024)
x[1] <- z[1]
for (i in 2:1024) {
	x[i] <- 0.4*x[i-1] + 0.8*x[i-1]*z[i-1] + z[i]
}
omad <- hampel(x, k=20)

## Not run: 
plot(1:1024, x, type="l")
points(omad$ind, x[omad$ind], pch=21, col="darkred")
grid()
## End(Not run)

Hankel Matrix

Description

Generate Hankel matrix from column and row vector

Usage

  hankel(a, b)

Arguments

Details

hankel(a) returns the square Hankel matrix whose first column is a and whose elements are zero below the secondary diagonal. (I.e.: b may be missing.)

hankel(a, b) returns a Hankel matrix whose first column is a and whose last row is b. If the first element of b differs from the last element of a it is overwritten by this one.

Value

matrix of size (length(a), length(b))

See Also

Toeplitz, hadamard

Examples

hankel(1:5, 5:1)

Hausdorff Distance

Description

Hausdorff distance (aka Hausdorff dimension)

Usage

hausdorff_dist(P, Q)

Arguments

Details

Calculates the Hausdorff Distance between two sets of points, P and Q. Sets P and Q must be matrices with the same number of columns (dimensions).

The ‘directional’ Hausdorff distance (dhd) is defined as:

dhd(P,Q) = max p in P [ min q in Q [ ||p-q|| ] ]

Intuitively dhd finds the point p from the set P that is farthest from any point in Q and measures the distance from p to its nearest neighbor in Q. The Hausdorff Distance is defined as max(dhd(P,Q),dhd(Q,P)).

Value

A single scalar, the Hausdorff distance (dimension).

References

Barnsley, M. (1993). Fractals Everywhere. Morgan Kaufmann, San Francisco.

See Also

distmat

Examples

P <- matrix(c(1,1,2,2, 5,4,5,4), 4, 2)
Q <- matrix(c(4,4,5,5, 2,1,2,1), 4, 2)
hausdorff_dist(P, Q)    # 4.242641 = sqrt(sum((c(4,2)-c(1,5))^2))

Haversine Formula

Description

Haversine formula to calculate the arc distance between two points on earth (i.e., along a great circle).

Usage

haversine(loc1, loc2, R = 6371.0)

Arguments

Details

The Haversine formula is more robust for the calculating the distance as with the spherical cosine formula. The user may want to assume a slightly different earth radius, so this can be provided as input.

The location can be input in two different formats, as latitude and longitude in a character string, e.g. for Frankfurt airport as '50 02 00N, 08 34 14E', or as a numerical two-vector in degrees (not radians).

Here for latitude 'N' and 'S' stand for North and South, and for longitude 'E' or 'W' stand for East and West. For the degrees format, South and West must be negative.

These two formats can be mixed.

Value

Returns the distance in km.

Author(s)

Hans W. Borchers

References

Entry 'Great_circle_distance' in Wikipedia.

See Also

Implementations of the Haversine formula can also be found in other R packages, e.g. 'geoPlot' or 'geosphere'.

Examples

FRA = '50 02 00N, 08 34 14E'  # Frankfurt Airport
ORD = '41 58 43N, 87 54 17W'  # Chicago O'Hare Interntl. Airport
fra <- c(50+2/60, 8+34/60+14/3600)
ord <- c(41+58/60+43/3600, -(87+54/60+17/3600))

dis <- haversine(FRA, ORD)    # 6971.059 km
fprintf('Flight distance Frankfurt-Chicago is %8.3f km.\n', dis)

dis <- haversine(fra, ord)
fprintf('Flight distance Frankfurt-Chicago is %8.3f km.\n', dis)

Hessenberg Matrix

Description

Generates the Hessenberg matrix for A.

Usage

hessenberg(A)

Arguments

Details

An (upper) Hessenberg matrix has zero entries below the first subdiagonal.

The function generates a Hessenberg matrix H and a unitary matrix P (a similarity transformation) such that A = P * H * t(P).

The Hessenberg matrix has the same eigenvalues. If A is symmetric, its Hessenberg form will be a tridiagonal matrix.

Value

Returns a list with two elements,

References

Press, Teukolsky, Vetterling, and Flannery (2007). Numerical Recipes: The Art of Scientific Computing. 3rd Edition, Cambridge University Press. (Section 11.6.2)

See Also

householder

Examples

A <- matrix(c(-149,   -50,  -154,
               537,   180,   546,
               -27,    -9,   -25), nrow = 3, byrow = TRUE)
hb  <- hessenberg(A)
hb
## $H
##           [,1]         [,2]        [,3]
## [1,] -149.0000  42.20367124 -156.316506
## [2,] -537.6783 152.55114875 -554.927153
## [3,]    0.0000   0.07284727    2.448851
## 
## $P
##      [,1]       [,2]      [,3]
## [1,]    1  0.0000000 0.0000000
## [2,]    0 -0.9987384 0.0502159
## [3,]    0  0.0502159 0.9987384

hb$P %*% hb$H %*% t(hb$P)
##      [,1] [,2] [,3]
## [1,] -149  -50 -154
## [2,]  537  180  546
## [3,]  -27   -9  -25

Hessian Matrix

Description

Numerically compute the Hessian matrix.

Usage

hessian(f, x0, h = .Machine$double.eps^(1/4), ...)

Arguments

Details

Computes the hessian matrix based on the three-point central difference formula, expanded to two variables.

Assumes that the function has continuous partial derivatives.

Value

An n-by-n matrix with \frac{\partial^2 f}{\partial x_i \partial x_j} as (i, j) entry.

References

Fausett, L. V. (2007). Applied Numerical Analysis Using Matlab. Second edition, Prentice Hall.

See Also

hessdiag, hessvec, laplacian

Examples

f <- function(x) cos(x[1] + x[2])
x0 <- c(0, 0)
hessian(f, x0)

f <- function(u) {
    x <- u[1]; y <- u[2]; z <- u[3]
    return(x^3 + y^2 + z^2 +12*x*y + 2*z)
}
x0 <- c(1,1,1)
hessian(f, x0)

Hessian utilities

Description

Fast multiplication of Hessian and vector where computation of the full Hessian is not needed. Or determine the diagonal of the Hessian when non-diagonal entries are not needed or are nearly zero.

Usage

  hessvec(f, x, v, csd = FALSE, ...)

  hessdiag(f, x, ...)

Arguments

Details

hessvec computes the product of a Hessian of a function times a vector without deriving the full Hessian by approximating the gradient (see the reference). If the function allows for the complex-step method, the gradient can be calculated much more accurate (see grad_csd).

hessdiag computes only the diagonal of the Hessian by applying the central difference formula of second order to approximate the partial derivatives.

Value

hessvec returns the product H(f,x) * v as a vector.

hessdiag returns the diagonal of the Hessian of f.

References

B.A. Pearlmutter, Fast Exact Multiplication by the Hessian, Neural Computation (1994), Vol. 6, Issue 1, pp. 147-160.

See Also

hessian

Examples

  ## Not run: 
    set.seed(1237); n <- 100
    a <- runif(n); b <- rnorm(n)
    fn <- function(x, a, b) sum(exp(-a*x)*sin(b*pi*x))
    x0 <- rep(1, n)
    v0 <- rexp(n, rate=0.1)
    
    # compute with full hessian
    h0 <- hessian(fn, x0, a = a, b = b)             # n=100 runtimes
    v1 <- c(h0 %*% v0)                              # 0.167   sec
    
    v2 <- hessvec(fn, x0, v0, a = a, b = b)         # 0.00209 sec
    v3 <- hessvec(fn, x0, v0, csd=TRUE,a=a, b=b)    # 0.00145 sec
    v4 <- hessdiag(fn, x0, a = a, b = b) * v0       # 0.00204 sec
    
    # compare with exact analytical Hessian
    hex <- diag((a^2-b^2*pi^2)*exp(-a*x0)*sin(b*pi*x0) - 
                 2*a*b*pi*exp(-a*x0)*cos(b*pi*x0))
    vex <- c(hex %*% v0)

    max(abs(vex - v1))          # 2.48e-05
    max(abs(vex - v2))          # 7.15e-05
    max(abs(vex - v3))          # 0.09e-05
    max(abs(vex - v4))          # 2.46e-05 
## End(Not run)

Hilbert Matrix

Description

Generate Hilbert matrix of dimension n

Usage

  hilb(n)

Arguments

Details

Generate the Hilbert matrix H of dimension n with elements H[i, j] = 1/(i+j-1).

(Note: This matrix is ill-conditioned, see e.g. det(hilb(6)).)

Value

matrix of dimension n

See Also

vander

Examples

  hilb(5)

Histogram Count (Matlab style)

Description

Histogram-like counting.

Usage

  histc(x, edges)

Arguments

Details

n = histc(x,edges) counts the number of values in vector x that fall between the elements in the edges vector (which must contain monotonically nondecreasing values). n is a length(edges) vector containing these counts.

If x is a matrix then cnt and bin are matrices too, and

for (j in (1:n)) cnt[k,j] <- sum(bin[, j] == k)

Value

returns a list with components cnt and bin. n(k) counts the number of values in x that lie between edges(k) <= x(i) < edges(k+1). The last counts any values of x that match edges(n). Values outside the values in edges are not counted. Use -Inf and Inf in edges to include all values.

bin[i] returns k if edges(k) <= x(i) < edges(k+1), and 0 if x[i] lies outside the grid.

See Also

hist, histss, findInterval

Examples

x <- seq(0.0, 1.0, by = 0.05)
e <- seq(0.1, 0.9, by = 0.10)
histc(x, e)
# $cnt
# [1] 2 2 2 2 2 2 2 2 1
# $bin
# [1] 0 0 1 1 2 2 3 3 4 4 5 5 6 6 7 7 8 8 9 0 0

## Not run: 
# Compare
findInterval(x, e)
# [1] 0 0 1 1 2 2 3 3 4 4 5 5 6 6 7 7 8 8 9 9 9
findInterval(x, e, all.inside = TRUE)
# [1] 1 1 1 1 2 2 3 3 4 4 5 5 6 6 7 7 8 8 8 8 8
# cnt[i] <- sum(findInterval(x, e) == i)
## End(Not run)

x <- matrix( c(0.5029, 0.2375, 0.2243, 0.8495,
               0.0532, 0.1644, 0.4215, 0.4135,
               0.7854, 0.0879, 0.1221, 0.6170), 3, 4, byrow = TRUE)
e <- seq(0.0, 1.0, by = 0.2)
histc(x, e)
# $cnt
#      [,1] [,2] [,3] [,4]
# [1,]    1    2    1    0
# [2,]    0    1    1    0
# [3,]    1    0    1    1
# [4,]    1    0    0    1
# [5,]    0    0    0    1
# [6,]    0    0    0    0
# 
# $bin
#      [,1] [,2] [,3] [,4]
# [1,]    3    2    2    5
# [2,]    1    1    3    3
# [3,]    4    1    1    4

Histogram Bin-width Optimization

Description

Method for selecting the bin size of time histograms.

Usage

histss(x, n = 100, plotting = FALSE)

Arguments

Details

Bin sizes of histograms are optimized in a way to best displays the underlying spike rate, for example in neurophysiological studies.

Value

Returns the same list as the hist function; the list is invisible if the histogram is plotted.

References

Shimazaki H. and S. Shinomoto. A method for selecting the bin size of a time histogram. Neural Computation (2007) Vol. 19(6), 1503-1527

See Also

hist, histc

Examples

x <- sin(seq(0, pi/2, length.out = 200))
H <- histss(x, n = 50, plotting = FALSE)
## Not run: 
plot(H, col = "gainsboro")  # Compare with hist(x), or
hist(x, breaks = H$breaks)  # the same 
## End(Not run)

Hooke-Jeeves Function Minimization Method

Description

An implementation of the Hooke-Jeeves algorithm for derivative-free optimization.

Usage

hooke_jeeves(x0, fn, ..., lb = NULL, ub = NULL, tol = 1e-08,
             maxfeval = 10000, target = Inf, info = FALSE)

Arguments

Details

This method computes a new point using the values of f at suitable points along the orthogonal coordinate directions around the last point.

Value

List with following components:

Note

Hooke-Jeeves is notorious for its number of function calls. Memoization is often suggested as a remedy.

For a similar implementation of Hooke-Jeeves see the ‘dfoptim’ package.

References

C.T. Kelley (1999), Iterative Methods for Optimization, SIAM.

Quarteroni, Sacco, and Saleri (2007), Numerical Mathematics, Springer-Verlag.

See Also

nelder_mead

Examples

##  Rosenbrock function
rosenbrock <- function(x) {
    n <- length(x)
    x1 <- x[2:n]
    x2 <- x[1:(n-1)]
    sum(100*(x1-x2^2)^2 + (1-x2)^2)
}

hooke_jeeves(c(0,0,0,0), rosenbrock)
## $xmin
## [1] 1.000002 1.000003 1.000007 1.000013
## $fmin
## [1] 5.849188e-11
## $count
## [1] 1691
## $convergence
## [1] 0
## $info
## $info$solver
## [1] "Hooke-Jeeves"
## $info$iterations
## [1] 26

hooke_jeeves(rep(0,4), lb=rep(-1,4), ub=0.5, rosenbrock)
## $xmin
## [1] 0.50000000 0.26221320 0.07797602 0.00608027
## $fmin
## [1] 1.667875
## $count
## [1] 536
## $convergence
## [1] 0
## $info
## $info$solver
## [1] "Hooke-Jeeves"
## $info$iterations
## [1] 26

Horner's Rule

Description

Compute the value of a polynomial via Horner's Rule.

Usage

horner(p, x)
hornerdefl(p, x)

Arguments

Details

horner utilizes the Horner scheme to evaluate the polynomial and its first derivative at the same time.

The polynomial p = p_1*x^n + p_2*x^{n-1} + ... + p_n*x + p_{n+1} is hereby represented by the vector p_1, p_2, ..., p_n, p_{n+1}, i.e. from highest to lowest coefficient.

hornerdefl uses a similar approach to return the value of p at x and a polynomial q that satisfies

p(t) = q(t) * (t - x) + r, r constant

which implies r=0 if x is a root of p. This will allow for a repeated root finding of polynomials.

Value

horner returns a list with two elements, list(y=..., dy=...) where the first list elements returns the values of the polynomial, the second the values of its derivative at the point(s) x.

hornerdefl returns a list list(y=..., dy=...) where q represents a polynomial, see above.

Note

For fast evaluation, there is no error checking for p and x, which both must be numerical vectors (x can be a matrix in horner).

References

Quarteroni, A., and Saleri, F. (2006) Scientific Computing with Matlab and Octave. Second Edition, Springer-Verlag, Berlin Heidelberg.

See Also

polyval

Examples

x <- c(-2, -1, 0, 1, 2)
p <- c(1, 0, 1)  # polynomial x^2 + x, derivative 2*x
horner(p, x)$y   #=>  5  2  1  2  5
horner(p, x)$dy  #=> -4 -2  0  2  4

p <- Poly(c(1, 2, 3))  # roots 1, 2, 3
hornerdefl(p, 3)          # q = x^2- 3 x + 2  with roots 1, 2

Householder Reflections

Description

Householder reflections and QR decomposition

Usage

  householder(A)

Arguments

Details

The Householder method applies a succession of elementary unitary matrices to the left of matrix A. These matrices are the so-called Householder reflections.

Value

List with two matrices Q and R, Q orthonormal and R upper triangular, such that A=Q%*%R.

References

Trefethen, L. N., and D. Bau III. (1997). Numerical Linear Algebra. SIAM, Society for Industrial and Applied Mathematics, Philadelphia.

See Also

givens

Examples

##  QR decomposition
A <- matrix(c(0,-4,2, 6,-3,-2, 8,1,-1), 3, 3, byrow=TRUE)
S <- householder(A)
(Q <- S$Q); (R <- S$R)
Q %*% R  # = A

##  Solve an overdetermined linear system of equations
A <- matrix(c(1:8,7,4,2,3,4,2,2), ncol=3, byrow=TRUE)
S <- householder(A); Q <- S$Q; R <- S$R
m <- nrow(A); n <- ncol(A)
b <- rep(6, 5)

x <- numeric(n)
b <- t(Q) %*% b
x[n] <- b[n] / R[n, n]
for (k in (n-1):1)
    x[k] <- (b[k] - R[k, (k+1):n] %*% x[(k+1):n]) / R[k, k]
qr.solve(A, rep(6, 5)); x

Matlab Test Functions

Description

Matlab test functions.

Usage

humps(x)
sinc(x)
psinc(x, n)

Arguments

Details

humps is a test function for finding zeros, for optimization and integration. Its root is at x = 1.2995, a (local) minimum at x = 0.6370, and the integral from 0.5 to 1.0 is 8.0715.

sinc is defined as sinc(t) = \frac{\sin(\pi t)}{\pi t}. It is the continuous inverse Fourier transform of the rectangular pulse of width 2\pi and height 1.

psinc is the 'periodic sinc function' and is defined as psinc(x,n) = \frac{\sin(x n/2)}{n \sin(x/2)}.

Value

Numeric scalar or vector.

Examples

## Not run: 
plot(humps(), type="l"); grid()

x <- seq(0, 10, length=101)
plot(x, sinc(x), type="l"); grid()

## End(Not run)

Hurst Exponent

Description

Calculates the Hurst exponent using R/S analysis.

Usage

  hurstexp(x, d = 50, display = TRUE)

Arguments

Details

hurstexp(x) calculates the Hurst exponent of a time series x using R/S analysis, after Hurst, with slightly different approaches, or corrects it with small sample bias, see for example Weron.

These approaches are a corrected R/S method, an empirical and corrected empirical method, and a try at a theoretical Hurst exponent. It should be mentioned that the results are sometimes very different, so providing error estimates will be highly questionable.

Optimal sample sizes are automatically computed with a length that possesses the most divisors among series shorter than x by no more than 1 percent.

Value

hurstexp(x) returns a list with the following components:

• Hs - simplified R over S approach

• Hrs - corrected R over S Hurst exponent

• He - empirical Hurst exponent

• Hal - corrected empirical Hurst exponent

• Ht - theoretical Hurst exponent

Note

Derived from Matlab code of R. Weron, published on Matlab Central.

References

H.E. Hurst (1951) Long-term storage capacity of reservoirs, Transactions of the American Society of Civil Engineers 116, 770-808.

R. Weron (2002) Estimating long range dependence: finite sample properties and confidence intervals, Physica A 312, 285-299.

See Also

fractal::hurstSpec, RoverS, hurstBlock and fArma::LrdModelling

Examples

##  Computing the Hurst exponent
data(brown72)
x72 <- brown72                          #  H = 0.72
xgn <- rnorm(1024)                      #  H = 0.50
xlm <- numeric(1024); xlm[1] <- 0.1     #  H = 0.43
for (i in 2:1024) xlm[i] <- 4 * xlm[i-1] * (1 - xlm[i-1])

hurstexp(brown72, d = 128)           # 0.72
# Simple R/S Hurst estimation:         0.6590931 
# Corrected R over S Hurst exponent:   0.7384611 
# Empirical Hurst exponent:            0.7068613 
# Corrected empirical Hurst exponent:  0.6838251 
# Theoretical Hurst exponent:          0.5294909

hurstexp(xgn)                        # 0.50
# Simple R/S Hurst estimation:         0.5518143 
# Corrected R over S Hurst exponent:   0.5982146 
# Empirical Hurst exponent:            0.6104621 
# Corrected empirical Hurst exponent:  0.5690305 
# Theoretical Hurst exponent:          0.5368124 

hurstexp(xlm)                        # 0.43
# Simple R/S Hurst estimation:         0.4825898 
# Corrected R over S Hurst exponent:   0.5067766 
# Empirical Hurst exponent:            0.4869625 
# Corrected empirical Hurst exponent:  0.4485892 
# Theoretical Hurst exponent:          0.5368124 


##  Compare with other implementations
## Not run: 
library(fractal)

x <- x72
hurstSpec(x)                    # 0.776   # 0.720
RoverS(x)                       # 0.717
hurstBlock(x, method="aggAbs")  # 0.648
hurstBlock(x, method="aggVar")  # 0.613
hurstBlock(x, method="diffvar") # 0.714
hurstBlock(x, method="higuchi") # 1.001

x <- xgn
hurstSpec(x)                    # 0.538   # 0.500
RoverS(x)                       # 0.663
hurstBlock(x, method="aggAbs")  # 0.463
hurstBlock(x, method="aggVar")  # 0.430
hurstBlock(x, method="diffvar") # 0.471
hurstBlock(x, method="higuchi") # 0.574

x <- xlm
hurstSpec(x)                    # 0.478   # 0.430
RoverS(x)                       # 0.622
hurstBlock(x, method="aggAbs")  # 0.316
hurstBlock(x, method="aggVar")  # 0.279
hurstBlock(x, method="diffvar") # 0.547
hurstBlock(x, method="higuchi") # 0.998

## End(Not run)

Hypotenuse Function

Description

Square root of sum of squares

Usage

hypot(x, y)

Arguments

Details

Element-by-element computation of the square root of the sum of squares of vectors resp. matrices x and y.

Value

Returns a vector or matrix of the same size.

Note

Returns c() if x or y is empty and the other one has length 1. If one input is scalar, the other a vector, the scalar will be extended to a vector of appropriate length. In all other cases, x and y have to be of the same size.

Examples

hypot(3,4)
hypot(1, c(3, 4, 5))
hypot(c(0, 0), c(3, 4))

Inverse Fast Fourier Transformation

Description

Performs the inverse Fast Fourier Transform.

Usage

ifft(x)

ifftshift(x)
fftshift(x)

Arguments

Details

ifft returns the value of the normalized discrete, univariate, inverse Fast Fourier Transform of the values in x.

ifftshift and fftshift shift the zero-component to the center of the spectrum, that is swap the left and right half of x.

Value

Real or complex vector of the same length.

Note

Almost an alias for R's fft(x, inverse=TRUE), but dividing by length(x).

See Also

fft

Examples

x <- c(1, 2, 3, 4)
(y <- fft(x))
ifft(x)
ifft(y)

##  Compute the derivative: F(df/dt) = (1i*k) * F(f)
#   hyperbolic secans f <- sech
df <- function(x) -sech(x) * tanh(x)
d2f <- function(x) sech(x) - 2*sech(x)^3
L <- 20                                 # domain [-L/2, L/2]
N <- 128                                # number of Fourier nodes
x <- linspace(-L/2, L/2, N+1)           # domain discretization
x <- x[1:N]                             # because of periodicity
dx <- x[2] - x[1]                       # finite difference
u <- sech(x)                            # hyperbolic secans
u1d <- df(x); u2d <- d2f(x)             # first and second derivative
ut <- fft(u)                            # discrete Fourier transform
k <- (2*pi/L)*fftshift((-N/2):(N/2-1))  # shifted frequencies
u1 <- Re(ifft((1i*k) * ut))             # inverse transform
u2 <- Re(ifft(-k^2 * ut))               # first and second derivative
## Not run: 
plot(x, u1d, type = "l", col = "blue")
points(x, u1)
grid()
figure()
plot(x, u2d, type = "l", col = "darkred")
points(x, u2)
grid()
## End(Not run)

Polygon Region

Description

Points inside polygon region.

Usage

inpolygon(x, y, xp, yp, boundary = FALSE)

Arguments

Details

For a polygon defined by points (xp, yp), determine if the points (x, y) are inside or outside the polygon. The boundary can be included or excluded (default) for the interior.

Value

Logical vector, the same length as x.

Note

Special care taken for points on the boundary.

References

Hormann, K., and A. Agathos (2001). The Point in Polygon Problem for Arbitrary Polygons. Computational Geometry, Vol. 20, No. 3, pp. 131–144.

See Also

polygon

Examples

xp <- c(0.5, 0.75, 0.75, 0.5, 0.5)
yp <- c(0.5, 0.5, 0.75, 0.75, 0.5)
x <- c(0.6, 0.75, 0.6, 0.5)
y <- c(0.5, 0.6, 0.75, 0.6)
inpolygon(x, y, xp, yp, boundary = FALSE)  # FALSE
inpolygon(x, y, xp, yp, boundary = TRUE)   # TRUE

## Not run: 
pg <- matrix(c(0.15, 0.75, 0.25, 0.45, 0.70,
               0.80, 0.35, 0.55, 0.20, 0.90), 5, 2)
plot(c(0, 1), c(0, 1), type="n")
polygon(pg[,1], pg[,2])
P <- matrix(runif(20000), 10000, 2)
R <- inpolygon(P[, 1], P[, 2], pg[, 1], pg[,2])
clrs <- ifelse(R, "red", "blue")
points(P[, 1], P[, 2], pch = ".", col = clrs)
## End(Not run)

Adaptive Numerical Integration

Description

Combines several approaches to adaptive numerical integration of functions of one variable.

Usage

integral(fun, xmin, xmax,
         method = c("Kronrod", "Clenshaw","Simpson"),
         no_intervals = 8, random = FALSE,
         reltol = 1e-8, abstol = 0, ...)

Arguments

Details

integral combines the following methods for adaptive numerical integration (also available as separate functions):

• Kronrod (Gauss-Kronrod)

• Clenshaw (Clenshaw-Curtis; not yet made adaptive)

• Simpson (adaptive Simpson)

Recommended default method is Gauss-Kronrod. Also try Clenshaw-Curtis that may be faster at times.

Most methods require that function f is vectorized. This will be checked and the function vectorized if necessary.

By default, the integration domain is subdivided into no_intervals subdomains to avoid 0 results if the support of the integrand function is small compared to the whole domain. If random is true, nodes will be picked randomly, otherwise forming a regular division.

If the interval is infinite, quadinf will be called that accepts the same methods as well. [If the function is array-valued, quadv is called that applies an adaptive Simpson procedure, other methods are ignored – not true anymore.]

Value

Returns the integral, no error terms given.

Note

integral does not provide ‘new’ functionality, everything is already contained in the functions called here. Other interesting alternatives are Gauss-Richardson (quadgr) and Romberg (romberg) integration.

References

Davis, Ph. J., and Ph. Rabinowitz (1984). Methods of Numerical Integration. Dover Publications, New York.

See Also

quadgk, quadgr, quadcc, simpadpt, romberg, quadv, quadinf

Examples

##  Very smooth function
fun <- function(x) 1/(x^4+x^2+0.9)
val <- 1.582232963729353
for (m in c("Kron", "Clen", "Simp")) {
    Q <- integral(fun, -1, 1, reltol = 1e-12, method = m)
    cat(m, Q, abs(Q-val), "\n")}
# Kron 1.582233 3.197442e-13 
# Rich 1.582233 3.197442e-13  # use quadgr()
# Clen 1.582233 3.199663e-13 
# Simp 1.582233 3.241851e-13 
# Romb 1.582233 2.555733e-13  # use romberg()

##  Highly oscillating function
fun <- function(x) sin(100*pi*x)/(pi*x)
val <- 0.4989868086930458
for (m in c("Kron", "Clen", "Simp")) {
    Q <- integral(fun, 0, 1, reltol = 1e-12, method = m)
    cat(m, Q, abs(Q-val), "\n")}
# Kron 0.4989868 2.775558e-16 
# Rich 0.4989868 4.440892e-16  # use quadgr()
# Clen 0.4989868 2.231548e-14
# Simp 0.4989868 6.328271e-15 
# Romb 0.4989868 1.508793e-13  # use romberg()

## Evaluate improper integral
fun <- function(x) log(x)^2 * exp(-x^2)
val <- 1.9475221803007815976
Q <- integral(fun, 0, Inf, reltol = 1e-12)
# For infinite domains Gauss integration is applied!
cat(m, Q, abs(Q-val), "\n")
# Kron 1.94752218028062 2.01587635473288e-11 

## Example with small function support
fun <- function(x)
            ifelse (x <= 0 | x >= pi, 0, sin(x))
integral(fun, -100, 100, no_intervals = 1)      # 0
integral(fun, -100, 100, no_intervals = 10)     # 1.99999999723
integral(fun, -100, 100, random=FALSE)          # 2
integral(fun, -100, 100, random=TRUE)           # 2 (sometimes 0 !)
integral(fun, -1000, 10000, random=FALSE)       # 0
integral(fun, -1000, 10000, random=TRUE)        # 0 (sometimes 2 !)

Numerically Evaluate Double and Triple Integrals

Description

Numerically evaluate a double integral, resp. a triple integral by reducing it to a double integral.

Usage

integral2(fun, xmin, xmax, ymin, ymax, sector = FALSE,
            reltol = 1e-6, abstol = 0, maxlist = 5000,
            singular = FALSE, vectorized = TRUE, ...)

integral3(fun, xmin, xmax, ymin, ymax, zmin, zmax,
            reltol = 1e-6, ...)

Arguments

Details

integral2 implements the ‘TwoD’ algorithm, that is Gauss-Kronrod with (3, 7)-nodes on 2D rectangles.

The borders of the domain of integration must be finite. The limits of y, that is ymin and ymax, can be constants or scalar functions of x that describe the lower and upper boundaries. These functions must be vectorized.

integral2 attempts to satisfy ERRBND <= max(AbsTol,RelTol*|Q|). This is absolute error control when |Q| is sufficiently small and relative error control when |Q| is larger.

The function fun itself must be fully vectorized: It must accept arrays X and Y and return an array Z = f(X,Y) of corresponding values. If option vectorized is set to FALSE the procedure will enforce this vectorized behavior.

With sector=TRUE the region is a generalized sector that is described in polar coordinates (r,theta) by

0 <= a <= theta <= b – a and b must be constants
c <= r <= d – c and d can be constants or ...

... functions of theta that describe the lower and upper boundaries. Functions must be vectorized.
NOTE Polar coordinates are used only to describe the region – the integrand is f(x,y) for both kinds of regions.

integral2 can be applied to functions that are singular on a boundary. With value singular=TRUE, this option causes integral2 to use transformations to weaken singularities for better performance.

integral3 also accepts functions for the inner interval limits. ymin, ymax must be constants or functions of one variable (x), zmin, zmax constants or functions of two variables (x, y), all functions vectorized.

The triple integral will be first integrated over the second and third variable with integral2, and then integrated over a single variable with integral.

Value

Returns a list with Q the integral and error the error term.

Note

To avoid recursion, a possibly large matrix will be used and passed between subprograms. A more efficient implementation may be possible.

Author(s)

Copyright (c) 2008 Lawrence F. Shampine for Matlab code and description of the program; adapted and converted to R by Hans W Borchers.

References

Shampine, L. F. (2008). MATLAB Program for Quadrature in 2D. Proceedings of Applied Mathematics and Computation, 2008, pp. 266–274.

See Also

integral, cubature:adaptIntegrate

Examples

fun <- function(x, y) cos(x) * cos(y)
integral2(fun, 0, 1, 0, 1, reltol = 1e-10)
# $Q:     0.708073418273571  # 0.70807341827357119350 = sin(1)^2
# $error: 8.618277e-19       # 1.110223e-16

##  Compute the volume of a sphere
f <- function(x, y) sqrt(1 -x^2 - y^2)
xmin <- 0; xmax <- 1
ymin <- 0; ymax <- function(x) sqrt(1 - x^2)
I <- integral2(f, xmin, xmax, ymin, ymax)
I$Q                             # 0.5236076 - pi/6 => 8.800354e-06

##  Compute the volume over a sector
I <- integral2(f, 0,pi/2, 0,1, sector = TRUE)
I$Q                             # 0.5236308 - pi/6 => 3.203768e-05

##  Integrate 1/( sqrt(x + y)*(1 + x + y)^2 ) over the triangle
##   0 <= x <= 1, 0 <= y <= 1 - x.  The integrand is infinite at (0,0).
f <- function(x,y) 1/( sqrt(x + y) * (1 + x + y)^2 )
ymax <- function(x) 1 - x
I <- integral2(f, 0,1, 0,ymax)
I$Q + 1/2 - pi/4                # -3.247091e-08

##  Compute this integral as a sector
rmax <- function(theta) 1/(sin(theta) + cos(theta))
I <- integral2(f, 0,pi/2, 0,rmax, sector = TRUE, singular = TRUE)
I$Q + 1/2 - pi/4                # -4.998646e-11

##  Examples of computing triple integrals
f0 <- function(x, y, z) y*sin(x) + z*cos(x)
integral3(f0, 0, pi, 0,1, -1,1) # - 2.0 => 0.0

f1 <- function(x, y, z) exp(x+y+z)
integral3(f1, 0, 1, 1, 2, 0, 0.5)
## [1] 5.206447                         # 5.20644655

f2 <- function(x, y, z) x^2 + y^2 + z
a <- 2; b <- 4
ymin <- function(x) x - 1
ymax <- function(x) x + 6
zmin <- -2
zmax <- function(x, y) 4 + y^2
integral3(f2, a, b, ymin, ymax, zmin, zmax)
## [1] 47416.75556                      # 47416.7555556

f3 <- function(x, y, z) sqrt(x^2 + y^2)
a <- -2; b <- 2
ymin <- function(x) -sqrt(4-x^2)
ymax <- function(x)  sqrt(4-x^2)
zmin <- function(x, y)  sqrt(x^2 + y^2)
zmax <- 2
integral3(f3, a, b, ymin, ymax, zmin, zmax)
## [1] 8.37758                          # 8.377579076269617

One-dimensional Interpolation

Description

One-dimensional interpolation of points.

Usage

interp1(x, y, xi = x,
        method = c("linear", "constant", "nearest", "spline", "cubic"))

Arguments

Details

Interpolation to find yi, the values of the underlying function at the points in the vector xi.

Methods can be:

Value

Numeric vector representing values at points xi.

Note

Method ‘spline’ uses the spline approach by Moler et al., and is identical with the Matlab option of the same name, but slightly different from R's spline function.

The Matlab option “cubic” seems to have no direct correspondence in R. Therefore, we simply use pchip here.

See Also

approx, spline

Examples

x <- c(0.8, 0.3, 0.1, 0.6, 0.9, 0.5, 0.2, 0.0, 0.7, 1.0, 0.4)
y <- x^2
xi <- seq(0, 1, len = 81)
yl <- interp1(x, y, xi, method = "linear")
yn <- interp1(x, y, xi, method = "nearest")
ys <- interp1(x, y, xi, method = "spline")

## Not run: 
plot(x, y); grid()
lines(xi, yl, col="blue", lwd = 2)
lines(xi, yn, col="black", lty = 2)
lines(xi, ys, col="red")
  
## End(Not run)

## Difference between spline (Matlab) and spline (R).
x <- 1:6
y <- c(16, 18, 21, 17, 15, 12)
xs <- linspace(1, 6, 51)
ys <- interp1(x, y, xs, method = "spline")
sp <- spline(x, y, n = 51, method = "fmm")

## Not run: 
plot(x, y, main = "Matlab and R splines")
grid()
lines(xs, ys, col = "red")
lines(sp$x, sp$y, col = "blue")
legend(4, 20, c("Matlab spline", "R spline"), 
              col = c("red", "blue"), lty = 1)
  
## End(Not run)

Two-dimensional Data Interpolation

Description

Two-dimensional data interpolation similar to a table look-up.

Usage

    interp2(x, y, Z, xp, yp, method = c("linear", "nearest", "constant"))

Arguments

Details

Computes a vector containing elements corresponding to the elements of xp and yp, determining by interpolation within the two-dimensional function specified by vectors x and y, and matrix Z.

x and y must be monotonically increasing. They specify the points at which the data Z is given. Therefore, length(x) = nrow(Z) and length(y) = ncol(Z) must be satisfied.

xp and yp must be of the same length.

The functions appears vectorized as xp, yp can be vectors, but internally they are treated in a for loop.

Value

Vector the length of xp of interpolated values.

For methods “constant” and “nearest” the intervals are considered closed from left and below. Out of range values are returned as NAs.

Note

The corresponding Matlab function has also the methods “cubic” and “spline”. If in need of a nonlinear interpolation, take a look at barylag2d in this package and the example therein.

See Also

interp1, barylag2d

Examples

## Not run: 
    x <- linspace(-1, 1, 11)
    y <- linspace(-1, 1, 11)
    mgrid <- meshgrid(x, y)
    Z <- mgrid$X^2 + mgrid$Y^2
    xp <- yp <- linspace(-1, 1, 101)

    method <- "linear"
    zp <- interp2(x, y, Z, xp, yp, method)
    plot(xp, zp, type = "l", col = "blue")

    method = "nearest"
    zp <- interp2(x, y, Z, xp, yp, method)
    lines(xp, zp, col = "red")
    grid()
## End(Not run)

Matrix Inverse (Matlab Style)

Description

Invert a numeric or complex matrix.

Usage

inv(a)

Arguments

Details

Computes the matrix inverse by calling solve(a) and catching the error if the matrix is nearly singular.

Value

square matrix that is the inverse of a.

Note

inv() is the function name used in Matlab/Octave.

See Also

solve

Examples

A <- hilb(6)
B <- inv(A)
B
# Compute the inverse matrix through Cramer's rule:
n <- nrow(A)
detA <- det(A) 
b <- matrix(NA, nrow = n, ncol = n)
for (i in 1:n) {
    for (j in 1:n) {
        b[i, j] <- (-1)^(i+j) * det(A[-j, -i]) / detA
    }
}
b

Inverse Laplacian

Description

Numerical inversion of Laplace transforms.

Usage

invlap(Fs, t1, t2, nnt, a = 6, ns = 20, nd = 19)

Arguments

Details

The transform Fs may be any reasonable function of a variable s^a, where a is a real exponent. Thus, the function invlap can solve fractional problems and invert functions Fs containing (ir)rational or transcendental expressions.

Value

Returns a list with components x the x-coordinates and y the y-coordinates representing the original function in the interval [t1,t2].

Note

Based on a presentation in the first reference. The function invlap on MatlabCentral (by ) served as guide. The Talbot procedure from the second reference could be an interesting alternative.

References

J. Valsa and L. Brancik (1998). Approximate Formulae for Numerical Inversion of Laplace Transforms. Intern. Journal of Numerical Modelling: Electronic Networks, Devices and Fields, Vol. 11, (1998), pp. 153-166.

L.N.Trefethen, J.A.C.Weideman, and T.Schmelzer (2006). Talbot quadratures and rational approximations. BIT. Numerical Mathematics, 46(3):653–670.

Examples

Fs <- function(s) 1/(s^2 + 1)           # sine function
Li <- invlap(Fs, 0, 2*pi, 100)

## Not run: 
plot(Li[[1]], Li[[2]], type = "l", col = "blue"); grid()

Fs <- function(s) tanh(s)/s             # step function
L1 <- invlap(Fs, 0.01, 20, 1000)
plot(L1[[1]], L1[[2]], type = "l", col = "blue")
L2 <- invlap(Fs, 0.01, 20, 2000, 6, 280, 59)
lines(L2[[1]], L2[[2]], col="darkred"); grid()

Fs <- function(s) 1/(sqrt(s)*s)
L1 <- invlap(Fs, 0.01, 5, 200, 6, 40, 20)
plot(L1[[1]], L1[[2]], type = "l", col = "blue"); grid()

Fs <- function(s) 1/(s^2 - 1)           # hyperbolic sine function
Li <- invlap(Fs, 0, 2*pi, 100)
plot(Li[[1]], Li[[2]], type = "l", col = "blue"); grid()

Fs <- function(s) 1/s/(s + 1)           # exponential approach
Li <- invlap(Fs, 0, 2*pi, 100)
plot(Li[[1]], Li[[2]], type = "l", col = "blue"); grid()

gamma <- 0.577215664901532              # Euler-Mascheroni constant
Fs <- function(s) -1/s * (log(s)+gamma) # natural logarithm
Li <- invlap(Fs, 0, 2*pi, 100)
plot(Li[[1]], Li[[2]], type = "l", col = "blue"); grid()

Fs <- function(s) (20.5+3.7343*s^1.15)/(21.5+3.7343*s^1.15+0.8*s^2.2+0.5*s^0.9)/s
L1 <- invlap(Fs, 0.01, 5, 200, 6, 40, 20)
plot(L1[[1]], L1[[2]], type = "l", col = "blue")
grid()
## End(Not run)

isempty Property

Description

Determine if an object is empty.

Usage

isempty(x)

Arguments

Details

An empty object has length zero.

Value

TRUE if x has length 0; otherwise, FALSE.

Examples

isempty(c(0))            # FALSE
isempty(matrix(0, 1, 0)) # TRUE

Positive Definiteness

Description

Test for positive definiteness.

Usage

isposdef(A, psd = FALSE, tol = 1e-10)

Arguments

Details

Whether matrix A is positive definite will be determined by applying the Cholesky decomposition. The matrix must be symmetric.

With psd=TRUE the matrix will be tested for being semi-positive definite. If not positive definite, still a warning will be generated.

Value

Returns TRUE or FALSE.

Examples

A <- magic(5)
# isposdef(A)
## [1] FALSE
## Warning message:
## In isposdef(A) : Matrix 'A' is not symmetric.
## FALSE

A <- t(A) %*% A
isposdef(A)
## [1] TRUE

A[5, 5] <- 0
isposdef(A)
## [1] FALSE

isprime Property

Description

Vectorized version, returning for a vector or matrix of positive integers a vector of the same size containing 1 for the elements that are prime and 0 otherwise.

Usage

  isprime(x)

Arguments

Details

Given an array of positive integers returns an array of the same size of 0 and 1, where the i indicates a prime number in the same position.

Value

array of elements 0, 1 with 1 indicating prime numbers

See Also

factors, primes

Examples

  x <- matrix(1:10, nrow=10, ncol=10, byrow=TRUE)
  x * isprime(x)

  # Find first prime number octett:
  octett <- c(0, 2, 6, 8, 30, 32, 36, 38) - 19
  while (TRUE) {
      octett <- octett + 210
      if (all(as.logical(isprime(octett)))) {
          cat(octett, "\n", sep="  ")
          break
      }
  }

Iterative Methods

Description

Iterative solutions of systems of linear equations.

Usage

itersolve(A, b, x0 = NULL, nmax = 1000, tol = .Machine$double.eps^(0.5),
            method = c("Gauss-Seidel", "Jacobi", "Richardson"))

Arguments

Details

Iterative methods are based on splitting the matrix A=(P-A)-A with a so-called ‘preconditioner’ matrix P. The methods differ in how to choose this preconditioner.

Value

Returns a list with components x the solution, iter the number of iterations, and method the name of the method applied.

Note

Richardson's method allows to specify a ‘preconditioner’; this has not been implemented yet.

References

Quarteroni, A., and F. Saleri (2006). Scientific Computing with MATLAB and Octave. Springer-Verlag, Berlin Heidelberg.

See Also

qrSolve

Examples

N <- 10
A <- Diag(rep(3,N)) + Diag(rep(-2, N-1), k=-1) + Diag(rep(-1, N-1), k=1)
b <- A %*% rep(1, N)
x0 <- rep(0, N)

itersolve(A, b, tol = 1e-8, method = "Gauss-Seidel")
# [1]  1  1  1  1  1  1  1  1  1  1
# [1]  87
itersolve(A, b, x0 = 1:10, tol = 1e-8, method = "Jacobi")
# [1]  1  1  1  1  1  1  1  1  1  1
# [1]  177

Jacobian Matrix

Description

Jacobian matrix of a function R^n –> R^m .

Usage

jacobian(f, x0, heps = .Machine$double.eps^(1/3), ...)

Arguments

Details

Computes the derivative of each funktion f_j by variable x_i separately, taking the discrete step h.

Value

Numeric m-by-n matrix J where the entry J[j, i] is \frac{\partial f_j}{\partial x_i}, i.e. the derivatives of function f_j line up in row i for x_1, \ldots, x_n.

Note

Obviously, this function is not vectorized.

References

Quarteroni, A., R. Sacco, and F. Saleri (2007). Numerical Mathematics. Second Edition, Springer-Verlag, Berlin Heidelberg.

See Also

gradient

Examples

##  Example function from Quarteroni & Saleri
f <- function(x) c(x[1]^2 + x[2]^2 - 1, sin(pi*x[1]/2) + x[2]^3)
jf <- function(x) 
          matrix( c(2*x[1], pi/2 * cos(pi*x[1]/2), 2*x[2], 3*x[2]^2), 2, 2)
all.equal(jf(c(1,1)), jacobian(f, c(1,1)))
# TRUE

Interpolation by Kriging

Description

Simple and ordinary Kriging interpolation and interpolating function.

Usage

kriging(u, v, u0, type = c("ordinary", "simple"))

Arguments