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 wellknown 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:  20240607 05:53:41 UTC 
Source:  CRAN 
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.
The package encompasses functions from all areas of numerical analysis, for example:
Root finding and minimization of univariate functions,
e.g. NewtonRaphson, BrentDekker, 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 GaussKronrod quadrature.
Solvers for ordinary differential equations and systems,
EulerHeun, classical RungeKutta, ode23, or predictorcorrector method
such as the AdamsBashfordMoulton.
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 fullblown 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.
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)
Hans Werner Borchers
Maintainer: Hans W Borchers <[email protected]>
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. SpringerVerlag, 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). EditorinChief: 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, SpringerVerlag, Berlin Heidelberg.
Skiena, St. S. (2008). The Algorithm Design Manual. Second Edition, SpringerVerlag, London. The Stony Brook Algorithm Repository: https://algorist.com/algorist.html.
Stoer, J., and R. Bulirsch (2002). Introduction to Numerical Analysis. Third Edition, SpringerVerlag, New York.
Strang, G. (2007). Computational Science and Engineering. WellesleyCambridge 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.
The R package ‘matlab’ contains some of the basic routines from Matlab, but unfortunately not any of the higher math routines.
## Not run:
## See examples in the help files for all functions.
## End(Not run)
Thirdorder AdamsBashfordMoulton predictorcorrector method.
abm3pc(f, a, b, y0, n = 50, ...)
f 
function in the differential equation 
a , b

endpoints of the interval. 
y0 
starting values at point 
n 
the number of steps from 
... 
additional parameters to be passed to the function. 
Combined AdamsBashford and AdamsMoulton (or: multistep) method of third order with corrections according to the predictorcorrector approach.
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.
This function serves demonstration purposes only.
Fausett, L. V. (2007). Applied Numerical Analysis Using Matlab. Second edition, Prentice Hall.
## Attempt on a nonstiff 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)
accumarray
groups elements from a data set and applies a function
to each group.
accumarray(subs, val, sz = NULL, func = sum, fillval = 0)
uniq(a, first = FALSE)
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. 
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]
.
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 
n 
vector of indices such that 
The Matlab function accumarray
can also handle sparse matrices.
## 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
The arithmeticgeometric mean of real or complex numbers.
agmean(a, b)
a , b

vectors of real or complex numbers of the same length (or scalars). 
The arithmeticgeometric mean is defined as the common limit of the two
sequences $a_{n+1} = (a_n + b_n)/2$
and $b_{n+1} = \sqrt(a_n b_n)$
.
When used for negative or complex numbers, the complex square root function is applied.
Returns a list with compoinents: agm
a vector of arithmeticgeometric
means, componentwise, niter
the number of iterations, and prec
the overall estimated precision.
Gauss discovered that elliptic integrals can be effectively computed via the arithmeticgeometric mean (see example below), for example:
$\int_0^{\pi/2} \frac{dt}{\sqrt{1  m^2 sin^2(t)}} = \frac{(a+b) \pi}{4 \cdot agm(a,b)}$
where $m = (ab)/(a+b)$
https://mathworld.wolfram.com/ArithmeticGeometricMean.html
Arithmetic, geometric, and harmonic mean.
## Accuracy test: Gauss constant
1/agmean(1, sqrt(2))$agm  0.834626841674073186 # 1.11e16 < eps = 2.22e16
## Gauss' AGMbased computation of \pi
a < 1.0
b < 1.0/sqrt(2)
s < 0.5
d < 1L
while (abs(ab) > eps()) {
t < a
a < (a + b)*0.5
b < sqrt(t*b)
c < (at)*(at)
d < 2L * d
s < s  d*c
}
approx_pi < (a+b)^2 / s / 2.0
abs(approx_pi  pi) # 8.881784e16 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*N1)
a < 1
b < a * (1m) / (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's acceleration method.
aitken(f, x0, nmax = 12, tol = 1e8, ...)
f 
Function with a fixpoint. 
x0 
Starting value. 
nmax 
Maximum number of iterations. 
tol 
Relative tolerance. 
... 
Additional variables passed to f. 
Aitken's acceleration method, or deltasquared 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.
The fixpoint (as found so far).
Sometimes used to accerate NewtonRaphson (Steffensen's method).
Quarteroni, A., and F. Saleri (2006). Scientific Computing with Matlab and Octave. Second Edition, SpringerVerlag, Berlin Heidelberg.
# 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
Interpolate smooth curve through given points on a plane.
akimaInterp(x, y, xi)
x , y

x/ycoordinates of (irregular) grid points defining the curve. 
xi 
xcoordinates of points where to interpolate. 
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.
Returns the interpolated values at the points xi
as a vector.
There is also a 2dimensional version in package ‘akima’.
Matlab code by H. Shamsundar under BSC License; reimplementation in R by Hans W Borchers.
Akima, H. (1970). A New Method of Interpolation and Smooth Curve Fitting Based on Local Procedures. Journal of the ACM, Vol. 17(4), pp 589602.
Hyman, J. (1983). Accurate Monotonicity Preserving Cubic Interpolation. SIAM J. Sci. Stat. Comput., Vol. 4(4), pp. 645654.
Akima, H. (1996). Algorithm 760: RectangularGridData Surface Fitting that Has the Accurancy of a Bicubic Polynomial. ACM TOMS Vol. 22(3), pp. 357361.
Akima, H. (1996). Algorithm 761: ScatteredData Surface Fitting that Has the Accuracy of a Cubic Polynomial. ACM TOMS, Vol. 22(3), pp. 362371.
kriging
, akima::aspline
, akima::interp
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)
and(l, k)
resp. or(l, k)
the same as (l & k) + 0
resp.
(l  k) + 0
.
and(l, k)
or(l, k)
l , k

Arrays. 
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.
Logical vector.
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)
Plots Andrews' curves in cartesian or polar coordinates.
andrewsplot(A, f, style = "pol", scaled = FALSE, npts = 101)
A 
numeric matrix with at least two columns. 
f 
factor or integer vector with 
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. 
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.
Generates a plot, no return value.
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.
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.
polar
, andrews::andrews
## 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 (Matlab style)
Real(z)
Imag(z)
angle(z)
z 
Vector or matrix of real or complex numbers 
These are just Matlab names for the corresponding functions in R. The
angle
function is simply defined as atan2(Im(z), Re(z))
.
returning real or complex values; angle
returns in radians.
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.
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)
An implementation of the NelderMead algorithm for derivativefree optimization / function minimization.
anms(fn, x0, ...,
tol = 1e10, maxfeval = NULL)
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. 
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 NelderMead. Of course, if the
optimum is near to the boundary, results will not be as accurate as when the
minimum is in the interior.
List with following components:
xmin 
minimum solution found. 
fmin 
value of 
nfeval 
number of function calls performed. 
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 NelderMead see the package ‘dfoptim’.
Nelder, J., and R. Mead (1965). A simplex method for function minimization. Computer Journal, Volume 7, pp. 308313.
O'Neill, R. (1971). Algorithm AS 47: Function Minimization Using a Simplex Procedure. Applied Statistics, Volume 20(3), pp. 338345.
J. C. Lagarias et al. (1998). Convergence properties of the NelderMead simplex method in low dimensions. SIAM Journal for Optimization, Vol. 9, No. 1, pp 112147.
Fuchang Gao and Lixing Han (2012). Implementing the NelderMead simplex algorithm with adaptive parameters. Computational Optimization and Applications, Vol. 51, No. 1, pp. 259277.
## Rosenbrock function
rosenbrock < function(x) {
n < length(x)
x1 < x[2:n]
x2 < x[1:(n1)]
sum(100*(x1x2^2)^2 + (1x2)^2)
}
anms(rosenbrock, c(0,0,0,0,0))
# $xmin
# [1] 1 1 1 1 1
# $fmin
# [1] 8.268732e21
# $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)
Calculates the approximate or sample entropy of a time series.
approx_entropy(ts, edim = 2, r = 0.2*sd(ts), elag = 1)
sample_entropy(ts, edim = 2, r = 0.2*sd(ts), tau = 1)
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 
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 selfmatching, and it does not depend that much on the length of the time series.
The approximate, or sample, entropy, a scalar value.
This code here derives from Matlab versions at Mathwork's File Exchange, “Fast Approximate Entropy” and “Sample Entropy” by Kijoon Lee under BSD license.
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.
RHRV::CalculateApEn
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)
Calculates the arc length of a parametrized curve.
arclength(f, a, b, nmax = 20, tol = 1e05, ...)
f 
parametrization of a curve in ndim. 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. 
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 onedimensional
function f:R>R
by defining F
(if f
is vectorized)
as F:t>c(t,f(t))
.
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.
If by chance certain equidistant points of the curve lie on a straight line,
the result may be wrong, then use polylength
below.
HwB <[email protected]>
## Example: parametrized 3Dcurve 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 = 1e10) #=> 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 1dimensional 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 = 1e12, 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, (i1)*0.5, i*0.5, tol = 1e10)$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 = 1e14)
print(L1, digits = 16)
# [1] 82.81020372882216
## End(Not run)
## Not run:
# 
# Arclength 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 generates an orthonormal basis of the Krylov space and a Hessenberg matrix.
arnoldi(A, q, m)
A 
a square nbyn matrix. 
q 
a vector of length n. 
m 
an integer. 
arnoldi(A, q, m)
carries out m
iterations of the
Arnoldi iteration with nbyn matrix A
and starting vector
q
(which need not have unit 2norm). For m < n
it
produces an nby(m+1) matrix Q
with orthonormal columns
and an (m+1)bym 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 mth column of the mbym identity matrix.
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)
.
Nicholas J. Higham (2008). Functions of Matrices: Theory and Computation, SIAM, Philadelphia.
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 in one dimension.
barylag(xi, yi, x)
xi , yi

x and ycoordinates of supporting nodes. 
x 
xcoordinates of interpolation points. 
barylag
interpolates the given data using the barycentric
Lagrange interpolation formula (vectorized to remove all loops).
Values of interpolated data at points x
.
Barycentric interpolation is preferred because of its numerical stability.
Berrut, J.P., and L. Nick Trefethen (2004). “Barycentric Lagrange Interpolation”. SIAM Review, Vol. 46(3), pp.501–517.
Lagrange or Newton interpolation.
## 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)
Twodimensional barycentric Lagrange interpolation.
barylag2d(F, xn, yn, xf, yf)
F 
matrix representing values of a function in two dimensions. 
xn , yn

x and ycoordinates of supporting nodes. 
xf , yf

