Package 'pracma'

Title: Practical Numerical Math Functions
Description: Provides a large number of functions from numerical analysis and linear algebra, numerical optimization, differential equations, time series, plus some well-known special mathematical functions. Uses 'MATLAB' function names where appropriate to simplify porting.
Authors: Hans W. Borchers [aut, cre]
Maintainer: Hans W. Borchers <[email protected]>
License: GPL (>= 3)
Version: 2.4.4
Built: 2024-06-07 05:53:41 UTC
Source: CRAN

Help Index


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 <[email protected]>

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

f

function in the differential equation y=f(x,y)y' = f(x, y).

a, b

endpoints of the interval.

y0

starting values at point a.

n

the number of steps from a to b.

...

additional parameters to be passed to the function.

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

subs

vector or matrix of positive integers, used as indices for the result vector.

val

numerical vector.

sz

size of the resulting array.

func

function to be applied to a vector of numbers.

fillval

value used to fill the array when there are no indices pointing to that component.

a

numerical vector.

first

logical, shall the first or last element encountered be used.

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

b

vector of unique elements of a.

m

vector of indices such that b = a[m]

n

vector of indices such that a = b[n]

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

a, b

vectors of real or complex numbers of the same length (or scalars).

Details

The arithmetic-geometric mean is defined as the common limit of the two sequences an+1=(an+bn)/2a_{n+1} = (a_n + b_n)/2 and bn+1=(anbn)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:

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

where m=(ab)/(a+b)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

f

Function with a fixpoint.

x0

Starting value.

nmax

Maximum number of iterations.

tol

Relative tolerance.

...

Additional variables passed to f.

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

x, y

x/y-coordinates of (irregular) grid points defining the curve.

xi

x-coordinates of points where to interpolate.

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

l, k

Arrays.

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

A

numeric matrix with at least two columns.

f

factor or integer vector with nrow(A) elements.

style

character variable, only possible values ‘cart’ or ‘pol’.

scaled

logical; if true scales each column to have mean 0 and standard deviation 1 (not yet implemented).

npts

number of points to plot.

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

z

Vector or matrix of real or complex numbers

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

fn

nonlinear function to be minimized.

x0

starting vector.

tol

relative tolerance, to be used as stopping rule.

maxfeval

maximum number of function calls.

...

additional arguments to be passed to the function.

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:

xmin

minimum solution found.

fmin

value of f at minimum.

nfeval

number of function calls performed.

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

ts

a time series.

edim

the embedding dimension, as for chaotic time series; a preferred value is 2.

r

filter factor; work on heart rate variability has suggested setting r to be 0.2 times the standard deviation of the data.

elag

embedding lag; defaults to 1, more appropriately it should be set to the smallest lag at which the autocorrelation function of the time series is close to zero. (At the moment it cannot be changed by the user.)

tau

delay time for subsampling, similar to elag.

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

f

parametrization of a curve in n-dim. space.

a, b

begin and end of the parameter interval.

nmax

maximal number of iterations.

tol

relative tolerance requested.

...

additional arguments to be passed to the function.

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 <[email protected]>

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

A

a square n-by-n matrix.

q

a vector of length n.

m

an integer.

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

xi, yi

x- and y-coordinates of supporting nodes.

x

x-coordinates of interpolation points.

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

F

matrix representing values of a function in two dimensions.

xn, yn

x- and y-coordinates of supporting nodes.

xf, yf

x- and y-coordinates of an interpolating grid..

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

n

the index, a whole number greater or equal to 0.

x

real number or vector of real numbers; if missing, the Bernoulli numbers will be given, otherwise the polynomial.

Details

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

The Bernoulli polynomials can be directly defined as

Bn(x)=k=0n(nk)bnkxkB_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 Bn=Bn(0)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

f

function to be approximated by Bernstein polynomials.

k

integer between 0 and n, the k-th Bernstein polynomial of order n.

n

order of the Bernstein polynomial(s).

x

numeric scalar or vector where the Bernstein polynomials will be calculated.

Details

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

Bk,n(x)=(nk)xk(1x)nkB_{k,n}(x) = {{n}\choose{k}} x^k (1-x)^{n-k}

