##' Differentiate a rational function
##'
##' Calculate the derivative of a rational function. The returned value
##' result is still an object of class "rationalfun".
##' @S3method deriv rationalfun
##' @method deriv rationalfun
##' @export
##' @param expr an object of class "rationalfun"
##' @param \dots not used in this function
##' @return An object of class "rationalfun" representing
##' the derivative of the original rational function.
##' @seealso \code{\link[polynom]{deriv.polynomial}},
##' \code{\link[stats]{deriv}}
##' @keywords symbolmath
##' @examples # (x + 1) / (x^2 + x + 1)
##' r <- rationalfun(c(1, 1), c(1, 1, 1))
##' deriv(r)
deriv.rationalfun <- function(expr, ...)
{
expr <- simplify(expr);
pp <- expr$numerator;
qq <- expr$denominator;
numer <- deriv(pp) * qq - pp * deriv(qq);
denom <- qq * qq;
return(rationalfun.poly(numer, denom));
}
integraterf.transcendent <- function(x, ...)
{
numer <- x$numerator;
denom <- x$denominator;
numercoef <- coef(numer);
if(numercoef %*% numercoef < 1e-20) return(substitute(0 * x));
ddenom <- deriv(denom);
coefrf <- rationalfun(coef(numer), coef(ddenom));
fac <- solve(denom);
fac.re <- Re(fac);
fac.im <- Im(fac);
real.index <- abs(fac.im) < 1e-10;
real.roots <- fac.re[real.index];
fac.re <- fac.re[!real.index];
fac.im <- fac.im[!real.index];
conju.index <- duplicated(fac.re) & duplicated(abs(fac.im));
conju.roots.re <- fac.re[conju.index];
conju.roots.im <- fac.im[conju.index];
conju.roots.inner <- conju.roots.re^2 + conju.roots.im^2;
real.coefs <- predict(coefrf, real.roots);
conju.coefs <- predict(coefrf, complex(real = conju.roots.re,
imaginary = conju.roots.im));
conju.coefs.re <- Re(conju.coefs);
conju.coefs.im <- Im(conju.coefs);
expr <- NULL;
if(length(real.coefs))
{
term0 <- substitute(log(abs(ex)),
list(ex = .linear(real.roots[1])));
term <- call("*", real.coefs[1], term0);
expr <- term;
for(i in seq_along(real.coefs)[-1])
{
term0 <- substitute(log(abs(ex)),
list(ex = .linear(real.roots[i])));
term <- call("*", abs(real.coefs[i]), term0);
op <- if(real.coefs[i] > 0) "+" else "-";
expr <- call(op, expr, term);
}
}
if(length(conju.coefs))
{
term01 <- substitute(log(ex),
list(ex = .quadratic(-2 * conju.roots.re[1],
conju.roots.inner[1])));
if(is.null(expr))
{
expr <- call("*", conju.coefs.re[1], term01);
} else {
term1 <- call("*", abs(conju.coefs.re[1]), term01);
op <- if(conju.coefs.re[1] > 0) "+" else "-";
expr <- call(op, expr, term1);
}
term02 <- substitute(atan(ex),
list(ex = .frac(conju.roots.re[1],
conju.roots.im[1])));
term2 <- call("*", abs(2 * conju.coefs.im[1]), term02);
op <- if(conju.coefs.im[1] > 0) "-" else "+";
expr <- call(op, expr, term2);
for(i in seq_along(conju.coefs)[-1])
{
term01 <- substitute(log(ex),
list(ex = .quadratic(-2 * conju.roots.re[i],
conju.roots.inner[i])));
term1 <- call("*", abs(conju.coefs.re[i]), term01);
op <- if(conju.coefs.re[i] > 0) "+" else "-";
expr <- call(op, expr, term1);
term02 <- substitute(atan(ex),
list(ex = .frac(conju.roots.re[i],
conju.roots.im[i])));
term2 <- call("*", abs(2 * conju.coefs.im[i]), term02);
op <- if(conju.coefs.im[i] > 0) "-" else "+";
expr <- call(op, expr, term2);
}
}
return(expr);
}
##' Integrate a rational function
##'
##' Calculate the integral of a rational function. See "Details".
##'
##' The returned value is a function call with argument named "x".
##' That is, the integral is an expression in R with an explicit form,
##' which could be evaluated directly by calling \code{\link{eval}()},
##' or indirectly using the \code{\link{int2fun}()} function.
##'
##' The algorithm is based on the Hermite-Ostrogradski formula which is
##' discussed in the reference. See the article for more details.
##' @S3method integral rationalfun
##' @method integral rationalfun
##' @export
##' @param expr an object of class "rationalfun"
##' @param \dots not used in this function
##' @return A function call representing the explicit form of the
##' integral.
##' @seealso \code{\link[polynom]{integral.polynomial}}
##' @references T. N. Subramaniam, and Donald E. G. Malm,
##' How to Integrate Rational Functions,
##' \emph{The American Mathematical Monthly},
##' Vol. 99, No.8 (1992), 762-772.
##' @keywords symbolmath
##' @examples # (x + 1) / (x^2 + x + 1)
##' r <- rationalfun(c(1, 1), c(1, 1, 1))
##' expr <- integral(r)
##' # Evaluate the call directly
##' eval(expr, list(x = 2))
##' # Use int2fun()
##' f <- int2fun(expr)
##' f(2)
integral.rationalfun <- function(expr, ...)
{
x <- expr;
pp <- x$numerator;
qq <- x$denominator;
poly.part <- pp / qq;
poly.expr <- parse(text = as.character(poly.part))[[1]];
pp <- pp %% qq;
qq1 <- .GCD(qq, deriv(qq));
if(.degree(qq1) == 0)
{
x <- rationalfun.poly(pp, qq);
x <- simplify(x);
transcendent.expr <- integraterf.transcendent(x);
expr <- call("+", transcendent.expr, poly.expr);
return(expr);
}
qq2 <- qq / qq1;
ss <- deriv(qq1) * qq2 / qq1;
# P = P1' * Q2 - P1 * S + P2 * Q1
np <- .degree(qq2);
nq <- .degree(qq1);
xx <- seq(0.295782039048272, 5.277857726515673, length.out = np + nq);
p1p <- outer(xx, 0:(nq - 1), function(x, k) k * x^(k - 1));
term1 <- sweep(p1p, 1, predict(qq2, xx), "*");
p1 <- outer(xx, 0:(nq - 1), function(x, k) x^k);
term2 <- sweep(p1, 1, predict(ss, xx), "*");
p2 <- outer(xx, 0:(np - 1), function(x, k) x^k);
term3 <- sweep(p2, 1, predict(qq1, xx), "*");
A <- cbind(term1 - term2, term3);
coefs <- solve(A, predict(pp, xx));
p1.coef <- coefs[1:nq];
p2.coef <- coefs[-(1:nq)];
rational.part <- rationalfun(p1.coef, coef(qq1));
rational.part <- simplify(rational.part);
rational.expr <- parse(text = as.character(rational.part))[[1]];
transcendent.part <- rationalfun(p2.coef, coef(qq2));
transcendent.part <- simplify(transcendent.part);
transcendent.expr <- integraterf.transcendent(transcendent.part);
expr <- call("+", transcendent.expr, rational.expr);
expr <- call("+", expr, poly.expr);
return(expr);
}
##' Convert a call to a function
##'
##' Convert a function call to a function in R. In this package, the
##' function is typically used to convert the result of
##' \code{\link{integral.rationalfun}()} to a function with one
##' argument.
##' @export
##' @param expr a function call, typically returned by
##' \code{\link{integral.rationalfun}()}.
##' @return A function with one argument which could be a real
##' or complex vector.
##' @seealso \code{\link[polynom]{integral.polynomial}}
##' @keywords symbolmath
##' @examples x <- rationalfun(c(-6, -1, -8, 15, -1, 8, -9, 2),
##' c(8, 12, 16, 4, 4))
##' int <- integral(x)
##' fun <- int2fun(int)
##' fun(c(0, 1))
int2fun <- function(expr)
{
f <- function(x) NULL;
body(f) <- expr;
return(f);
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.