x and ycoordinates of an interpolating grid.. 
Wellknown Lagrange interpolation using barycentric coordinates, here extended to two dimensions. The function is completely vectorized.
xcoordinates run downwards in F, ycoordinates to the right. That conforms to the usage in image or contour plots, see the example below.
Matrix of size length(xf)
bylength(yf)
giving the interpolated
values at al the grid points (xf, yf)
.
Copyright (c) 2004 Greg von Winckel of a Matlab function under BSD license; translation to R by Hans W Borchers with permission.
Berrut, J.P., and L. Nick Trefethen (2004). “Barycentric Lagrange Interpolation”. SIAM Review, Vol. 46(3), pp.501–517.
## Example from Rhelp
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) # "NelderMead"
# $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)
The Bernoulli numbers are a sequence of rational numbers that play an important role for the series expansion of hyperbolic functions, in the EulerMacLaurin formula, or for certain values of Riemann's function at negative integers.
bernoulli(n, x)
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. 
The calculation of the Bernoulli numbers uses the values of the zeta function
at negative integers, i.e. $B_n = n \, zeta(1n)$
. Bernoulli numbers
$B_n$
for odd n
are 0 except $B_1$
which is set to 0.5 on
purpose.
The Bernoulli polynomials can be directly defined as
$B_n(x) = \sum_{k=0}^n {n \choose k} b_{nk}\, x^k$
and it is immediately clear that the Bernoulli numbers are then given as
$B_n = B_n(0)$
.
Returns the first n+1
Bernoulli numbers, if x
is missing, or
the value of the Bernoulli polynomial at point(s) x
.
The definition uses B_1 = 1/2
in accordance with the definition of
the Bernoulli polynomials.
See the entry on Bernoulli numbers in the Wikipedia.
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 base polynomials and approximations.
bernstein(f, n, x)
bernsteinb(k, n, x)
f 
function to be approximated by Bernstein polynomials. 
k 
integer between 0 and n, the kth Bernstein polynomial of order n. 
n 
order of the Bernstein polynomial(s). 
x 
numeric scalar or vector where the Bernstein polynomials will be calculated. 
The Bernstein basis polynomials $B_{k,n}(x)$
are defined as
$B_{k,n}(x) = {{n}\choose{k}} x^k (1x)^{nk}$
and form a basis for the vector space of polynomials of degree
$n$
over the interval $[0,1]$
.
bernstein(f, n, x)
computes the approximation of function
f
through Bernstein polynomials of degree n
, resp.
computes the value of this approximation at x
. The function
is vectorized and applies a brute force calculation.
But if x
is a scalar, the value will be calculated using
De Casteljau's algorithm for higher accuracy. For bigger n
the binomial coefficients may be in for problems.
Returns a scalar or vector of function values.
See https://en.wikipedia.org/wiki/Bernstein_polynomial
## 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)
Finding roots of univariate functions in bounded intervals.
bisect(fun, a, b, maxiter = 500, tol = NA, ...)
secant(fun, a, b, maxiter = 500, tol = 1e08, ...)
regulaFalsi(fun, a, b, maxiter = 500, tol = 1e08, ...)
fun 
Function or its name as a string. 
a , b

interval end points. 
maxiter 
maximum number of iterations; default 100. 
tol 
absolute tolerance; default 
... 
additional arguments passed to the function. 
“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 socalled ‘Illinois’ improvement is used here.
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
Quarteroni, A., R. Sacco, and F. Saleri (2007). Numerical Mathematics. Second Edition, SpringerVerlag, Berlin Heidelberg.
bisect(sin, 3.0, 4.0)
# $root $f.root $iter $estim.prec
# 3.1415926536 1.2246467991e16 52 4.4408920985e16
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
Literal bit representation.
bits(x, k = 54, pos_sign = FALSE, break0 = FALSE)
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. 
The literal bit/binary representation of a floating point number is computed by subtracting powers of 2.
Returns a string containing the binary representation.
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"
Create a string of blank characters.
blanks(n)
n 
integer greater or equal to 0. 
blanks(n)
is a string of n
blanks.
String of n
blanks.
blanks(6)
Build a block diagonal matrix.
blkdiag(...)
... 
sequence of nonempty, numeric matrices 
Generate a block diagonal matrix from A, B, C, .... All the arguments must be numeric and nonempty matrices.
a numeric matrix
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.
a1 < matrix(c(1,2), 1)
a2 < as.matrix(c(1,2))
blkdiag(a1, diag(1, 2, 2), a2)
Find root of continuous function of one variable.
brentDekker(fun, a, b, maxiter = 500, tol = 1e12, ...)
brent(fun, a, b, maxiter = 500, tol = 1e12, ...)
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. 
brentDekker
implements a version of the BrentDekker algorithm,
a well known root finding algorithms for real, univariate, continuous
functions. The BrentDekker approach is a clever combination of secant
and bisection with quadratic interpolation.
brent
is simply an alias for brentDekker
.
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. 
Quarteroni, A., R. Sacco, and F. Saleri (2007). Numerical Mathematics. Second Edition, SpringerVerlag, Berlin Heidelberg.
# 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
The Brown72 data set represents a fractal Brownian motion with a prescribed Hurst exponent 0f 0.72 .
data(brown72)
The format is: one column.
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.
## Not run:
data(brown72)
plot(brown72, type = "l", col = "blue")
grid()
## End(Not run)
Broyden's method for the numerical solution of nonlinear systems of
n
equations in n
variables.
broyden(Ffun, x0, J0 = NULL, ...,
maxiter = 100, tol = .Machine$double.eps^(1/2))
Ffun 

x0 
Numeric vector of length 
J0 
Jacobian of the function at 
... 
additional parameters passed to the function. 
maxiter 
Maximum number of iterations. 
tol 
Tolerance, relative accuracy. 
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 rankone update thereafter, applying the socalled ShermanMorrison 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'.
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.
Applied to a system of n
linear equations it will stop in
2n
steps
Quarteroni, A., R. Sacco, and F. Saleri (2007). Numerical Mathematics. Second Edition, SpringerVerlag, Berlin Heidelberg.
## 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.092626e09; 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.34325e08; 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.284374e10
# 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)
Apply a binary function elementwise.
bsxfun(func, x, y)
arrayfun(func, ...)
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. 
bsxfun
applies elementbyelement 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.
The result will be a vector or matrix of the same size as x, y
.
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.
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)
BulirschStoer algorithm for solving Ordinary Differential Equations (ODEs) very accurately.
bulirsch_stoer(f, t, y0, ..., tol = 1e07)
midpoint(f, t0, tfinal, y0, tol = 1e07, kmax = 25)
f 
function describing the differential equation 
t 
vector of 
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. 
The BulirschStoer algorithm is a wellknown method to obtain highaccuracy 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.
BulirschStoer and midpoint shall not be used with nonsmooth 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.
bulirsch_stoer returns a list with x
the grid points input, and
y
a vector of function values at the se points.
Will be extended to become a fullblown BulirschStoer implementation.
Copyright (c) 2014 Hans W Borchers
J. Stoer and R. Bulirsch (2002). Introduction to Numerical Analysis. Third Edition, Texts in Applied Mathematics 12, Springer Science + Business, LCC, New York.
## 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 equallyspaced grid points
yy[nrow(yy), 1] # 1.1e11
## 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.5e11
S$y[nrow(S$y), 1] # 7.1e04
## 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)
Solves boundary value problems of linear second order differential equations.
bvp(f, g, h, x, y, n = 50)
f , g , h

functions on the right side of the differential equation.
If 
x 

y 
boundary conditions such that

n 
number of intermediate grid points; default 50. 
Solves the twopoint boundary value problem given as a linear differential equation of second order in the form:
$y'' = f(x) y' + g(x) y + h(x)$
with the finite element method. The solution $y(x)$
shall exist
on the interval $[a, b]$
with boundary conditions $y(a) = y_a$
and $y(b) = y_b$
.
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.
Uses a tridiagonal equation solver that may be faster then qr.solve
for large values of n
.
Kutz, J. N. (2005). Practical Scientific Computing. Lecture Notes 981952420, University of Washington, Seattle.
## 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)
Transforms between cartesian, spherical, polar, and cylindrical coordinate systems in two and three dimensions.
cart2sph(xyz)
sph2cart(tpr)
cart2pol(xyz)
pol2cart(prz)
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. 
The spherical coordinate system used here consists of
 theta
, azimuth angle relative to the positive xaxis
 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.
All functions return a (2 or 3dimensional) 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.
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.
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.123234e17 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))
Displays or changes working directory, or lists files therein.
cd(dname)
pwd()
what(dname = getwd())
dname 
(relative or absolute) directory path. 
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.
Name of the current working directory.
# cd()
# pwd()
# what()
Functions for rounding and truncating numeric values towards near integer values.
ceil(n)
Fix(n)
n 
a numeric vector or matrix 
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.
integer values
x < c(1.2, 0.8, 0, 0.5, 1.1, 2.9)
ceil(x)
Fix(x)
Computes the characteristic polynomial (and the inverse of the matrix, if requested) using the FaddeewLeverrier method.
charpoly(a, info = FALSE)
a 
quadratic matrix; size should not be much larger than 100. 
info 
logical; if true, the inverse matrix will also be reported. 
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).
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.
Hou, S.H. (1998). Classroom Note: A Simple Proof of the Leverrier–Faddeev Characteristic Polynomial Algorithm, SIAM Review, 40(3), pp. 706–709.
a < magic(5)
A < charpoly(a, info = TRUE)
A$cp
roots(A$cp)
A$det
zapsmall(A$inv %*% a)
Function approximation through Chebyshev polynomials (of the first kind).
chebApprox(x, fun, a, b, n)
x 
Numeric vector of points within interval 
fun 
Function to be approximated. 
a , b

Endpoints of the interval. 
n 
An integer 
Return approximate ycoordinates of points at x by computing the
Chebyshev approximation of degree n for fun
on the interval
[a, b]
.
A numeric vector of the same length as x
.
TODO: Evaluate the Chebyshev approximative polynomial by using the Clenshaw recurrence formula. (Not yet vectorized, that's why we still use the Horner scheme.)
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.
# 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(ysyp)) # 0.006925271
max(abs(ysyc)) # 1.151207e05
## 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 Coefficients for Chebyshev polynomials of the first kind.
chebCoeff(fun, a, b, n)
fun 
function to be approximated. 
a , b

endpoints of the interval. 
n 
an integer 
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).
Vector of coefficients for the Chebyshev polynomials, from low to high degrees (see the example).
See the “Chebfun Project” <https://www.chebfun.org/> by Nick Trefethen.
Weisstein, Eric W. “Chebyshev Polynomial of the First Kind." From MathWorld — A Wolfram Web Resource. https://mathworld.wolfram.com/ChebyshevPolynomialoftheFirstKind.html
## 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 and their values.
chebPoly(n, x = NULL)
n 
an integer 
x 
a numeric vector, possibly empty; default 
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.
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
.
See the “Chebfun Project” <https://www.chebfun.org/> by Nick Trefethen.
Carothers, N. L. (1998). A Short Course on Approximation Theory. Bowling Green State University.
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 from points in the plane
circlefit(xp, yp)
xp , yp