and form a basis for the vector space of polynomials of degree nn over the interval [0,1][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

fun

Function or its name as a string.

a, b

interval end points.

maxiter

maximum number of iterations; default 100.

tol

absolute tolerance; default eps^(1/2)

...

additional arguments passed to the function.

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

x

a positive or negative floating point number.

k

number of binary digits after the decimal point

pos_sign

logical; shall the '+' sign be included.

break0

logical; shall trailing zeros be included.

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

n

integer greater or equal to 0.

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

...

sequence of non-empty, numeric matrices

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

fun

function whose root is to be found.

a, b

left and right end points of an interval; function values need to be of different sign at the endpoints.

maxiter

maximum number of iterations.

tol

relative tolerance.

...

additional arguments to be passed to the function.

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

root

location of the root.

f.root

funtion value at the root.

f.calls

number of function calls.

estim.prec

estimated relative precision.

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

Ffun

n functions of n variables.

x0

Numeric vector of length n.

J0

Jacobian of the function at x0.

...

additional parameters passed to the function.

maxiter

Maximum number of iterations.

tol

Tolerance, relative accuracy.

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

func

function with two or more input parameters.

x, y

two vectors, matrices, or arrays of the same size.

...

list of arrays of the same size.

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

f

function describing the differential equation y=f(t,y)y' = f(t, y).

t

vector of x-values where the values of the ODE function will be computed; needs to be increasingly sorted.

y0

starting values as column vector.

...

additional parameters to be passed to the function.

tol

relative tolerance in the Ricardson extrapolation.

t0, tfinal

start and end point of the interval.

kmax

maximal number of steps in the Richardson extrapolation.

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

f, g, h

functions on the right side of the differential equation. If f, g or h is a scalar instead of a function, it is assumed to be a constant coefficient in the differential equation.

x

x[1], x[2] are the interval borders where the solution shall be computed.

y

boundary conditions such that y(x[1]) = y[1], y(x[2]) = y[2].

n

number of intermediate grid points; default 50.

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)y'' = f(x) y' + g(x) y + h(x)