Vectors representing the x and y coordinates of plane points 
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.
Returns x and ycoordinates 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.
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.
# 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)
List or remove items from workspace, or display system information.
clear(lst)
ver()
who()
whos()
lst 
Character vector of names of variables in the global environment. 
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.
Invisibly NULL.
ls
, rm
, sessionInfo
# clear() # DON'T
# who()
# whos()
# ver()
ClenshawCurtis Quadrature Formula
clenshaw_curtis(f, a = 1, b = 1, n = 1024, ...)
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 
ClenshawCurtis quadrature is based on sampling the integrand on Chebyshev points, an operation that can be implemented using the Fast Fourier Transform.
Numerical scalar, the value of the integral.
Trefethen, L. N. (2008). Is Gauss Quadrature Better Than ClenshawCurtis? SIAM Review, Vol. 50, No. 1, pp 67–87.
## 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.3e10
Generates all combinations of length m
of a vector a
.
combs(a, m)
a 
numeric vector of some length 
m 
integer with 
combs
generates combinations of length n
of the elements
of the vector a
.
matrix representing combinations of the elements of a
combs(seq(2, 10, by=2), m = 3)
Computes the companion matrix of a real or complex vector.
compan(p)
p 
vector representing a polynomial 
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.
A square matrix of length(p)1
rows and columns
p < c(1, 0, 7, 6)
compan(p)
# 0 7 6
# 1 0 0
# 0 1 0
Complex step derivatives of realvalued functions, including gradients, Jacobians, and Hessians.
complexstep(f, x0, h = 1e20, ...)
grad_csd(f, x0, h = 1e20, ...)
jacobian_csd(f, x0, h = 1e20, ...)
hessian_csd(f, x0, h = 1e20, ...)
laplacian_csd(f, x0, h = 1e20, ...)
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 
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
.
complexstep(f, x0)
returns the derivative $f'(x_0)$
of $f$
at $x_0$
. The function is vectorized in x0
.
This surprising approach can be easily deduced from the complexanalytic Taylor formula.
HwB <[email protected]>
Martins, J. R. R. A., P. Sturdza, and J. J. Alonso (2003). The Complexstep Derivative Approximation. ACM Transactions on Mathematical Software, Vol. 29, No. 3, pp. 245–262.
## 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 1e10
# numDeriv::grad(f, 1.5) # 4.05342789388197, error 1e12 Richardson
# pracma::numderiv # 4.05342789389868, error 5e14 Richardson
complexstep(f, 1.5) # 4.05342789389862, error 1e15
# 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
Condition number of a matrix.
cond(M, p = 2)
M 
Numeric matrix; vectors will be considered as column vectors. 
p 
Indicates the 
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
wellconditioned matrix.
cond(M)
returns the 2norm 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.)
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.
Trefethen, L. N., and D. Bau III. (1997). Numerical Linear Algebra. SIAM, Philadelphia.
cond(hilb(8))
Convolution and polynomial multiplication.
conv(x, y)
x , y

real or complex vectors. 
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
.
Another vector.
conv
utilizes fast Fourier transformation.
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 not available in R.
cot(z)
csc(z)
sec(z)
acot(z)
acsc(z)
asec(z)
z 
numeric or complex scalar or vector. 
The usual trigonometric cotangens, cosecans, and secans functions and their inverses, computed through the other well known – in R – sine, cosine, and tangens functions.
Result vector of numeric or complex values.
These function names are available in Matlab, that is the reason they have been added to the ‘pracma’ package.
Trigonometric and hyperbolic functions in R.
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
Closed composite NewtonCotes formulas of degree 2 to 8.
cotes(f, a, b, n, nodes, ...)
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 NewtonCotes formula. 
... 
additional parameters to be passed to the function. 
2 to 8 point closed and summed NewtonCotes 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.
The integral as a scalar.
It is generally recommended not to apply NewtonCotes formula of degrees
higher than 6, instead increase the number n
of subintervals used.
Standard NewtonCotes 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.
Quarteroni, A., R. Sacco, and F. Saleri (2007). Numerical Mathematics. Second Edition, SpringerVerlag, Berlin Heidelberg.
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 not available in R.
coth(z)
csch(z)
sech(z)
acoth(z)
acsch(z)
asech(z)
z 
numeric or complex scalar or vector. 
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.
Result vector of numeric or complex values.
These function names are available in Matlab, that is the reason they have been added to the ‘pracma’ package.
Trigonometric and hyperbolic functions in R.
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
The CrankNicolson method for solving ordinary differential equations is a combination of the generic steps of the forward and backward Euler methods.
cranknic(f, t0, t1, y0, ..., N = 100)
f 
function in the differential equation 
t0 , t1

start and end points of the interval. 
y0 
starting values as row or column vector;
for 
N 
number of steps. 
... 
Additional parameters to be passed to the function. 
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.
List with components t
for grid (or ‘time’) points between t0
and t1
, and y
an nbym matrix with solution variables in
columns, i.e. each row contains one time stamp.
This is for demonstration purposes only; for real problems or applications
please use ode23
or rkf54
.
Quarteroni, A., and F. Saleri (2006). Scientific Computing With MATLAB and Octave. Second Edition, SpringerVerlag, Berlin Heidelberg.
## 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 or cross product
cross(x, y)
x 
numeric vector or matrix 
y 
numeric vector or matrix 
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 n1
vectors in ndimensional
space is realized as crossn
.
3dim. vector if x
and <
are vectors, a matrix of
3dim. vectors if x
and y
are matrices themselves.
cross(c(1, 2, 3), c(4, 5, 6)) # 3 6 3
Vector cross product of n1
vectors in ndimensional space
crossn(A)
A 
matrix of size 
The rows of the matrix A
are taken as(n1)
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.
a vector of length n
The ‘scalar triple product’ in $R^3$
can be defined as
spatproduct < function(a, b, c) dot(a, crossn(b, c))
It represents the volume of the parallelepiped spanned by the three vectors.
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)))
Computes the natural interpolation cubic spline.
cubicspline(x, y, xi = NULL, endp2nd = FALSE, der = c(0, 0))
x , y

x and ycoordinates of points to be interpolated. 
xi 
xcoordinates of points at which the interpolation is to be performed. 
endp2nd 
logical; if true, the derivatives at the endpoints are
prescribed by 
der 
a twocomponents vector prescribing derivatives at endpoints. 
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.
Returns either the interpolated values at the points xi
or, if
is.null(xi)
, the piecewise polynomial that represents the spline.
From the piecewise polynomial returned one can easily generate the spline function, see the examples.
Quarteroni, Q., and F. Saleri (2006). Scientific Computing with Matlab and Octave. SpringerVerlag Berlin Heidelberg.
## 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)
Polynomial fitting of parametrized points on 2D curves, also requiring to meet some points exactly.
curvefit(u, x, y, n, U = NULL, V = NULL)
u 
the parameter vector. 
x , y

x, ycoordinates for each parameter value. 
n 
order of the polynomials, the same in x and ydirction. 
U 
parameter values where points will be fixed. 
V 
matrix with two columns and 
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.
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 ydirection.
In the same manner, derivatives/directions could be prescribed at certain points.
## 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)
Finds cutting points for vector s of real numbers.
cutpoints(x, nmax = 8, quant = 0.95)
x 
vector of real values. 
nmax 
the maximum number of cutting points to choose 
quant 
quantile of the gaps to consider for cuts. 
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.
Returns a list with components cutp
, the cutting points selected,
and cutd
, the gap between values of x
at this cutting point.
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 RForge project “Rulebased Models”) is quite cumbersome to use.
Witten, I. H., and E. Frank (2005). Data Mining: Practical Machine Learning Tools and Techniques. Morgan Kaufmann Publishers, San Francisco.
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 kmeans clustering
km < kmeans(x, n)
points(x, rep(0, N), col = km$cluster, pch = "+")
## A 2dimensional 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)
Numerically evaluate double integral over rectangle.
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, ...)
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 
... 
additional parameters to be passed to the integrand. 
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.
Numerical scalar, the value of the integral.
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 = 2e8
f4 < function(x, y, z) y*sin(x)+z*cos(x)
triplequad(f4, 0,pi, 0,1, 1,1) #  2.0 => 2.220446e16
Deconvolution and polynomial division.
deconv(b, a)
b , a

real or complex vectors. 
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.
List with elements named q
and r
.
TODO: Base deconv
on some filter1d
function.
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
Detect events in solutions of a differential equation.
deeve(x, y, yv = 0, idx = NULL)
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. 
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.
A (time) point x0
at which the event happens.
The interpolation is linear only for the moment.
## 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()
Transforms between angles in degrees and radians.
deg2rad(deg)
rad2deg(rad)
deg 
(array of) angles in degrees. 
rad 
(array of) angles in radians. 
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.
The angle in degrees or radians.
deg2rad(c(0, 10, 20, 30, 40, 50, 60, 70, 80, 90))
rad2deg(seq(pi/2, pi/2, length = 19))
Removes the mean value or (piecewise) linear trend from a vector or from each column of a matrix.
detrend(x, tt = 'linear', bp = c())
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 
detrend
computes the leastsquares 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 straightline fit, use polyfit
.
removes the mean or (piecewise) linear trend from x
and returns it
in y=detrend(x)
, that is xy
is the linear trend.
Detrending is often used for FFT processing.
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 solution of a differential equation solver.
deval(x, y, xp, idx = NULL)
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. 
Determines where the points xp
lie within the vector x
and interpolates linearly.
An length(xp)
bylength(idx)
matrix of values at points
xp
.
The interpolation is linear only for the moment.
## 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()
Generate diagonal matrices or return diagonal of a matrix
Diag(x, k = 0)
x 
vector or matrix 
k 
integer indicating a secondary diagonal 
If x
is a vector, Diag(x, k)
generates a matrix with x
as the (kth 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
.
matrix or vector
In Matlab/Octave this function is called diag()
and has a different
signature than the corresponding function in R.
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)
Display text or array, or produce beep sound.
disp(...)
beep()
... 
any R object that can be printed. 
Display text or array, or produces the computer's default beep sound using ‘cat’ with closing newline.
beep() returns NULL invisibly, disp() displays with newline.
disp("Some text, and numbers:", pi, exp(1))
# beep()
Computes the Euclidean distance between rows of two matrices.
distmat(X, Y)
pdist(X)
pdist2(X, Y)
X 
matrix of some size 
Y 
matrix of some size 
Computes Euclidean distance between two vectors A and B as:
AB = 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)
.
matrix of size m x n
if x
is of size m x k
and
y
is of size n x k
.
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 10100 times faster, utilizing the similarity between Euclidean distance and matrix operations.
Copyright (c) 1999 Roland Bunschoten for a Matlab version on MatlabCentral
under the name distance.m
. Translated to R by Hans W Borchers.
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))
}
'dot' or 'scalar' product of vectors or pairwise columns of matrices.
dot(x, y)
x 
numeric vector or matrix 
y 
numeric vector or matrix 
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.
A scalar or vector of length the number of columns of x
and
y
.
dot(1:5, 1:5) #=> 55
# Length of space diagonal in 3dim cube:
sqrt(dot(c(1,1,1), c(1,1,1))) #=> 1.732051
Eigenvalues of a matrix
eig(a)
a 
real or complex square matrix 
Computes the eigenvalues of a square matrix of real or complex numbers,
using the R routine eigen
without computing the eigenvectors.
Vector of eigenvalues
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's iteration method for eigenvalues and eigenvectors.
eigjacobi(A, tol = .Machine$double.eps^(2/3))
A 
a real symmetric matrix. 
tol 
requested tolerance. 
The Jacobi eigenvalue method repeatedly performs (Givens) transformations until the matrix becomes almost diagonal.
Returns a list with components V
, a matrix containing the
eigenvectors as columns, and D
a vector of the eigenvalues.
This R implementation works well up to 50x50matrices.
Mathews, J. H., and K. D. Fink (2004). Numerical Methods Using Matlab. Fourth edition, Pearson education, Inc., New Jersey.
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.
einsteinF(d, x)
x 
numeric or complex vector. 
d 
parameter to select one of the Einstein functions E1, E2, E3, E4. 
The Einstein functions are sometimes used for the PlanckEinstein oscillator in one degree of freedom.
The functions are defined as:
$E1(x) = \frac{x^2 e^x}{(e^x  1)^2}$
$E2(x) = \frac{x}{e^x  1}$
$E3(x) = ln(1  e^{x})$
$E4(x) = \frac{x}{e^x  1}  ln(1  e^{x})$
E1
has an inflection point as x=2.34694130...
.
Numeric/complex scalar or vector.
## 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)
Complete elliptic integrals of the first and second kind, and Jacobi elliptic integrals.
ellipke(m, tol = .Machine$double.eps)
ellipj(u, m, tol = .Machine$double.eps)
u 
numeric vector. 
m 
input vector, all input elements must satisfy 
tol 
tolerance; default is machine precision. 
ellipke
computes the complete elliptic integrals to accuracy
tol
, based on the algebraicgeometric mean.
ellipj
computes the Jacobi elliptic integrals sn
, cn
,
and dn
. For instance, $sn$
is the inverse function for
$u = \int_0^\phi dt / \sqrt{1  m \sin^2 t}$
with $sn(u) = \sin(\phi)$
.
Some definitions of the elliptic functions use the modules k
instead
of the parameter m
. They are related by k^2=m=sin(a)^2
where
a
is the ‘modular angle’.
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
.
Abramowitz, M., and I. A. Stegun (1965). Handbook of Mathematical Functions. Dover Publications, New York.
elliptic::sn,cn,dn
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
Distance from 1.0 to the next largest doubleprecision number.
eps(x = 1.0)
x 
scalar or numerical vector or matrix. 
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)))
.
Returns a scalar.
for (i in 5:5) cat(eps(10^i), "\n")
# 1.694066e21
# 1.355253e20
# 2.168404e19
# 1.734723e18
# 1.387779e17
# 2.220446e16
# 1.776357e15
# 1.421085e14
# 1.136868e13
# 1.818989e12
# 1.455192e11
The error or Phi function is a variant of the cumulative normal (or Gaussian) distribution.
erf(x)
erfinv(y)
erfc(x)
erfcinv(y)
erfcx(x)
erfz(z)
erfi(z)
x , y

vector of real numbers. 
z 
real or complex number; must be a scalar. 
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.
Real or complex number(s), the value(s) of the function.
For the complex error function we used Fortran code from the book S. Zhang & J. Jin “Computation of Special Functions” (Wiley, 1996).
First version by Hans W Borchers;
vectorized version of erfz
by Michael Lachmann.
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
Draws symmetric error bars in x and/or ydirection.
errorbar(x, y, xerr = NULL, yerr = NULL,
bar.col = "red", bar.len = 0.01,
grid = TRUE, with = TRUE, add = FALSE, ...)
x , y

x, ycoordinates 
xerr , yerr

length of the error bars, relative to the x, yvalues. 
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

errorbar
plots y
versus x
with symmetric error bars,
with a length determined by xerr
resp. yerr
in x and/or
ydirection. 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))
.
Generates a plot, no return value.
plotrix::plotCI
, Hmisc::errbar
## 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's eta function valid in the entire complex plane.
eta(z)
z 
Real or complex number or a numeric or complex vector. 
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.
Returns a complex vector of function values.
Copyright (c) 2001 Paul Godfrey for a Matlab version available on Mathwork's Matlab Central under BSD license.
Zhang, Sh., and J. Jin (1996). Computation of Special Functions. WileyInterscience, New York.
z < 0.5 + (1:5)*1i
eta(z)
z < c(0, 0.5+1i, 1, 1i, 2+2i, 1, 2, 11i)
eta(z)
Euler and EulerHeun ODE solver.
euler_heun(f, a, b, y0, n = 100, improved = TRUE, ...)
f 
function in the differential equation 
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. 
euler_heun
is an integration method for ordinary differential
equations using the simple Euler resp. the improved EulerHeun Method.
List with components t
for grid (or ‘time’) points, and y
the vector of predicted values at those grid points.
Quarteroni, A., and F. Saleri (). Scientific Computing with MATLAB and Octave. Second Edition, SpringerVerlag, Berlin Heidelberg, 2006.
## Flameup 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)
The exponential integral functions E1 and Ei and the logarithmic integral Li.
The exponential integral is defined for $x > 0$
as
$\int_x^\infty \frac{e^{t}}{t} dt$
and by analytic continuation in the complex plane. It can also be defined as the Cauchy principal value of the integral
$\int_{\infty}^x \frac{e^t}{t} dt$
This is denoted as $Ei(x)$
and the relationship between Ei
and
expint(x)
for x real, x > 0 is as follows:
$Ei(x) =  E1(x) i \pi$
The logarithmic integral $li(x)$
for real $x, x > 0$
, is defined as
$li(x) = \int_0^x \frac{dt}{log(t)}$
and the Eulerian logarithmic integral as $Li(x) = li(x)  li(2)$
.
The integral $Li$
approximates the prime number function $\pi(n)$
,
i.e., the number of primes below or equal to n (see the examples).
expint(x)
expint_E1(x)
expint_Ei(x)
li(x)
x 
vector of real or complex numbers. 
For x
in [38, 2]
we use a series expansion,
otherwise a continued fraction, see the references below, chapter 5.
Returns a vector of real or complex numbers, the vectorized exponential integral, resp. the logarithmic integral.
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”.
Abramowitz, M., and I.A. Stegun (1965). Handbook of Mathematical Functions. Dover Publications, New York.
gsl::expint_E1,expint_Ei
, primes
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)
Computes the exponential of a matrix.
expm(A, np = 128)
logm(A)
A 
numeric square matrix. 
np 
number of points to use on the unit circle. 
For an analytic function $f$
and a matrix $A$
the expression
$f(A)$
can be computed by the Cauchy integral
$f(A) = (2 \pi i)^{1} \int_G (zIA)^{1} f(z) dz$
where $G$
is a closed contour around the eigenvalues of $A$
.
Here this is achieved by taking G to be a circle and approximating the integral by the trapezoid rule.
logm
is a fake at the moment as it computes the matrix logarithm
through taking the logarithm of its eigenvalues; will be replaced by an
approach using Pade interpolation.
Another more accurate and more reliable approach for computing these functions can be found in the R package ‘expm’.
Matrix of the same size as A
.
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).
Idea and Matlab code for a cubic root by Nick Trefethen in his “10 digits 1 page” project.
Moler, C., and Ch. Van Loan (2003). Nineteen Dubious Ways to Compute the Exponential of a Matrix, TwentyFive 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.
expm::expm
## 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
Create basic matrices.
eye(n, m = n)
ones(n, m = n)
zeros(n, m = n)
m , n

numeric scalars specifying size of the matrix 
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.
Diag
,
eye(3)
ones(3, 1)
zeros(1, 3)
Easytouse contour and 3D surface resp mesh plotter.
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, ...)
f 
2D function to be plotted, must accept 
xlim , ylim