with the finite element method. The solution y(x)y(x) shall exist on the interval [a,b][a, b] with boundary conditions y(a)=yay(a) = y_a and y(b)=yby(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

xyz

cartesian coordinates x, y, z as vector or matrix.

tpr

spherical coordinates theta, phi, and r as vector or matrix.

prz

polar coordinates phi, r or cylindrical coordinates phi, r, z as vector or matrix.

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

dname

(relative or absolute) directory path.

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

n

a numeric vector or matrix

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

a

quadratic matrix; size should not be much larger than 100.

info

logical; if true, the inverse matrix will also be reported.

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

x

Numeric vector of points within interval [a, b].

fun

Function to be approximated.

a, b

Endpoints of the interval.

n

An integer >= 0.

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

fun

function to be approximated.

a, b

endpoints of the interval.

n

an integer >= 0.

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

n

an integer >= 0.

x

a numeric vector, possibly empty; default NULL.

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

xp, yp

Vectors representing the x and y coordinates of plane points

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

lst

Character vector of names of variables in the global environment.

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

f

function, the integrand, without singularities.

a, b

lower and upper limit of the integral; must be finite.

n

Number of Chebyshev nodes to account for.

...

Additional parameters to be passed to the function

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

a

numeric vector of some length n

m

integer with 0 <= m <= n

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

p

vector representing a polynomial

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

f

Function that is to be differentiated.

x0

Point at which to differentiate the function.

h

Step size to be applied; shall be very small.

...

Additional variables to be passed to f.

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(x0)f'(x_0) of ff at x0x_0. The function is vectorized in x0.

Note

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

Author(s)

HwB <[email protected]>

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

M

Numeric matrix; vectors will be considered as column vectors.

p

Indicates the p-norm. At the moment, norms other than p=2 are not implemented.

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

x, y

real or complex vectors.

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

z

numeric or complex scalar or vector.

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

f

the integrand as function of two variables.

a, b

lower and upper limit of the integral.

n

number of subintervals (grid points).

nodes

number of nodes in the Newton-Cotes formula.

...

additional parameters to be passed to the function.

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

z

numeric or complex scalar or vector.

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

f

function in the differential equation y=f(x,y)y' = f(x, y);
defined as a function R×RmRmR \times R^m \rightarrow R^m, where mm is the number of equations.

t0, t1

start and end points of the interval.

y0

starting values as row or column vector; for mm equations y0 needs to be a vector of length m.

N

number of steps.

...

Additional parameters to be passed to the function.

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

x

numeric vector or matrix

y

numeric vector or matrix

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

A

matrix of size (n-1) x n where n >= 2.

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 R3R^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

x, y

x- and y-coordinates of points to be interpolated.

xi

x-coordinates of points at which the interpolation is to be performed.

endp2nd

logical; if true, the derivatives at the endpoints are prescribed by der.

der

a two-components vector prescribing derivatives at endpoints.

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

u

the parameter vector.

x, y

x-, y-coordinates for each parameter value.

n

order of the polynomials, the same in x- and y-dirction.

U

parameter values where points will be fixed.

V

matrix with two columns and lemgth(U) rows; first column contains the x-, the second the y-values of those points kept fixed.

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

x

vector of real values.

nmax

the maximum number of cutting points to choose

quant

quantile of the gaps to consider for cuts.

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

f

function of two variables, the integrand.

xa, xb

left and right endpoint for first variable.

ya, yb

left and right endpoint for second variable.

za, zb

left and right endpoint for third variable.

dim

which variable to integrate first.

subdivs

number of subdivisions to use.

tol

relative tolerance to use in integrate.

...

additional parameters to be passed to the integrand.

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

b, a

real or complex vectors.

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

x

vector of (time) points at which the differential equation has been solved.

y

values of the function(s) that have been computed for the given (time) points.

yv

point or numeric vector at which the solution is wanted.

idx

index of functions whose vales shall be returned.

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

deg

(array of) angles in degrees.

rad

(array of) angles in radians.

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

x

vector or matrix, columns considered as the time series.

tt

trend type, ‘constant’ or ‘linear’, default is ‘linear’.

bp

break points, indices between 1 and nrow(x).

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

x

vector of (time) points at which the differential equation has been solved.

y

values of the function(s) that have been computed for the given (time) points.

xp

point or numeric vector at which the solution is wanted; must be sorted.

idx

index of functions whose vales shall be returned.

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

x

vector or matrix

k

integer indicating a secondary diagonal

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

...

any R object that can be printed.

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

X

matrix of some size m x k; vector will be taken as row matrix.

Y

matrix of some size n x k; vector will be taken as row matrix.

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

x

numeric vector or matrix

y

numeric vector or matrix

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

a

real or complex square matrix

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

A

a real symmetric matrix.

tol

requested tolerance.

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

x

numeric or complex vector.

d

parameter to select one of the Einstein functions E1, E2, E3, E4.

Details

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

The functions are defined as:

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

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

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

E4(x)=xex1ln(1ex)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

u

numeric vector.

m

input vector, all input elements must satisfy 0 <= x <= 1.

tol

tolerance; default is machine precision.

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, snsn is the inverse function for

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

with sn(u)=sin(ϕ)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

x

scalar or numerical vector or matrix.

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

x, y

vector of real numbers.

z

real or complex number; must be a scalar.

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

x, y

x-, y-coordinates

xerr, yerr

length of the error bars, relative to the x-, y-values.

bar.col

color of the error bars; default: red

bar.len

length of the cross bars orthogonal to the error bars; default: 0.01.

grid

logical; should the grid be plotted?; default: true

with

logical; whether to end the error bars with small cross bars.

add

logical; should the error bars be added to an existing plot?; default: false.

...

additional plotting parameters that will be passed to the plot function.

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

z

Real or complex number or a numeric or complex vector.

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

f

function in the differential equation y=f(x,y)y' = f(x, y).

a, b

start and end points of the interval.

y0

starting value at a.

n

number of grid points.

improved

logical; shall the Heun method be used; default TRUE.

...

additional parameters to be passed to the function.

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>0x > 0 as

xettdt\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

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

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

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

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

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

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

The integral LiLi approximates the prime number function π(n)\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

x

vector of real or complex numbers.

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

A

numeric square matrix.

np

number of points to use on the unit circle.

Details

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

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

where GG is a closed contour around the eigenvalues of AA.

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

m, n

numeric scalars specifying size of the matrix

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

f

2-D function to be plotted, must accept (x,y) as a vector.

xlim, ylim

defines x- and y-ranges as intervals.

n

number of grid points in each direction.

col

colour of isolines lines, resp. the surface color.

filled

logical; shall the contour plot be

...

parameters to be passed to the persp function.

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

f

Function to be plotted.

interval

interval [a, b] to plot the function in

a, b

Left and right endpoint for the plot.

n

Number of points to plot.

col

Color of the function graph.

add

logical; shall the polt be added to an existing plot.

lty

line type; default 1.

lwd

line width; default 1.

marker

no. of markers to be added to the curve; defailt: none.

pch

poimt character; default circle.

grid

Logical; shall a grid be plotted?; default TRUE.

gridcol

Color of grid points.

fill

Logical; shall the area between function and axis be filled?; default: FALSE.

fillcol

Color of fill area.

xlab

Label on the x-axis.

ylab

Label on the y-axis.

main

Title of the plot

...

More parameters to be passed to plot.

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

fun

function to be plotted.

interv

left and right endpoint for the plot.

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

n

Vector of integers, for fact, resp. a single integer for factorial2.

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

n

nonnegative integer

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

f

function to be differentiated.

x

point(s) where differentiation will take place.

n

order of derivative, should only be between 1 and 8; for n=0 function values will be returned.

h

step size: if h=0 step size will be set automatically.

method

one of “central”, “forward”, or “backward”.

...

more variables to be passed to function f.

Details

Derivatives are computed applying central difference formulas that stem from the Taylor series approximation. These formulas have a convergence rate of O(h2)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

f

Function or its name as a string.

a, b

endpoints of the interval

endp

logical; shall the endpoints be considered as possible minima?

tol

absolute tolerance; default eps^(1/2).

...

Additional arguments to be passed to f.

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

figno

(single) number of plot device.

title

title of the plot device; not yet used.

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

x

single number.

xs

numeric vector, not necessarily sorted.

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

f

functions whose minima shall be found.

a, b

endpoints of the interval.

n

number of subintervals to generate and search.

tol

has no effect at this moment.

...

Additional parameters to be passed to the function.

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

x

numerical vector taken as a time series (no NAs allowed)

nups

minimum number of increasing steps before a peak is reached

ndowns

minimum number of decreasing steps after the peak

zero

can be ‘+’, ‘-’, or ‘0’; how to interprete succeeding steps of the same value: increasing, decreasing, or special

peakpat

define a peak as a regular pattern, such as the default pattern [+]{1,}[-]{1,}; if a pattern is provided, parameters nups and ndowns are not taken into account

minpeakheight

the minimum (absolute) height a peak has to have to be recognized as such

minpeakdistance

the minimum distance (in indices) peaks have to have to be counted

threshold

the minimum

npeaks

the number of peaks to return

sortstr

logical; should the peaks be returned sorted in decreasing oreder of their maximum value

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

v

logical or numeric vector or array

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

f

functions whose roots shall be found.

a, b

endpoints of the interval.

n

number of subintervals to generate and search.

tol

tolerance for identifying zeros.

...

Additional parameters to be passed to the function.

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

x0

start value.

f

function to be minimized.

g

gradient function of f; if NULL, a numerical gradient will be calculated.

maxiter

max. number of iterations.

tol

relative tolerance, to be used as stopping rule.

Details

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

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

where HH is the Hessian and gradgrad 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:

xmin

minimum solution found.

fmin

value of f at minimum.

niter

number of iterations performed.

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

a

numeric or complex matrix

dim

flipping dimension; can only be 1 (default) or 2

sz

integer vector of length 1 or 2.

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

f

function whose minimum or maximum is to be found.

a, b

endpoints of the interval to be searched.

maxiter

maximal number of iterations.

maximum

logical; shall maximum or minimum be found; default FALSE.

tol

relative tolerance; left over for compatibility.

rel.tol, abs.tol

relative and absolute tolerance.

...

additional variables to be passed to the function.

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

xmin

location of the minimum resp. maximum.

fmin

function value at the optimum.

niter

number of iterations used.

estim.prec

estimated precision.

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

x0

starting point.

fn

objective function to be minimized.

gr

gradient function of the objective; not used for SQP method.

...

additional parameters to be passed to the function.

method

method options 'SQP', 'auglag'; only 'SQP is implemented.

A, b

linear ineqality constraints of the form A x <= b .

Aeq, beq

linear eqality constraints of the form Aeq x = beq .

lb, ub

bounds constraints of the form lb <= x <= ub .

hin

nonlinear inequality constraints of the form hin(x) <= 0 .

heq

nonlinear equality constraints of the form heq(x) = 0 .

tol

relative tolerance.

maxiter

maximum number of iterations.

maxfeval

maximum number of function evaluations.

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:

par

the best minimum found.

value

function value at the minimum.

convergence

integer indicating the terminating situation.

info

parameter list describing the final situation.

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

fn

function whose minimum or maximum is to be found.

x0

point considered near to the optimum.

...

additional variables to be passed to the function.

lower, upper

lower and upper bounds constraints.

method

"Nelder-Mead" (default) or "Hooke-Jeeves"; can be abbreviated.

minimize

logical; shall a minimum or a maximum be found.

maxiter

maximal number of iterations

tol

relative tolerance.

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

xopt

location of the location of minimum resp. maximum.

fmin

function value at the optimum.

count

number of function calls.

convergence

info about convergence: not used at the moment.

info

special information from the solver.

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

x0

starting point.

fn

objective function to be minimized.

gr

gradient function of the objective.

...

additional parameters to be passed to the function.

tol

relative tolerance.

maxiter

maximum number of iterations.

maxfeval

maximum number of function evaluations.

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:

par

the best minimum found.

value

function value at the minimum.

counts

number of function and gradient calls.

convergence

integer indicating the terminating situation.

message

description of the final situation.

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

f, g

functions given by name or string.

x1, x2

endpoints of the interval.

p

Numeric scalar or Inf, -Inf; default is 2.

npoints

number of points to be considered in the interval.

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

x

grid points on the x-axis, must be distinct.

y

discrete values of the function at the grid points.

xs

point at which to approximate (not vectorized).

k

order of derivative, k<=length(x)-1 required.

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

fmt

a character vector of format strings.

...

values passed to the format string.

file

a connection or a character string naming the file to print to; default is "" which means standard output.

append

logical; shall the output be appended to the file; default is FALSE.

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

n

integer, the ‘order’ of the curve

which

character string, which curve to cumpute.

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

x

numeric vector.

Details

The normalized Fresnel integrals are defined as

S(x)=0xsin(π/2t2)dtS(x) = \int_0^x \sin(\pi/2 \, t^2) dt

C(x)=0xcos(π/2t2)dtC(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

f

function describing the system of equations.

x0

point near to the root.

J

Jacobian function of f, or NULL.

maxiter

maximum number of iterations in gaussNewton.

tol

tolerance to be used in Gauss-Newton.

...

additional variables to be passed to the function.

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

x

location of the solution.

fval

function value at the solution.

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

fun

function whose root is sought.

x

a point near the root or an interval giving end points.

maxiter

maximum number of iterations.

tol

relative tolerance.

...

additional arguments to be passed to the function.

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

x

location of the root.

fval

function value at the root.

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

fz

complex(-analytic) function.

z0

complex point near the assumed root.

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

x

positive real number.

a

real number.

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

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

and

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

while the regularized incomplete gamma function is γ(x,a)/Γ(a)\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

z

Real or complex number or a numeric or complex vector.

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

f

function to be integrated.

a, b

end points of the interval.

...

variables to be passed to the function.

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

n

Number of nodes in the interval ]-Inf, Inf[.

Details

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

f(x)ex2dx\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

n

Number of nodes in the interval [0, Inf[.

a

exponent of x in the integrand: must be greater or equal to 0, otherwise the integral would not converge.

Details

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

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

over the infinite interval ]0,[]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

n

Number of nodes in the interval [a,b].

a, b

lower and upper limit of the integral; must be finite.

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 f1(x)2++fm(x)2f_1(x)^2 + \ldots + f_m(x)^2 or FFF' F where F=(f1,,fm)F = (f_1, \ldots, f_m) is a multivariate function of nn variables, not necessarily n=mn = m.

Usage

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

Arguments

Ffun

m functions of n variables.

Jfun

function returning the Jacobian matrix of Ffun; if NULL, the default, the Jacobian will be computed numerically. The gradient of f will be computed internally from the Jacobian (i.e., cannot be supplied).

x0

Numeric vector of length n.

maxiter

Maximum number of iterations.

tol

Tolerance, relative accuracy.

...

Additional parameters to be passed to f.

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

a, b

vectors of integers.

extended

logical; if TRUE the extended Euclidean algorithm will be applied.

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

P

matrix of points, x_i-coordinates in the ith column.

tol

relative tolerance.

maxiter

maximum number of iterations.

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

x

numeric vector, matrix, or array.

dim

dimension along which to take the mean; dim=1 means along columns, dim=2 along rows, the result will still be a row vector, not a column vector as in Matlab.

percent

percentage, between 0 and 100, of trimmed values.

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

A

numeric square matrix.

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

A

square matrix.

b

numerical vector or column vector.

x0

initial iterate.

errtol

relative residual reduction factor.

kmax

maximum number of iterations

reorth

reorthogonalization method, see Details.

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

f

Function or its name as a string.

a, b

endpoints of the interval.

maxiter

maximum number of iterations.

tol

absolute tolerance; default sqrt(eps).

...

Additional arguments to be passed to f.

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

f

function of several variables.

x0

point where the gradient is to build.

heps

step size.

...

more variables to be passed to function f.

Details

Computes the gradient

(fx1,,fxn)(\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

F

vector of function values, or a matrix of values of a function of two variables.

h1

x-coordinates of grid points, or one value for the difference between grid points in x-direction.

h2

y-coordinates of grid points, or one value for the difference between grid points in y-direction.

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:

X

numerical gradient/slope in x-direction.

Y

numerical gradient/slope in x-direction.

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

A

numeric matrix with nrow(A)>=ncol(A).

tol

numerical tolerance for being equal to zero.

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

n

An integer of the form 2^e, 12*2^e, or 20*2^e

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

fun

function whose root is to be found.

x0

starting value for the iteration.

maxiter

maximum number of iterations.

tol

absolute tolerance; default eps^(1/2)

...

additional arguments to be passed to the function.

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

x

numeric vector representing a time series

k

window length 2*k+1 in indices

t0

threshold, default is 3 (Pearson's rule), see below.

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

a

vector that will be the first column

b

vector that if present will form the last row.

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

P, Q

numerical matrices, representing points in an m-dim. space.

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

loc1, loc2

Locations on earth; for format see Details.

R

Average earth radius R = 6371 km, can be changed on input.

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

A

square matrix

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,

H

the upper Hessenberg Form of matrix A.

H

a unitary matrix.

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

f

univariate function of several variables.

x0

point in RnR^n.

h

step size.

...

variables to be passed to f.

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 2fxixj\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

f

function whose hessian is to be computed.

x

point in R^n.

v

vector of length n.

csd

logocal, shall complex-step be applied.

...

more arguments to be passed to the function.

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

n

positive integer specifying the dimension of the Hilbert matrix

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

x

numeric vector or matrix.

edges

numeric vector of grid points, must be monotonically non-decreasing.

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

x

numeric vector or matrix.

n

maximum number of bins.

plotting

logical; shall a histogram be plotted.

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

x0

starting vector.

fn

nonlinear function to be minimized.

...

additional arguments to be passed to the function.

lb, ub

lower and upper bounds.

tol

relative tolerance, to be used as stopping rule.

maxfeval

maximum number of allowed function evaluations.

target

iteration stops when this value is reached.

info

logical, whether to print information during the main loop.

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:

xmin

minimum solution found so far.

fmin

value of f at minimum.

count

number of function evaluations.

convergence

NOT USED at the moment.

info

special info from the solver.

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

p

Numeric vector representing a polynomial.

x

Numeric scalar, vector or matrix at which to evaluate the polynomial.

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

A

numeric matrix with nrow(A)>=ncol(A).

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

x

numeric scalar or vector.

n

positive integer.

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)=sin(πt)πtsinc(t) = \frac{\sin(\pi t)}{\pi t}. It is the continuous inverse Fourier transform of the rectangular pulse of width 2π2\pi and height 11.

psinc is the 'periodic sinc function' and is defined as psinc(x,n)=sin(xn/2)nsin(x/2)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

x

a time series.

d

smallest box size; default 50.

display

logical; shall the results be printed to the console?

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

x, y

Vectors of real or complex numbers of the same size

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

x

a real or complex vector

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

x, y

x-, y-coordinates of points to be tested for being inside the polygon region.

xp, yp

coordinates of the vertices specifying the polygon.

boundary

Logical; does the boundary belong to the interior.

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

fun

integrand, univariate (vectorized) function.

xmin, xmax

endpoints of the integration interval.

method

integration procedure, see below.

no_intervals

number of subdivisions at at start.

random

logical; shall the length of subdivisions be random.

reltol

relative tolerance.

abstol

absolute tolerance; not used.

...

additional parameters to be passed to the function.

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

fun

function

xmin, xmax

lower and upper limits of x.

ymin, ymax

lower and upper limits of y.

zmin, zmax

lower and upper limits of z.

sector

logical.

reltol

relative tolerance.

abstol

absolute tolerance.

maxlist

maximum length of the list of rectangles.

singular

logical; are there singularities at vertices.

vectorized

logical; is the function fully vectorized.

...

additional parameters to be passed to the function.

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

x

Numeric vector; points on the x-axis; at least two points require; will be sorted if necessary.

y

Numeric vector; values of the assumed underlying function; x and y must be of the same length.

xi

Numeric vector; points at which to compute the interpolation; all points must lie between min(x) and max(x).

method

One of “constant", “linear", “nearest", “spline", or “cubic"; default is “linear"

Details

Interpolation to find yi, the values of the underlying function at the points in the vector xi.

Methods can be:

linear linear interpolation (default)
constant constant between points
nearest nearest neighbor interpolation
spline cubic spline interpolation
cubic cubic Hermite interpolation

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

x, y

vectors with monotonically increasing elements, representing x- and y-coordinates of the data values in Z.

Z

numeric length(y)-by-length(x) matrix.

xp, yp

x-, y-coordinates of points at which interpolated values will be computed.

method

interpolation method, “linear” the most useful.

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

a

real or complex square matrix

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

Fs

function representing the function to be inverse-transformed.

t1, t2

end points of the interval.

nnt

number of grid points between t1 and t2.

a

shift parameter; it is recommended to preserve value 6.

ns, nd

further parameters, increasing them leads to lower error.

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

x

an R object

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

A

symmetric matrix

psd

logical, shall semi-positive definiteness be tested?

tol

tolerance to check symmetry and Cholesky decomposition.

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

x

vector or matrix of nonnegative integers

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

A

numerical matrix, square and non-singular.

b

numerical vector or column vector.

x0

starting solution for iteration; defaults to null vector.

nmax

maximum number of iterations.

tol

relative tolerance.

method

iterative method, Gauss-Seidel, Jacobi, or Richardson.

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

f

m functions of n variables.

x0

Numeric vector of length n.

heps

This is h in the derivative formula.

...

parameters to be passed to f.

Details

Computes the derivative of each funktion fjf_j by variable xix_i separately, taking the discrete step hh.

Value

Numeric m-by-n matrix J where the entry J[j, i] is fjxi\frac{\partial f_j}{\partial x_i}, i.e. the derivatives of function fjf_j line up in row ii for x1,,xnx_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

u

an nxm-matrix of n points in the m-dimensional space.

v

an n-dim. (column) vector of interpolation values.

u0

a kxm-matrix of k points in R^m to be interpolated.

type

character; values ‘simple’ or ‘ordinary’; no partial matching.

Details

Kriging is a geo-spatial estimation procedure that estimates points based on the variations of known points in a non-regular grid. It is especially suited for surfaces.

Value

kriging returns a k-dim. vektor of interpolation values.

Note

In the literature, different versions and extensions are discussed.

References

Press, W. H., A. A. Teukolsky, W. T. Vetterling, and B. P. Flannery (2007). Numerical recipes: The Art of Scientific Computing (3rd Ed.). Cambridge University Press, New York, Sect. 3.7.4, pp. 144-147.

See Also

akimaInterp, barylag2d, package kriging

Examples

##  Interpolate the Saddle Point function
f <- function(x) x[1]^2 - x[2]^2       # saddle point function

set.seed(8237)
n <- 36
x <- c(1, 1, -1, -1, runif(n-4, -1, 1)) # add four vertices
y <- c(1, -1, 1, -1, runif(n-4, -1, 1))
u <- cbind(x, y)
v <- numeric(n)
for (i in 1:n) v[i] <- f(c(x[i], y[i]))

kriging(u, v, c(0, 0))                      #=>  0.006177183
kriging(u, v, c(0, 0), type = "simple")     #=>  0.006229557

## Not run: 
xs <- linspace(-1, 1, 101)              # interpolation on a diagonal
u0 <- cbind(xs, xs)

yo <- kriging(u, v, u0, type = "ordinary")  # ordinary kriging
ys <- kriging(u, v, u0, type = "simple")    # simple kriging
plot(xs, ys, type = "l", col = "blue", ylim = c(-0.1, 0.1),
             main = "Kriging interpolation along the diagonal")
lines(xs, yo, col = "red")
legend( -1.0, 0.10, c("simple kriging", "ordinary kriging", "function"),
        lty = c(1, 1, 1), lwd = c(1, 1, 2), col=c("blue", "red", "black"))
grid()
lines(c(-1, 1), c(0, 0), lwd