defines x and yranges 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 
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.
Plots the function graph and invisibly returns NULL
.
Mimicks Matlab functions of the same names; Matlab's ezcontourf
can
be generated with filled=TRUE
.
## Not run:
f < function(xy) {
x < xy[1]; y < xy[2]
3*(1x)^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 w/o the need to define x, y
coordinates.
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", ...)
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 
gridcol 
Color of grid points. 
fill 
Logical; shall the area between function and axis be filled?;
default: 
fillcol 
Color of fill area. 
xlab 
Label on the 
ylab 
Label on the 
main 
Title of the plot 
... 
More parameters to be passed to 
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.
Plots the function graph and invisibly returns NULL
.
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.
## 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 function plot w/o the need to define x, y
coordinates.
ezpolar(fun, interv = c(0, 2*pi))
fun 
function to be plotted. 
interv 
left and right endpoint for the plot. 
Calculates the x, y
coordinates of points to be plotted and
calls the polar
function.
Plots the function graph and invisibly returns NULL
.
Mimick the Matlab function of the same name.
## Not run:
fun < function(x) 1 + cos(x)
ezpolar(fun)
## End(Not run)
Factorial for nonnegative integers n <= 170
.
fact(n)
factorial2(n)
n 
Vector of integers, for 
The factorial is computed by brute force; factorials for n >= 171
are not representable as ‘double’ anymore.
fact
returns the factorial of each element in n
.
If n < 0
the value is NaN
, and for n > 170
it is Inf
.
Nonintegers 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.
The R core function factorial
uses the gamma
function,
whose implementation is not accurate enough for larger input values.
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)
Returns a vector containing the prime factors of n
.
factors(n)
n 
nonnegative integer 
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.
Vector containing the prime factors of n
.
## 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 function differentiation for orders n=1..4
using
finite difference approximations.
fderiv(f, x, n = 1, h = 0,
method = c("central", "forward", "backward"), ...)
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 
h 
step size: if 
method 
one of “central”, “forward”, or “backward”. 
... 
more variables to be passed to function 
Derivatives are computed applying central difference formulas that stem
from the Taylor series approximation. These formulas have a convergence
rate of $O(h^2)$
.
Use the ‘forward’ (right side) or ‘backward’ (left side) method if the function can only be computed or is only defined on one side. Otherwise, always use the central difference formulas.
Optimal step sizes depend on the accuracy the function can be computed with.
Assuming internal functions with an accuracy 2.2e16, appropriate step
sizes might be 5e6, 1e4, 5e4, 2.5e3
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”.
Vector of the same length as x
.
Numerical differentiation suffers from the conflict between roundoff and truncation errors.
Kiusalaas, J. (2005). Numerical Methods in Engineering with Matlab. Cambridge University Press.
## 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 for function minimum.
fibsearch(f, a, b, ..., endp = FALSE, tol = .Machine$double.eps^(1/2))
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 
... 
Additional arguments to be passed to f. 
Fibonacci search for a univariate function minimum in a bounded interval.
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
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
Open, activate, and close grahics devices.
figure(figno, title = "")
figno 
(single) number of plot device. 
title 
title of the plot device; not yet used. 
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.
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.
Does not bring the activated graphics device in front.
dev.set, dev.off, dev.list
## Not run:
figure()
figure(2)
## End(Not run)
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]
.
findintervals(x, xs)
x 
single number. 
xs 
numeric vector, not necessarily sorted. 
Contrary to findInterval
, the vector xs
in
findintervals
need not be sorted.
Vector of indices in 1..length(xs)
.
If none is found, returns integer(0)
.
If x
is equal to the last element in xs
, the index
length(xs)
will also be returned.
xs < zapsmall(sin(seq(0, 10*pi, len=100)))
findintervals(0, xs)
# 1 10 20 30 40 50 60 70 80 90 100
Finding all local(!) minima of a unvariate function in an interval by splitting the interval in many small subintervals.
findmins(f, a, b, n = 100, tol = .Machine$double.eps^(2/3), ...)
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. 
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.
Numeric vector with the xpositions of all minima found in the interval.
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 (maxima) in a time series.
findpeaks(x, nups = 1, ndowns = nups, zero = "0", peakpat = NULL,
minpeakheight = Inf, minpeakdistance = 1,
threshold = 0, npeaks = 0, sortstr = FALSE)
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 
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 
This function is quite general as it relies on regular patterns to determine where a peak is located, from beginning to end.
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.
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.
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)
Finds indices of nonzero elements.
finds(v)
v 
logical or numeric vector or array 
Finds indices of true or nonzero elements of argument v
;
can be used with a logical expression.
Indices of elements matching the expression x
.
finds(3:3 >= 0)
finds(c(0, 1, 0, 2, 3))
Finding all roots of a unvariate function in an interval by splitting the interval in many small subintervals.
findzeros(f, a, b, n = 100, tol = .Machine$double.eps^(2/3), ...)
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. 
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.
Numeric vector with the xpositions of all roots found in the interval.
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)
Conjugate Gradient (CG) minimization through the DavidonFletcherPowell approach for function minimization.
The DavidonFletcherPowell (DFP) and the BroydenFletcherGoldfarbShanno (BFGS) methods are the first quasiNewton minimization methods developed. These methods differ only in some details; in general, the BFGS approach is more robust.
fletcher_powell(x0, f, g = NULL,
maxiter = 1000, tol = .Machine$double.eps^(2/3))
x0 
start value. 
f 
function to be minimized. 
g 
gradient function of 
maxiter 
max. number of iterations. 
tol 
relative tolerance, to be used as stopping rule. 
The starting point is Newton's method in the multivariate case, when the estimate of the minimum is updated by the following equation
$x_{new} = x  H^{1}(x) grad(g)(x)$
where $H$
is the Hessian and $grad$
the gradient.
The basic idea is to generate a sequence of good approximations to the inverse Hessian matrix, in such a way that the approximations are again positive definite.
List with following components:
xmin 
minimum solution found. 
fmin 
value of 
niter 
number of iterations performed. 
Used some Matlab code as described in the book “Applied Numerical Analysis Using Matlab” by L. V.Fausett.
J. F. Bonnans, J. C. Gilbert, C. Lemarechal, and C. A. Sagastizabal. Numerical Optimization: Theoretical and Practical Aspects. Second Edition, SpringerVerlag, Berlin Heidelberg, 2006.
## Rosenbrock function
rosenbrock < function(x) {
n < length(x)
x1 < x[2:n]
x2 < x[1:(n1)]
sum(100*(x1x2^2)^2 + (1x2)^2)
}
fletcher_powell(c(0, 0), rosenbrock)
# $xmin
# [1] 1 1
# $fmin
# [1] 1.774148e27
# $niter
# [1] 14
Flip matrices up and down or left and right; or circulating indices per dimension.
flipdim(a, dim)
flipud(a)
fliplr(a)
circshift(a, sz)
a 
numeric or complex matrix 
dim 
flipping dimension; can only be 1 (default) or 2 
sz 
integer vector of length 1 or 2. 
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).
the original matrix somehow flipped or circularly shifted.
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)
Find minimum of singlevariable function on fixed interval.
fminbnd(f, a, b, maxiter = 1000, maximum = FALSE,
tol = 1e07, rel.tol = tol, abs.tol = 1e15, ...)
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. 
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.
List with
xmin 
location of the minimum resp. maximum. 
fmin 
function value at the optimum. 
niter 
number of iterations used. 
estim.prec 
estimated precision. 
fminbnd
mimics the Matlab function of the same name.
R. P. Brent (1973). Algorithms for Minimization Without Derivatives. Dover Publications, reprinted 2002.
## 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) # complexstep derivative
xs < findzeros(g, 1, 1) # local minima and maxima
ys < f(xs); n0 < which.min(ys) # index of global minimum
fminbnd(f, xs[n01], 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)
Find minimum of multivariable functions with nonlinear constraints.
fmincon(x0, fn, gr = NULL, ..., method = "SQP",
A = NULL, b = NULL, Aeq = NULL, beq = NULL,
lb = NULL, ub = NULL, hin = NULL, heq = NULL,
tol = 1e06, maxfeval = 10000, maxiter = 5000)
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. 
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.
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. 
fmincon
mimics the Matlab function of the same name.
Xianyan Chen for the package NlcOptim.
J. Nocedal and S. J. Wright (2006). Numerical Optimization. Second Edition, Springer Science+Business Media, New York.
# 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
Find minimum of multivariable functions using derivativefree methods.
fminsearch(fn, x0, ..., lower = NULL, upper = NULL,
method = c("NelderMead", "HookeJeeves"),
minimize = TRUE, maxiter = 1000, tol = 1e08)
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 
"NelderMead" (default) or "HookeJeeves"; can be abbreviated. 
minimize 
logical; shall a minimum or a maximum be found. 
maxiter 
maximal number of iterations 
tol 
relative tolerance. 
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 "NelderMead" and "HookeJeeves" are available. Only
HookeJeeves can handle bounds constraints. For nonlinear constraints see
fmincon
, and for methods using gradients see fminunc
.
Important: fminsearch
may only give local solutions.
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. 
fminsearch
mimics the Matlab function of the same name.
Nocedal, J., and S. Wright (2006). Numerical Optimization. Second Edition, SpringerVerlag, New York.
# Rosenbrock function
rosena < function(x, a) 100*(x[2]x[1]^2)^2 + (ax[1])^2 # min: (a, a^2)
fminsearch(rosena, c(1.2, 1), a = sqrt(2), method="NelderMead")
## $xmin $fmin
## [1] 1.414292 2.000231 [1] 1.478036e08
fminsearch(rosena, c(1.2, 1), a = sqrt(2), method="HookeJeeves")
## $xmin $fmin
## [1] 1.414215 2.000004 [1] 1.79078e12
Find minimum of unconstrained multivariable functions.
fminunc(x0, fn, gr = NULL, ...,
tol = 1e08, maxiter = 0, maxfeval = 0)
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. 
The method used here for unconstrained minimization is a variant of a "variable metric" resp. quasiNewton approach.
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. 
fminunc
mimics the Matlab function of the same name.
The "variable metric" code provided by John Nash (package Rvmmin), strippeddown version by Hans W. Borchers.
J. Nocedal and S. J. Wright (2006). Numerical Optimization. Second Edition, Springer Science+Business Media, New York.
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
The fnorm
function calculates several different types of function
norms for depending on the argument p
.
fnorm(f, g, x1, x2, p = 2, npoints = 100)
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. 
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.
Numeric scalar (or Inf
), or NA
if one of these functions
returns NA
.
Another kind of ‘mean’ distance could be calculated by integrating the
difference fg
and dividing through the length of the interval.
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
Finite difference approximation using Fornberg's method for the derivatives of order 1 to k based on irregulat grid values.
fornberg(x, y, xs, k = 1)
x 
grid points on the xaxis, 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, 
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
.
Returns a matrix of size (length(xs))
, where the (k+1)th column
gives the value of the kth derivative. Especially the first column returns
the polynomial interpolation of the function.
Fornberg's method is considered to be numerically more stable than applying Vandermonde's matrix.
LeVeque, R. J. (2007). Finite Difference Methods for Ordinary and Partial Differential Equations. Society for Industrial and Applied Mathematics (SIAM), Philadelphia.
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 to stdout or a file.
fprintf(fmt, ..., file = "", append = FALSE)
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 
fprintf
applies the format string fmt
to all input
data ...
and writes the result to standard output or a file.
The usual Cstyle string formatting commands are used
Returns invisibly the number of bytes printed (using nchar
).
## Examples:
nbytes < fprintf("Results are:\n", file = "")
for (i in 1:10) {
fprintf("%4d %15.7f\n", i, exp(i), file = "")
}
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.
fractalcurve(n, which = c("hilbert", "sierpinski", "snowflake",
"dragon", "triangle", "arrowhead", "flowsnake", "molecule"))
n 
integer, the ‘order’ of the curve 
which 
character string, which curve to cumpute. 
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.
Returns a list with x, y
the x resp. ycoordinates of the
generated points describing the fractal curve.
Copyright (c) 2011 Jonas Lundgren for the Matlab toolbox fractal
curves
available on MatlabCentral under BSD license;
here reimplemented in R with explicit allowance from the author.
Peitgen, H.O., H. Juergens, and D. Saupe (1993). Fractals for the Classroom. SpringerVerlag Berlin Heidelberg.
## The Hilbert curve transforms a 2dim. 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)
(Normalized) Fresnel integrals S(x) and C(x)
fresnelS(x)
fresnelC(x)
x 
numeric vector. 
The normalized Fresnel integrals are defined as
$S(x) = \int_0^x \sin(\pi/2 \, t^2) dt$
$C(x) = \int_0^x \cos(\pi/2 \, t^2) dt$
This program computes the Fresnel integrals S(x) and C(x) using Fortran code by Zhang and Jin. The accuracy is almost up to Machine precision.
The functions are not (yet) truly vectorized, but use a call to ‘apply’.
The underlying function .fresnel
(not exported) computes single
values of S(x)
and C(x)
at the same time.
Numeric vector of function values.
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.
Zhang, S., and J. Jin (1996). Computation of Special Functions. WileyInterscience.
## Compute Fresnel integrals through GaussLegendre 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 a system of m
nonlinear equations of n
variables.
fsolve(f, x0, J = NULL,
maxiter = 100, tol = .Machine$double.eps^(0.5), ...)
f 
function describing the system of equations. 
x0 
point near to the root. 
J 
Jacobian function of 
maxiter 
maximum number of iterations in 
tol 
tolerance to be used in GaussNewton. 
... 
additional variables to be passed to the function. 
fsolve
tries to solve the components of function f
simultaneously and uses the GaussNewton method with numerical gradient
and Jacobian. If m = n
, it uses broyden
. Not applicable
for univariate root finding.
List with
x 
location of the solution. 
fval 
function value at the solution. 
fsolve
mimics the Matlab function of the same name.
Antoniou, A., and W.S. Lu (2007). Practical Optimization: Algorithms and Engineering Applications. Springer Science+Business Media, New York.
## 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)
Find root of continuous function of one variable.
fzero(fun, x, maxiter = 500, tol = 1e12, ...)
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. 
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 nontrivially: 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.
fzero
returns a list with
x 
location of the root. 
fval 
function value at the root. 
fzero
mimics the Matlab function of the same name, but is translated
from Octave's fzero
function, copyrighted (c) 2009 by Jaroslav Hajek.
Alefeld, Potra and Shi (1995). Enclosing Zeros of Continuous Functions. ACM Transactions on Mathematical Software, Vol. 21, No. 3.
fzero(sin, 3) # 3.141593
fzero(cos,c(1, 2)) # 1.570796
fzero(function(x) x^32*x5, 2) # 2.094551
Find the root of a complex function
fzsolve(fz, z0)
fz 
complex(analytic) function. 
z0 
complex point near the assumed root. 
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.
Complex point with sufficiently small function value.
fz < function(z) sin(z)^2 + sqrt(z)  log(z)
fzsolve(fz, 1+1i)
# 0.2555197+0.8948303i
Lower and upper incomplete gamma function.
gammainc(x, a)
incgam(x, a)
x 
positive real number. 
a 
real number. 
gammainc
computes the lower and upper incomplete gamma
function, including the regularized gamma function. The lower and
upper incomplete gamma functions are defined as
$\gamma(x, a) = \int_0^x e^{t} \, t^{a1} \, dt$
and
$\Gamma(x, a) = \int_x^{\infty} e^{t} \, t^{a1} \, dt$
while the regularized incomplete gamma function is
$\gamma(x, a)/\Gamma(a)$
.
incgam
(a name used in Pari/GP) computes the upper incomplete
gamma function alone, applying the R function pgamma
. The
accuracy is thus much higher. It works for a >= 1
, for even
smaller values a recursion will give the result.
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.
Directly converting Fortran code is often easier than translating Matlab code generated with f2matlab.
Zhang, Sh., and J. Jin (1996). Computation of Special Functions. WileyInterscience, New York.
gammainc( 1.5, 2)
gammainc(1.5, 2)
incgam(3, 1.2)
incgam(3, 0.5); incgam(3, 0.5)
Gamma function valid in the entire complex plane.
gammaz(z)
z 
Real or complex number or a numeric or complex vector. 
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))
.
Returns a complex vector of function values.
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.
Zhang, Sh., and J. Jin (1996). Computation of Special Functions. WileyInterscience, New York.
gamma
, gsl::lngamma_complex
max(gamma(1:10)  gammaz(1:10))
gammaz(1)
gammaz(c(22i, 11i, 0, 1+1i, 2+2i))
# Euler's reflection formula
z < 1+1i
gammaz(1z) * gammaz(z) # == pi/sin(pi*z)
Simple GaussianKronrod quadrature formula.
gauss_kronrod(f, a, b, ...)
f 
function to be integrated. 
a , b

end points of the interval. 
... 
variables to be passed to the function. 
Gaussian quadrature of degree 7 with GaussKronrod of degree 15 for error
estimation, the quadQK15
procedure in the QUADPACK library.
List of value and relative error.
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.
Fausett, L. V. (2007). Applied Numerical Analysis Using Matlab. Second edition, Prentice Hall.
gauss_kronrod(sin, 0, pi) # 2.000000000000000 , rel.error: 1.14e12
gauss_kronrod(exp, 0, 1) # 1.718281828459045 , rel.error: 0
# 1.718281828459045 , i.e. exp(1)  1
Nodes and weights for the npoint GaussHermite quadrature formula.
gaussHermite(n)
n 
Number of nodes in the interval 
GaussHermite quadrature is used for integrating functions of the form
$\int_{\infty}^{\infty} f(x) e^{x^2} dx$
over the infinite interval $]\infty, \infty[$
.
x
and w
are obtained from a tridiagonal eigenvalue problem.
The value of such an integral is then sum(w*f(x))
.
List with components x
, the nodes or points in]Inf, Inf[
, and
w
, the weights applied at these nodes.
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.
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.
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
Nodes and weights for the npoint GaussLaguerre quadrature formula.
gaussLaguerre(n, a = 0)
n 
Number of nodes in the interval 
a 
exponent of 
GaussLaguerre quadrature is used for integrating functions of the form
$\int_0^{\infty} f(x) x^a e^{x} dx$
over the infinite interval $]0, \infty[$
.
x
and w
are obtained from a tridiagonal eigenvalue problem.
The value of such an integral is then sum(w*f(x))
.
List with components x
, the nodes or points in[0, Inf[
, and
w
, the weights applied at these nodes.
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.
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.
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
Nodes and weights for the npoint GaussLegendre quadrature formula.
gaussLegendre(n, a, b)
n 
Number of nodes in the interval 
a , b

lower and upper limit of the integral; must be finite. 
x
and w
are obtained from a tridiagonal eigenvalue problem.
List with components x
, the nodes or points in[a,b]
, and
w
, the weights applied at these nodes.
Gauss quadrature is not suitable for functions with singularities.
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.
## Quadrature with GaussLegendre 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: < 1e15
# 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 = 1e15)
# Use GaussKronrod approach for error estimation
cc < gaussLegendre(103, 1, 1)
abs(Q  sum(cc$w * f(cc$x))) # rel.error < 1e10
# Use GaussHermite for vectorvalued 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 < 1e15
GaussNewton method of minimizing a term $f_1(x)^2 + \ldots + f_m(x)^2$
or $F' F$
where $F = (f_1, \ldots, f_m)$
is a multivariate function
of $n$
variables, not necessarily $n = m$
.
gaussNewton(x0, Ffun, Jfun = NULL,
maxiter =100, tol = .Machine$double.eps^(1/2), ...)
Ffun 

Jfun 
function returning the Jacobian matrix of 
x0 
Numeric vector of length 
maxiter 
Maximum number of iterations. 
tol 
Tolerance, relative accuracy. 
... 
Additional parameters to be passed to f. 
Solves the system of equations applying the GaussNewton's method. It is especially designed for minimizing a sumofsquares 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.
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, andrelerr
the absoulte distance between the last two solutions.
If n=m
then directly applying the newtonsys
function might
be a better alternative.
Antoniou, A., and W.S. Lu (2007). Practical Optimization: Algorithms and Engineering Applications. Springer Business+Science, New York.
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
Greatest common divisor and least common multiple
gcd(a, b, extended = FALSE)
Lcm(a, b)
a , b

vectors of integers. 
extended 
logical; if 
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.
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
.
The following relation is always true:
n * m = gcd(n, m) * lcm(n, m)
numbers::extGCD
gcd(12, 1:24)
gcd(46368, 75025) # Fibonacci numbers are relatively prime to each other
Lcm(12, 1:24)
Lcm(46368, 75025) # = 46368 * 75025
Compute the “geometric median” of points in ndimensional space, that is the point with the least sum of (Euclidean) distances to all these points.
geo_median(P, tol = 1e07, maxiter = 200)
P 
matrix of points, 
tol 
relative tolerance. 
maxiter 
maximum number of iterations. 
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.
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.
This is also known as the “1median problem” and can be generalized to the
“kmedian problem” for k cluster centers;
see kcca
in the ‘flexclust’ package.
See Wikipedia's entry on “Geometric median”.
# Generate 100 points on the unit sphere in the 10dim. 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.069063e08
# $niter
# [1] 10
Geometric and harmonic mean along a dimension of a vector, matrix, or
array.trimmean
is almost the same as mean
in R.
geomean(x, dim = 1)
harmmean(x, dim = 1)
trimmean(x, percent = 0)
x 
numeric vector, matrix, or array. 
dim 
dimension along which to take the mean; 
percent 
percentage, between 0 and 100, of trimmed values. 
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).
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.
To have an exact analogue of mean(x)
in Matlab,
apply trimmean(x)
.
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 Rotations and QR decomposition
givens(A)
A 
numeric square matrix. 
givens(A)
returns a QR decomposition (or factorization) of the
square matrix A
by applying unitary 2by2 matrices U
such
that U * [xk;xl] = [x,0]
where x=sqrt(xk^2+xl^2)
List with two matrices Q
and R
, Q
orthonormal and
R
upper triangular, such that A=Q%*%R
.
Golub, G. H., and Ch. F. van Loan (1996). Matrix Computations. Third edition, John Hopkins University Press, Baltimore.
## 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))
gmres(A,b)
attempts to solve the system of linear equations
A*x=b
for x
.
gmres(A, b, x0 = rep(0, length(b)),
errtol = 1e6, kmax = length(b)+1, reorth = 1)
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. 
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!)
Returns a list with components x
the solution, error
the
vector of residual norms, and niter
the number of iterations.
Based on Matlab code from C. T. Kelley's book, see references.
C. T. Kelley (1995). Iterative Methods for Linear and Nonlinear Equations. SIAM, Society for Industrial and Applied Mathematics, Philadelphia, USA.
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.49173e01 1.22147e01 1.39901e02 1.37817e02 2.81713e31
#
# $niter
# [1] 5
Golden Ratio search for a univariate function minimum in a bounded interval.
golden_ratio(f, a, b, ..., maxiter = 100, tol = .Machine$double.eps^0.5)
f 
Function or its name as a string. 
a , b

endpoints of the interval. 
maxiter 
maximum number of iterations. 
tol 
absolute tolerance; default 
... 
Additional arguments to be passed to f. 
‘Golden ratio’ search for a univariate function minimum in a bounded interval.
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
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 function gradient.
grad(f, x0, heps = .Machine$double.eps^(1/3), ...)
f 
function of several variables. 
x0 
point where the gradient is to build. 
heps 
step size. 
... 
more variables to be passed to function 
Computes the gradient
$(\frac{\partial f}{\partial x_1}, \ldots, \frac{\partial f}{\partial x_n})$
numerically using the “central difference formula”.
Vector of the same length as x0
.
Mathews, J. H., and K. D. Fink (1999). Numerical Methods Using Matlab. Third Edition, Prentice Hall.
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 numerical gradient.
gradient(F, h1 = 1, h2 = 1)
F 
vector of function values, or a matrix of values of a function of two variables. 
h1 
xcoordinates of grid points, or one value for the difference between grid points in xdirection. 
h2 
ycoordinates of grid points, or one value for the difference between grid points in ydirection. 
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 ydirection (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.
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 xdirection. 
Y 
numerical gradient/slope in xdirection. 
where each matrix is of the same size as F
.
TODO: If h2
is missing, it will not automatically be adapted.
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:
# Onedimensional 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()
# Twodimensional 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)
Modified GramSchmidt Process
gramSchmidt(A, tol = .Machine$double.eps^0.5)
A 
numeric matrix with 
tol 
numerical tolerance for being equal to zero. 
The modified GramSchmidt process uses the classical orthogonalization process to generate step by step an orthonoral basis of a vector space. The modified GramSchmidt iteration uses orthogonal projectors in order ro make the process numerically more stable.
List with two matrices Q
and R
, Q
orthonormal and
R
upper triangular, such that A=Q%*%R
.
Trefethen, L. N., and D. Bau III. (1997). Numerical Linear Algebra. SIAM, Society for Industrial and Applied Mathematics, Philadelphia.
## 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
Generate Hadamard matrix of a certain size.
hadamard(n)
n 
An integer of the form 2^e, 12*2^e, or 20*2^e 
An n
byn
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.
Matrix of size n
byn
of orthogonal columns consisting of
1 and 1 only.
Hadamard matrices have applications in combinatorics, signal processing, and numerical analysis.
hadamard(4)
H < hadamard(8)
t(H)
Finding roots of univariate functions using the Halley method.
halley(fun, x0, maxiter = 500, tol = 1e08, ...)
fun 
function whose root is to be found. 
x0 
starting value for the iteration. 
maxiter 
maximum number of iterations. 
tol 
absolute tolerance; default 
... 
additional arguments to be passed to the function. 
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.
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
https://mathworld.wolfram.com/HalleysMethod.html
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
Median absolute deviation (MAD) outlier in Time Series
hampel(x, k, t0 = 3)
x 
numeric vector representing a time series 
k 
window length 
t0 
threshold, default is 3 (Pearson's rule), see below. 
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.
Returning a list L
with L$y
the corrected time series and
L$ind
the indices of outliers in the ‘median absolut deviation’
sense.
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.
Pearson, R. K. (1999). “Data cleaning for dynamic modeling and control”. European Control Conference, ETH Zurich, Switzerland.
set.seed(8421)
x < numeric(1024)
z < rnorm(1024)
x[1] < z[1]
for (i in 2:1024) {
x[i] < 0.4*x[i1] + 0.8*x[i1]*z[i1] + 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)
Generate Hankel matrix from column and row vector
hankel(a, b)
a 
vector that will be the first column 
b 
vector that if present will form the last row. 
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.
matrix of size (length(a), length(b))
hankel(1:5, 5:1)
Hausdorff distance (aka Hausdorff dimension)
hausdorff_dist(P, Q)
P , Q

numerical matrices, representing points in an mdim. space. 
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 [ pq ] ]
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)).
A single scalar, the Hausdorff distance (dimension).
Barnsley, M. (1993). Fractals Everywhere. Morgan Kaufmann, San Francisco.
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 to calculate the arc distance between two points on earth (i.e., along a great circle).
haversine(loc1, loc2, R = 6371.0)
loc1 , loc2

Locations on earth; for format see Details. 
R 
Average earth radius R = 6371 km, can be changed on input. 
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 twovector 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.
Returns the distance in km.
Hans W. Borchers
Entry 'Great_circle_distance' in Wikipedia.
Implementations of the Haversine formula can also be found in other R packages, e.g. 'geoPlot' or 'geosphere'.
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 FrankfurtChicago is %8.3f km.\n', dis)
dis < haversine(fra, ord)
fprintf('Flight distance FrankfurtChicago is %8.3f km.\n', dis)
Generates the Hessenberg matrix for A.
hessenberg(A)
A 
square matrix 
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.
Returns a list with two elements,
H 
the upper Hessenberg Form of matrix A. 
H 
a unitary matrix. 
Press, Teukolsky, Vetterling, and Flannery (2007). Numerical Recipes: The Art of Scientific Computing. 3rd Edition, Cambridge University Press. (Section 11.6.2)
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
Numerically compute the Hessian matrix.
hessian(f, x0, h = .Machine$double.eps^(1/4), ...)
f 
univariate function of several variables. 
x0 
point in 
h 
step size. 
... 
variables to be passed to 
Computes the hessian matrix based on the threepoint central difference formula, expanded to two variables.
Assumes that the function has continuous partial derivatives.
An nbyn matrix with $\frac{\partial^2 f}{\partial x_i \partial x_j}$
as (i, j) entry.
Fausett, L. V. (2007). Applied Numerical Analysis Using Matlab. Second edition, Prentice Hall.
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)
Fast multiplication of Hessian and vector where computation of the full Hessian is not needed. Or determine the diagonal of the Hessian when nondiagonal entries are not needed or are nearly zero.
hessvec(f, x, v, csd = FALSE, ...)
hessdiag(f, x, ...)
f 
function whose hessian is to be computed. 
x 
point in 
v 
vector of length 
csd 
logocal, shall complexstep be applied. 
... 
more arguments to be passed to the function. 
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
complexstep 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.
hessvec
returns the product H(f,x) * v
as a vector.
hessdiag
returns the diagonal of the Hessian of f
.
B.A. Pearlmutter, Fast Exact Multiplication by the Hessian, Neural Computation (1994), Vol. 6, Issue 1, pp. 147160.
## 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^2b^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.48e05
max(abs(vex  v2)) # 7.15e05
max(abs(vex  v3)) # 0.09e05
max(abs(vex  v4)) # 2.46e05
## End(Not run)
Generate Hilbert matrix of dimension n
hilb(n)
n 
positive integer specifying the dimension of the Hilbert matrix 
Generate the Hilbert matrix H
of dimension n
with elements
H[i, j] = 1/(i+j1)
.
(Note: This matrix is illconditioned, see e.g. det(hilb(6))
.)
matrix of dimension n
hilb(5)
Histogramlike counting.
histc(x, edges)
x 
numeric vector or matrix. 
edges 
numeric vector of grid points, must be monotonically nondecreasing. 
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)
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.
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
Method for selecting the bin size of time histograms.
histss(x, n = 100, plotting = FALSE)
x 
numeric vector or matrix. 
n 
maximum number of bins. 
plotting 
logical; shall a histogram be plotted. 
Bin sizes of histograms are optimized in a way to best displays the underlying spike rate, for example in neurophysiological studies.
Returns the same list as the hist
function; the list is invisible
if the histogram is plotted.
Shimazaki H. and S. Shinomoto. A method for selecting the bin size of a time histogram. Neural Computation (2007) Vol. 19(6), 15031527
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)
An implementation of the HookeJeeves algorithm for derivativefree optimization.
hooke_jeeves(x0, fn, ..., lb = NULL, ub = NULL, tol = 1e08,
maxfeval = 10000, target = Inf, info = FALSE)
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. 
This method computes a new point using the values of f
at suitable
points along the orthogonal coordinate directions around the last point.
List with following components:
xmin 
minimum solution found so far. 
fmin 
value of 
count 
number of function evaluations. 
convergence 
NOT USED at the moment. 
info 
special info from the solver. 
HookeJeeves is notorious for its number of function calls. Memoization is often suggested as a remedy.
For a similar implementation of HookeJeeves see the ‘dfoptim’ package.
C.T. Kelley (1999), Iterative Methods for Optimization, SIAM.
Quarteroni, Sacco, and Saleri (2007), Numerical Mathematics, SpringerVerlag.
## Rosenbrock function
rosenbrock < function(x) {
n < length(x)
x1 < x[2:n]
x2 < x[1:(n1)]
sum(100*(x1x2^2)^2 + (1x2)^2)
}
hooke_jeeves(c(0,0,0,0), rosenbrock)
## $xmin
## [1] 1.000002 1.000003 1.000007 1.000013
## $fmin
## [1] 5.849188e11
## $count
## [1] 1691
## $convergence
## [1] 0
## $info
## $info$solver
## [1] "HookeJeeves"
## $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] "HookeJeeves"
## $info$iterations
## [1] 26
Compute the value of a polynomial via Horner's Rule.
horner(p, x)
hornerdefl(p, x)
p 
Numeric vector representing a polynomial. 
x 
Numeric scalar, vector or matrix at which to evaluate the polynomial. 
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^{n1} + ... + 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.
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.
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
).
Quarteroni, A., and Saleri, F. (2006) Scientific Computing with Matlab and Octave. Second Edition, SpringerVerlag, Berlin Heidelberg.
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 and QR decomposition
householder(A)
A 
numeric matrix with 
The Householder method applies a succession of elementary unitary
matrices to the left of matrix A
. These matrices are the socalled
Householder reflections.
List with two matrices Q
and R
, Q
orthonormal and
R
upper triangular, such that A=Q%*%R
.
Trefethen, L. N., and D. Bau III. (1997). Numerical Linear Algebra. SIAM, Society for Industrial and Applied Mathematics, Philadelphia.
## 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 (n1):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.
humps(x)
sinc(x)
psinc(x, n)
x 
numeric scalar or vector. 
n 
positive integer. 
humps
is a test function for finding zeros, for optimization
and integration. Its root is at x = 1.2995
, a (local) minimum
at x = 0.6370
, and the integral from 0.5
to 1.0
is 8.0715
.
sinc
is defined as $sinc(t) = \frac{\sin(\pi t)}{\pi t}$
.
It is the continuous inverse Fourier transform of the rectangular pulse
of width $2\pi$
and height $1$
.
psinc
is the 'periodic sinc function' and is defined as
$psinc(x,n) = \frac{\sin(x n/2)}{n \sin(x/2)}$
.
Numeric scalar or vector.
## Not run:
plot(humps(), type="l"); grid()
x < seq(0, 10, length=101)
plot(x, sinc(x), type="l"); grid()
## End(Not run)
Calculates the Hurst exponent using R/S analysis.
hurstexp(x, d = 50, display = TRUE)
x 
a time series. 
d 
smallest box size; default 50. 
display 
logical; shall the results be printed to the console? 
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.
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
Derived from Matlab code of R. Weron, published on Matlab Central.
H.E. Hurst (1951) Longterm storage capacity of reservoirs, Transactions of the American Society of Civil Engineers 116, 770808.
R. Weron (2002) Estimating long range dependence: finite sample properties and confidence intervals, Physica A 312, 285299.
fractal::hurstSpec, RoverS, hurstBlock
and fArma::LrdModelling
## 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[i1] * (1  xlm[i1])
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)
Square root of sum of squares
hypot(x, y)
x , y

Vectors of real or complex numbers of the same size 
Elementbyelement computation of the square root of the sum of squares
of vectors resp. matrices x
and y
.
Returns a vector or matrix of the same size.
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.
hypot(3,4)
hypot(1, c(3, 4, 5))
hypot(c(0, 0), c(3, 4))
Performs the inverse Fast Fourier Transform.
ifft(x)
ifftshift(x)
fftshift(x)
x 
a real or complex vector 
ifft
returns the value of the normalized discrete, univariate,
inverse Fast Fourier Transform of the values in x
.
ifftshift
and fftshift
shift the zerocomponent to the center
of the spectrum, that is swap the left and right half of x
.
Real or complex vector of the same length.
Almost an alias for R's fft(x, inverse=TRUE)
, but dividing by
length(x)
.
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/21)) # 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)
Points inside polygon region.
inpolygon(x, y, xp, yp, boundary = FALSE)
x , y

x, ycoordinates 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. 
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.
Logical vector, the same length as x
.
Special care taken for points on the boundary.
Hormann, K., and A. Agathos (2001). The Point in Polygon Problem for Arbitrary Polygons. Computational Geometry, Vol. 20, No. 3, pp. 131–144.
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)
Combines several approaches to adaptive numerical integration of functions of one variable.
integral(fun, xmin, xmax,
method = c("Kronrod", "Clenshaw","Simpson"),
no_intervals = 8, random = FALSE,
reltol = 1e8, abstol = 0, ...)
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. 
integral
combines the following methods for adaptive
numerical integration (also available as separate functions):
Kronrod (GaussKronrod)
Clenshaw (ClenshawCurtis; not yet made adaptive)
Simpson (adaptive Simpson)
Recommended default method is GaussKronrod. Also try ClenshawCurtis 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 arrayvalued,
quadv
is called that applies an adaptive Simpson procedure,
other methods are ignored – not true anymore.]
Returns the integral, no error terms given.
integral
does not provide ‘new’ functionality, everything is
already contained in the functions called here. Other interesting
alternatives are GaussRichardson (quadgr
) and Romberg
(romberg
) integration.
Davis, Ph. J., and Ph. Rabinowitz (1984). Methods of Numerical Integration. Dover Publications, New York.
quadgk
, quadgr
, quadcc
,
simpadpt
, romberg
,
quadv
, quadinf
## 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 = 1e12, method = m)
cat(m, Q, abs(Qval), "\n")}
# Kron 1.582233 3.197442e13
# Rich 1.582233 3.197442e13 # use quadgr()
# Clen 1.582233 3.199663e13
# Simp 1.582233 3.241851e13
# Romb 1.582233 2.555733e13 # 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 = 1e12, method = m)
cat(m, Q, abs(Qval), "\n")}
# Kron 0.4989868 2.775558e16
# Rich 0.4989868 4.440892e16 # use quadgr()
# Clen 0.4989868 2.231548e14
# Simp 0.4989868 6.328271e15
# Romb 0.4989868 1.508793e13 # use romberg()
## Evaluate improper integral
fun < function(x) log(x)^2 * exp(x^2)
val < 1.9475221803007815976
Q < integral(fun, 0, Inf, reltol = 1e12)
# For infinite domains Gauss integration is applied!
cat(m, Q, abs(Qval), "\n")
# Kron 1.94752218028062 2.01587635473288e11
## 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 a double integral, resp. a triple integral by reducing it to a double integral.
integral2(fun, xmin, xmax, ymin, ymax, sector = FALSE,
reltol = 1e6, abstol = 0, maxlist = 5000,
singular = FALSE, vectorized = TRUE, ...)
integral3(fun, xmin, xmax, ymin, ymax, zmin, zmax,
reltol = 1e6, ...)
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. 
integral2
implements the ‘TwoD’ algorithm, that is GaussKronrod
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 constantsc <= 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
.
Returns a list with Q
the integral and error
the error term.
To avoid recursion, a possibly large matrix will be used and passed between subprograms. A more efficient implementation may be possible.
Copyright (c) 2008 Lawrence F. Shampine for Matlab code and description of the program; adapted and converted to R by Hans W Borchers.
Shampine, L. F. (2008). MATLAB Program for Quadrature in 2D. Proceedings of Applied Mathematics and Computation, 2008, pp. 266–274.
integral
, cubature:adaptIntegrate
fun < function(x, y) cos(x) * cos(y)
integral2(fun, 0, 1, 0, 1, reltol = 1e10)
# $Q: 0.708073418273571 # 0.70807341827357119350 = sin(1)^2
# $error: 8.618277e19 # 1.110223e16
## 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.800354e06
## Compute the volume over a sector
I < integral2(f, 0,pi/2, 0,1, sector = TRUE)
I$Q # 0.5236308  pi/6 => 3.203768e05
## 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.247091e08
## 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.998646e11
## 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(4x^2)
ymax < function(x) sqrt(4x^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
Onedimensional interpolation of points.
interp1(x, y, xi = x,
method = c("linear", "constant", "nearest", "spline", "cubic"))
x 
Numeric vector; points on the xaxis; at least two points require; will be sorted if necessary. 
y 
Numeric vector; values of the assumed underlying function;

xi 
Numeric vector; points at which to compute the interpolation;
all points must lie between 
method 
One of “constant", “linear", “nearest", “spline", or “cubic"; default is “linear" 
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 
Numeric vector representing values at points xi
.
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.
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)
Twodimensional data interpolation similar to a table lookup.
interp2(x, y, Z, xp, yp, method = c("linear", "nearest", "constant"))
x , y

vectors with monotonically increasing elements, representing
x and ycoordinates of the data values in 
Z 
numeric 
xp , yp

x, ycoordinates of points at which interpolated values will be computed. 
method 
interpolation method, “linear” the most useful. 
Computes a vector containing elements corresponding to the elements of
xp
and yp
, determining by interpolation within the
twodimensional 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.
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.
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.
interp1
, barylag2d
## 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)
Invert a numeric or complex matrix.
inv(a)
a 
real or complex square matrix 
Computes the matrix inverse by calling solve(a)
and catching the error
if the matrix is nearly singular.
square matrix that is the inverse of a
.
inv()
is the function name used in Matlab/Octave.
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
Numerical inversion of Laplace transforms.
invlap(Fs, t1, t2, nnt, a = 6, ns = 20, nd = 19)
Fs 
function representing the function to be inversetransformed. 
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. 
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.
Returns a list with components x
the xcoordinates and y
the ycoordinates representing the original function in the interval
[t1,t2]
.
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.
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. 153166.
L.N.Trefethen, J.A.C.Weideman, and T.Schmelzer (2006). Talbot quadratures and rational approximations. BIT. Numerical Mathematics, 46(3):653–670.
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 # EulerMascheroni 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)
Determine if an object is empty.
isempty(x)
x 
an R object 
An empty object has length zero.
TRUE
if x
has length 0; otherwise, FALSE
.
isempty(c(0)) # FALSE
isempty(matrix(0, 1, 0)) # TRUE
Test for positive definiteness.
isposdef(A, psd = FALSE, tol = 1e10)
A 
symmetric matrix 
psd 
logical, shall semipositive definiteness be tested? 
tol 
tolerance to check symmetry and Cholesky decomposition. 
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 semipositive
definite. If not positive definite, still a warning will be generated.
Returns TRUE
or FALSE
.
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
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.
isprime(x)
x 
vector or matrix of nonnegative integers 
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.
array of elements 0, 1 with 1 indicating prime numbers
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 solutions of systems of linear equations.
itersolve(A, b, x0 = NULL, nmax = 1000, tol = .Machine$double.eps^(0.5),
method = c("GaussSeidel", "Jacobi", "Richardson"))
A 
numerical matrix, square and nonsingular. 
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, GaussSeidel, Jacobi, or Richardson. 
Iterative methods are based on splitting the matrix A=(PA)A
with a socalled ‘preconditioner’ matrix P. The methods differ in how
to choose this preconditioner.
Returns a list with components x
the solution, iter
the
number of iterations, and method
the name of the method applied.
Richardson's method allows to specify a ‘preconditioner’; this has not been implemented yet.
Quarteroni, A., and F. Saleri (2006). Scientific Computing with MATLAB and Octave. SpringerVerlag, Berlin Heidelberg.
N < 10
A < Diag(rep(3,N)) + Diag(rep(2, N1), k=1) + Diag(rep(1, N1), k=1)
b < A %*% rep(1, N)
x0 < rep(0, N)
itersolve(A, b, tol = 1e8, method = "GaussSeidel")
# [1] 1 1 1 1 1 1 1 1 1 1
# [1] 87
itersolve(A, b, x0 = 1:10, tol = 1e8, method = "Jacobi")
# [1] 1 1 1 1 1 1 1 1 1 1
# [1] 177
Jacobian matrix of a function R^n –> R^m .
jacobian(f, x0, heps = .Machine$double.eps^(1/3), ...)
f 

x0 
Numeric vector of length 
heps 
This is 
... 
parameters to be passed to f. 
Computes the derivative of each funktion $f_j$
by variable $x_i$
separately, taking the discrete step $h$
.
Numeric m
byn
matrix J
where the entry J[j, i]
is $\frac{\partial f_j}{\partial x_i}$
, i.e. the derivatives of function
$f_j$
line up in row $i$
for $x_1, \ldots, x_n$
.
Obviously, this function is not vectorized.
Quarteroni, A., R. Sacco, and F. Saleri (2007). Numerical Mathematics. Second Edition, SpringerVerlag, Berlin Heidelberg.
gradient
## 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
Simple and ordinary Kriging interpolation and interpolating function.
kriging(u, v, u0, type = c("ordinary", "simple"))
u 
an 
v 
an 
u0 
a 
type 
character; values ‘simple’ or ‘ordinary’; no partial matching. 
Kriging is a geospatial estimation procedure that estimates points based on the variations of known points in a nonregular grid. It is especially suited for surfaces.
kriging
returns a k
dim. vektor of interpolation values.
In the literature, different versions and extensions are discussed.
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. 144147.
akimaInterp
, barylag2d
, package kriging
## 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(n4, 1, 1)) # add four vertices
y < c(1, 1, 1, 1, runif(n4, 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