#' @title Basis functions for the quintic spline
#'
#' @author Thomas Blanchet, Juliette Fournier, Thomas Piketty
#'
#' @description Basis functions and their derivatives for the quintic Hermite
#' spline.
#'
#' @param x A value in [0, 1], at which the functions are evaluated.
#'
#' @return The value of the function.
#'
#' @rdname basis_functions_quintic_spline
#' @name basis_functions_quintic_spline
#' @aliases h00 h01 h10 h11 h20 h21 h00d1 h01d1 h10d1 h11d1 h20d1 h21d1
#' h00d2 h01d2 h10d2 h11d2 h20d2 h21d2
# Basis functions
h00 <- function(x) {
return(1 - 10*x^3 + 15*x^4 - 6*x^5)
}
h01 <- function(x) {
return(10*x^3 - 15*x^4 + 6*x^5)
}
h10 <- function(x) {
return(x - 6*x^3 + 8*x^4 - 3*x^5)
}
h11 <- function(x) {
return(-4*x^3 + 7*x^4 - 3*x^5)
}
h20 <- function(x) {
return(x^2/2 - (3*x^3)/2 + (3*x^4)/2 - x^5/2)
}
h21 <- function(x) {
return(x^3/2 - x^4 + x^5/2)
}
# First derivatives
h00d1 <- function(x) {
return(-30*x^2 + 60*x^3 - 30*x^4)
}
h01d1 <- function(x) {
return(30*x^2 - 60*x^3 + 30*x^4)
}
h10d1 <- function(x) {
return(1 - 18*x^2 + 32*x^3 - 15*x^4)
}
h11d1 <- function(x) {
return(-12*x^2 + 28*x^3 - 15*x^4)
}
h20d1 <- function(x) {
return(x - (9*x^2)/2 + 6*x^3 - (5*x^4)/2)
}
h21d1 <- function(x) {
return((3*x^2)/2 - 4*x^3 + (5*x^4)/2)
}
# Second derivatives
h00d2 <- function(x) {
return(-60*x + 180*x^2 - 120*x^3)
}
h01d2 <- function(x) {
return(60*x - 180*x^2 + 120*x^3)
}
h10d2 <- function(x) {
return(-36*x + 96*x^2 - 60*x^3)
}
h11d2 <- function(x) {
return(-24*x + 84*x^2 - 60*x^3)
}
h20d2 <- function(x) {
return(1 - 9*x + 18*x^2 - 10*x^3)
}
h21d2 <- function(x) {
return(3*x - 12*x^2 + 10*x^3)
}
# Third derivatives
h00d3 <- function(x) {
return(-60 + 360*x - 360*x^2)
}
h01d3 <- function(x) {
return(60 - 360*x + 360*x^2)
}
h10d3 <- function(x) {
return(-36 + 192*x - 180*x^2)
}
h11d3 <- function(x) {
return(-24 + 168*x - 180*x^2)
}
h20d3 <- function(x) {
return(-9 + 36*x - 30*x^2)
}
h21d3 <- function(x) {
return(3 - 24*x + 30*x^2)
}
#' @title Quintic spline function
#'
#' @author Thomas Blanchet, Juliette Fournier, Thomas Piketty
#'
#' @description Quintic spline over an interval \code{[x0, x1]}. It is such that:
#' \itemize{
#' \item{\code{h(x0) == y0}}
#' \item{\code{h(x1) == y1}}
#' \item{\code{h’(x0) == s0}}
#' \item{\code{h’(x1) == s1}}
#' \item{\code{h’’(x0) == a0}}
#' \item{\code{h’’(x1) == a1}}
#' }
#'
#' @param x A number between \code{x0} and \code{x1}, at which the spline is
#' evaluated.
#' @param x0 Lower bound of the interval.
#' @param x1 Upper bound of the interval.
#' @param y0 Value of the spline at \code{x0}.
#' @param y1 Value of the spline at \code{x1}.
#' @param s0 Slope of the spline at \code{x0}.
#' @param s1 Slope of the spline at \code{x1}.
#' @param a0 Second derivative of the spline at \code{x0}.
#' @param a1 Second derivative of the spline at \code{x1}.
#'
#' @return The value of the interpolation at \code{x}.
h <- function(x, x0, x1, y0, y1, s0, s1, a0, a1) {
# Normalize x in the [0, 1] interval
x <- (x - x0)/(x1 - x0)
return(
y0*h00(x) +
y1*h01(x) +
s0*(x1 - x0)*h10(x) +
s1*(x1 - x0)*h11(x) +
a0*(x1 - x0)^2*h20(x) +
a1*(x1 - x0)^2*h21(x)
)
}
#' @title First derivative of the quintic spline function
#'
#' @author Thomas Blanchet, Juliette Fournier, Thomas Piketty
#'
#' @description First derivative of the quintic spline over an interval
#' \code{[x0, x1]}.
#'
#' @param x A number between \code{x0} and \code{x1}, at which the function is
#' evaluated.
#' @param x0 Lower bound of the interval.
#' @param x1 Upper bound of the interval.
#' @param y0 Value of the spline at \code{x0}.
#' @param y1 Value of the spline at \code{x1}.
#' @param s0 Slope of the spline at \code{x0}.
#' @param s1 Slope of the spline at \code{x1}.
#' @param a0 Second derivative of the spline at \code{x0}.
#' @param a1 Second derivative of the spline at \code{x1}.
#'
#' @return The value of first derivative of the interpolation function at \code{x}.
hd1 <- function(x, x0, x1, y0, y1, s0, s1, a0, a1) {
# Normalize x in the [0, 1] interval
x <- (x - x0)/(x1 - x0)
return(
y0*h00d1(x)/(x1 - x0) +
y1*h01d1(x)/(x1 - x0) +
s0*h10d1(x) +
s1*h11d1(x) +
a0*(x1 - x0)*h20d1(x) +
a1*(x1 - x0)*h21d1(x)
)
}
#' @title Second derivative of the quintic spline function
#'
#' @author Thomas Blanchet, Juliette Fournier, Thomas Piketty
#'
#' @description Second derivative of the quintic spline over an interval
#' \code{[x0, x1]}.
#'
#' @param x A number between \code{x0} and \code{x1}, at which the function is
#' evaluated.
#' @param x0 Lower bound of the interval.
#' @param x1 Upper bound of the interval.
#' @param y0 Value of the spline at \code{x0}.
#' @param y1 Value of the spline at \code{x1}.
#' @param s0 Slope of the spline at \code{x0}.
#' @param s1 Slope of the spline at \code{x1}.
#' @param a0 Second derivative of the spline at \code{x0}.
#' @param a1 Second derivative of the spline at \code{x1}.
#'
#' @return The value of second derivative of the interpolation function at \code{x}.
hd2 <- function(x, x0, x1, y0, y1, s0, s1, a0, a1) {
# Normalize x in the [0, 1] interval
x <- (x - x0)/(x1 - x0)
return(
y0*h00d2(x)/(x1 - x0)^2 +
y1*h01d2(x)/(x1 - x0)^2 +
s0*h10d2(x)/(x1 - x0) +
s1*h11d2(x)/(x1 - x0) +
a0*h20d2(x) +
a1*h21d2(x)
)
}
#' @title Third derivative of the quintic spline function
#'
#' @author Thomas Blanchet, Juliette Fournier, Thomas Piketty
#'
#' @description Third derivative of the quintic spline over an interval
#' \code{[x0, x1]}.
#'
#' @param x A number between \code{x0} and \code{x1}, at which the function is
#' evaluated.
#' @param x0 Lower bound of the interval.
#' @param x1 Upper bound of the interval.
#' @param y0 Value of the spline at \code{x0}.
#' @param y1 Value of the spline at \code{x1}.
#' @param s0 Slope of the spline at \code{x0}.
#' @param s1 Slope of the spline at \code{x1}.
#' @param a0 Second derivative of the spline at \code{x0}.
#' @param a1 Second derivative of the spline at \code{x1}.
#'
#' @return The value of third derivative of the interpolation function at \code{x}.
hd3 <- function(x, x0, x1, y0, y1, s0, s1, a0, a1) {
# Normalize x in the [0, 1] interval
x <- (x - x0)/(x1 - x0)
return(
y0*h00d3(x)/(x1 - x0)^3 +
y1*h01d3(x)/(x1 - x0)^3 +
s0*h10d3(x)/(x1 - x0)^2 +
s1*h11d3(x)/(x1 - x0)^2 +
a0*h20d3(x)/(x1 - x0) +
a1*h21d3(x)/(x1 - x0)
)
}
#' @title Quintic spline function
#'
#' @author Thomas Blanchet, Juliette Fournier, Thomas Piketty
#'
#' @description Quintic spline interpolation function over several intervals.
#'
#' @param x A vector of values between \code{min(xk)} and \code{max(xk)},
#' at which the interpolation function is evaluated.
#' @param xk A vector of interpolation points.
#' @param yk A vector of values at each interpolation point.
#' @param sk A vector of slopes at each interpolation point.
#' @param ak A vector of second derivatives at each interpolation point.
#'
#' @return The values of the interpolation function at each point \code{x}.
#'
#' @export
quintic_spline <- function(x, xk, yk, sk, ak) {
# Find the interval in which each value of x falls
k <- cut(x, breaks=xk, include.lowest=TRUE, labels=FALSE)
return(h(x,
xk[k], xk[k + 1],
yk[k], yk[k + 1],
sk[k], sk[k + 1],
ak[k], ak[k + 1]
))
}
#' @title First derivative of the quintic spline function
#'
#' @author Thomas Blanchet, Juliette Fournier, Thomas Piketty
#'
#' @description First derivative of the quintic spline interpolation function
#' over several intervals.
#'
#' @param x A vector of values between \code{min(xk)} and \code{max(xk)},
#' at which the function is evaluated.
#' @param xk A vector of interpolation points.
#' @param yk A vector of values at each interpolation point.
#' @param sk A vector of slopes at each interpolation point.
#' @param ak A vector of second derivatives at each interpolation point.
#'
#' @return The values of the first derivative of the interpolation function at
#' each point \code{x}.
#'
#' @export
deriv_quintic_spline <- function(x, xk, yk, sk, ak) {
# Find the interval in which each value of x falls
k <- cut(x, breaks=xk, include.lowest=TRUE, labels=FALSE)
return(hd1(x,
xk[k], xk[k + 1],
yk[k], yk[k + 1],
sk[k], sk[k + 1],
ak[k], ak[k + 1]
))
}
#' @title Second derivative of the quintic spline function
#'
#' @author Thomas Blanchet, Juliette Fournier, Thomas Piketty
#'
#' @description Second derivative of the quintic spline interpolation function
#' over several intervals.
#'
#' @param x A vector of values between \code{min(xk)} and \code{max(xk)},
#' at which the function is evaluated.
#' @param xk A vector of interpolation points.
#' @param yk A vector of values at each interpolation point.
#' @param sk A vector of slopes at each interpolation point.
#' @param ak A vector of second derivatives at each interpolation point.
#'
#' @return The values of the second derivative of the interpolation function at
#' each point \code{x}.
deriv2_quintic_spline <- function(x, xk, yk, sk, ak) {
# Find the interval in which each value of x falls
k <- cut(x, breaks=xk, include.lowest=TRUE, labels=FALSE)
return(hd2(x,
xk[k], xk[k + 1],
yk[k], yk[k + 1],
sk[k], sk[k + 1],
ak[k], ak[k + 1]
))
}
#' @title Third derivative of the quintic spline function
#'
#' @author Thomas Blanchet, Juliette Fournier, Thomas Piketty
#'
#' @description Third derivative of the quintic spline interpolation function
#' over several intervals.
#'
#' @param x A vector of values between \code{min(xk)} and \code{max(xk)},
#' at which the function is evaluated.
#' @param xk A vector of interpolation points.
#' @param yk A vector of values at each interpolation point.
#' @param sk A vector of slopes at each interpolation point.
#' @param ak A vector of second derivatives at each interpolation point.
#'
#' @return The values of the third derivative of the interpolation function at
#' each point \code{x}.
deriv3_quintic_spline <- function(x, xk, yk, sk, ak) {
# Find the interval in which each value of x falls
k <- cut(x, breaks=xk, include.lowest=TRUE, labels=FALSE)
return(hd3(x,
xk[k], xk[k + 1],
yk[k], yk[k + 1],
sk[k], sk[k + 1],
ak[k], ak[k + 1]
))
}
#' @title Tension of several splines
#'
#' @author Thomas Blanchet, Juliette Fournier, Thomas Piketty
#'
#' @description Tension (as mesured by the integral of the square of the
#' third derivative) of several splines put together.
#'
#' @param xk A vector of interpolation points.
#' @param yk A vector of values at each interpolation point.
#' @param sk A vector of slopes at each interpolation point.
#' @param ak A vector of second derivatives at each interpolation point.
#'
#' @return The tension of the spline over [\code{min(xk)}, \code{max(xk)}].
tension <- function(xk, yk, sk, ak) {
n <- length(xk)
return(sum(tension_spline(
xk[1:(n - 1)], xk[2:n],
yk[1:(n - 1)], yk[2:n],
sk[1:(n - 1)], sk[2:n],
ak[1:(n - 1)], ak[2:n]
)))
}
#' @title Gradient of the tension of several splines
#'
#' @author Thomas Blanchet, Juliette Fournier, Thomas Piketty
#'
#' @description The gradient of the tension of several splines with respect
#' to the value of the second derivative at each knot.
#'
#' @param xk A vector of interpolation points.
#' @param yk A vector of values at each interpolation point.
#' @param sk A vector of slopes at each interpolation point.
#' @param ak A vector of second derivatives at each interpolation point.
#'
#' @return A vector of length \code{length(ak)} with the derivatives.
tension_grad <- function(xk, yk, sk, ak) {
n <- length(xk)
grad_a0k <- tension_spline_deriv_a0(
xk[1:(n - 1)], xk[2:n],
yk[1:(n - 1)], yk[2:n],
sk[1:(n - 1)], sk[2:n],
ak[1:(n - 1)], ak[2:n]
)
grad_a1k <- tension_spline_deriv_a1(
xk[1:(n - 1)], xk[2:n],
yk[1:(n - 1)], yk[2:n],
sk[1:(n - 1)], sk[2:n],
ak[1:(n - 1)], ak[2:n]
)
return(c(grad_a0k, 0) + c(0, grad_a1k))
}
#' @title Estimate natural quintic spline with known first derivatives
#'
#' @author Thomas Blanchet, Juliette Fournier, Thomas Piketty
#'
#' @description Estimate a natural quintic spline from known values of the
#' function and its derivatives: that is, estimate the second derivative of
#' the spline at each knot to ensure a continuous third derivative and
#' a zero third derivative at the extremities.
#'
#' @param xk A vector of interpolation points.
#' @param yk A vector of values at each interpolation point.
#' @param sk A vector of slopes at each interpolation point.
#'
#' @return The vector \code{ak} of second derivatives at each knot.
#'
#' @importFrom MASS ginv
natural_quintic_spline <- function(xk, yk, sk) {
n <- length(xk)
x1 <- xk[1:(n - 2)]
x2 <- xk[2:(n - 1)]
x3 <- xk[3:(n - 0)]
y1 <- yk[1:(n - 2)]
y2 <- yk[2:(n - 1)]
y3 <- yk[3:(n - 0)]
s1 <- sk[1:(n - 2)]
s2 <- sk[2:(n - 1)]
s3 <- sk[3:(n - 0)]
# Tridiagonal matrix representing the system of equations
A <- matrix(data=0, nrow=n, ncol=n)
A[cbind(1:(n - 1), 2:n)] <- c(
-3/(xk[2] - xk[1]),
-3/(x3 - x2)
)
A[cbind(1:n, 1:n)] <- c(
9/(xk[2] - xk[1]),
9/(x2 - x1) + 9/(x3 - x2),
9/(xk[n] - xk[n - 1])
)
A[cbind(2:n, 1:(n - 1))] <- c(
-3/(x2 - x1),
-3/(xk[n] - xk[n - 1])
)
b <- c(
- yk[1]*60/(xk[2] - xk[1])^3
+ yk[2]*60/(xk[2] - xk[1])^3
- sk[1]*36/(xk[2] - xk[1])^2
- sk[2]*24/(xk[2] - xk[1])^2,
+ y1*60/(x2 - x1)^3
- y2*60/(x2 - x1)^3
- y2*60/(x3 - x2)^3
+ y3*60/(x3 - x2)^3
+ s1*24/(x2 - x1)^2
+ s2*36/(x2 - x1)^2
- s2*36/(x3 - x2)^2
- s3*24/(x3 - x2)^2,
+ yk[n - 1]*60/(xk[n] - xk[n - 1])^3
- yk[n]*60/(xk[n] - xk[n - 1])^3
+ sk[n - 1]*24/(xk[n] - xk[n - 1])^2
+ sk[n]*36/(xk[n] - xk[n - 1])^2
)
# Solve and return solution
return(as.vector(ginv(A) %*% b))
}
#' @title Estimate natural quintic spline with unknown first derivatives
#'
#' @author Thomas Blanchet, Juliette Fournier, Thomas Piketty
#'
#' @description Estimate a natural quintic spline from known values of the
#' function only. That is, it estimates the first and second derivatives at
#' each point by ensuring that the third and fourth derivatives are
#' continuous at each knot, and equal to zero at the first and the last point.
#'
#' @param xk A vector of interpolation points.
#' @param yk A vector of values at each interpolation point.
#'
#' @return A list with two components:
#' \itemize{
#' \item{\code{sk}} A vector of first derivatives at each point.
#' \item{\code{ak}} A vector of second derivatives at each point.
#' }
#'
#' @importFrom MASS ginv
natural_quintic_spline_noderiv <- function(xk, yk) {
n <- length(xk)
x1 <- xk[1:(n - 2)]
x2 <- xk[2:(n - 1)]
x3 <- xk[3:(n - 0)]
y1 <- yk[1:(n - 2)]
y2 <- yk[2:(n - 1)]
y3 <- yk[3:(n - 0)]
# Matrix representing the system of equations
A <- matrix(data=0, nrow=2*n, ncol=2*n)
# sk, third derivative
A[cbind(1:(n - 1), 2:n)] <- c(
-24/(xk[2] - xk[1])^2,
24/(x3 - x2)^2
)
A[cbind(1:n, 1:n)] <- c(
-36/(xk[2] - xk[1])^2,
36/(x3 - x2)^2 - 36/(x2 - x1)^2,
-36/(xk[n] - xk[n - 1])^2
)
A[cbind(2:n, 1:(n - 1))] <- c(
-24/(x2 - x1)^2,
-24/(xk[n] - xk[n - 1])^2
)
# ak, third derivative
A[cbind(1:(n - 1), (n + 2):(2*n))] <- c(
3/(xk[2] - xk[1]),
-3/(x3 - x2)
)
A[cbind(1:n, (n + 1):(2*n))] <- c(
-9/(xk[2] - xk[1]),
9/(x2 - x1) + 9/(x3 - x2),
9/(xk[n] - xk[n - 1])
)
A[cbind(2:n, (n + 1):(2*n - 1))] <- c(
-3/(x2 - x1),
-3/(xk[n] - xk[n - 1])
)
# sk, fourth derivative
A[cbind((n + 1):(2*n - 1), 2:n)] <- c(
168/(xk[2] - xk[1])^3,
-168/(x3 - x2)^3
)
A[cbind((n + 1):(2*n), 1:n)] <- c(
192/(xk[2] - xk[1])^3,
-192/(x2 - x1)^3 - 192/(x3 - x2)^3,
-192/(xk[n] - xk[n - 1])^3
)
A[cbind((n + 2):(2*n), 1:(n - 1))] <- c(
-168/(x2 - x1)^3,
-168/(xk[n] - xk[n - 1])^3
)
# ak, fourth derivative
A[cbind((n + 1):(2*n - 1), (n + 2):(2*n))] <- c(
-24/(xk[2] - xk[1])^2,
24/(x3 - x2)^2
)
A[cbind((n + 1):(2*n), (n + 1):(2*n))] <- c(
36/(xk[2] - xk[1])^2,
36/(x2 - x1)^2 - 36/(x3 - x2)^2,
36/(xk[n] - xk[n - 1])^2
)
A[cbind((n + 2):(2*n), (n + 1):(2*n - 1))] <- c(
-24/(x2 - x1)^2,
-24/(xk[n] - xk[n - 1])^2
)
b <- c(
-60*(yk[2] - yk[1])/(xk[2] - xk[1])^3,
(60*y1)/(x2 - x1)^3 - (60/(x2 - x1)^3 + 60/(x3 - x2)^3)*y2 + (60*y3)/(x3 - x2)^3,
-60*(yk[n] - yk[n - 1])/(xk[n] - xk[n - 1])^3,
360*(yk[2] - yk[1])/(xk[2] - xk[1])^4,
(360*y1)/(x2 - x1)^4 - (360/(x2 - x1)^4 - 360/(x3 - x2)^4)*y2 - (360*y3)/(x3 - x2)^4,
-360*(yk[n] - yk[n - 1])/(xk[n] - xk[n - 1])^4
)
v <- as.vector(ginv(A) %*% b)
# Solve and return solution
return(list(sk=v[1:n], ak=v[(n + 1):(2*n)]))
}
#' @title Estimate natural quintic spline with unknown first derivatives
#' with second derivative clamped at the last point
#'
#' @author Thomas Blanchet, Juliette Fournier, Thomas Piketty
#'
#' @description Estimate a natural quintic spline from known values of the
#' function and the value of the first derivative at the first knot. That is,
#' it estimates the first and second derivatives at each point by ensuring
#' that the third and fourth derivatives are continuous at each knot.
#'
#' @param xk A vector of interpolation points.
#' @param yk A vector of values at each interpolation point.
#' @param an The value of the second derivative at the last point.
#'
#' @return A list with two components:
#' \itemize{
#' \item{\code{sk}} A vector of first derivatives at each point.
#' \item{\code{ak}} A vector of second derivatives at each point.
#' }
#'
#' @importFrom MASS ginv
clamped1_quintic_spline_noderiv <- function(xk, yk, an) {
n <- length(xk)
x1 <- xk[1:(n - 2)]
x2 <- xk[2:(n - 1)]
x3 <- xk[3:(n - 0)]
y1 <- yk[1:(n - 2)]
y2 <- yk[2:(n - 1)]
y3 <- yk[3:(n - 0)]
# Matrix representing the system of equations
A <- matrix(data=0, nrow=2*n, ncol=2*n)
# sk, third derivative
A[cbind(1:(n - 1), 2:n)] <- c(
-24/(xk[2] - xk[1])^2,
24/(x3 - x2)^2
)
A[cbind(1:n, 1:n)] <- c(
-36/(xk[2] - xk[1])^2,
36/(x3 - x2)^2 - 36/(x2 - x1)^2,
-36/(xk[n] - xk[n - 1])^2
)
A[cbind(2:n, 1:(n - 1))] <- c(
-24/(x2 - x1)^2,
-24/(xk[n] - xk[n - 1])^2
)
# ak, third derivative
A[cbind(1:(n - 1), (n + 2):(2*n))] <- c(
3/(xk[2] - xk[1]),
-3/(x3 - x2)
)
A[cbind(1:n, (n + 1):(2*n))] <- c(
-9/(xk[2] - xk[1]),
9/(x2 - x1) + 9/(x3 - x2),
9/(xk[n] - xk[n - 1])
)
A[cbind(2:n, (n + 1):(2*n - 1))] <- c(
-3/(x2 - x1),
-3/(xk[n] - xk[n - 1])
)
# sk, fourth derivative
A[cbind((n + 1):(2*n - 1), 2:n)] <- c(
168/(xk[2] - xk[1])^3,
-168/(x3 - x2)^3
)
A[cbind((n + 1):(2*n), 1:n)] <- c(
192/(xk[2] - xk[1])^3,
-192/(x2 - x1)^3 - 192/(x3 - x2)^3,
-192/(xk[n] - xk[n - 1])^3
)
A[cbind((n + 2):(2*n), 1:(n - 1))] <- c(
-168/(x2 - x1)^3,
-168/(xk[n] - xk[n - 1])^3
)
# ak, fourth derivative
A[cbind((n + 1):(2*n - 1), (n + 2):(2*n))] <- c(
-24/(xk[2] - xk[1])^2,
24/(x3 - x2)^2
)
A[cbind((n + 1):(2*n), (n + 1):(2*n))] <- c(
36/(xk[2] - xk[1])^2,
36/(x2 - x1)^2 - 36/(x3 - x2)^2,
36/(xk[n] - xk[n - 1])^2
)
A[cbind((n + 2):(2*n), (n + 1):(2*n - 1))] <- c(
-24/(x2 - x1)^2,
-24/(xk[n] - xk[n - 1])^2
)
# Clamping
A[2*n, 1:(2*n)] <- c(rep(0, 2*n - 1), 1)
b <- c(
-60*(yk[2] - yk[1])/(xk[2] - xk[1])^3,
(60*y1)/(x2 - x1)^3 - (60/(x2 - x1)^3 + 60/(x3 - x2)^3)*y2 + (60*y3)/(x3 - x2)^3,
-60*(yk[n] - yk[n - 1])/(xk[n] - xk[n - 1])^3,
360*(yk[2] - yk[1])/(xk[2] - xk[1])^4,
(360*y1)/(x2 - x1)^4 - (360/(x2 - x1)^4 - 360/(x3 - x2)^4)*y2 - (360*y3)/(x3 - x2)^4,
an
)
v <- as.vector(ginv(A) %*% b)
# Solve and return solution
return(list(sk=v[1:n], ak=v[(n + 1):(2*n)]))
}
#' @title Estimate natural quintic spline with unknown first derivatives
#' with first derivative clamped at the first point and second derivative
#' clamped at the last
#'
#' @author Thomas Blanchet, Juliette Fournier, Thomas Piketty
#'
#' @description Estimate a natural quintic spline from known values of the
#' function and the value of the first derivative at the first knot. That is,
#' it estimates the first and second derivatives at each point (except the first,
#' for which only the second derivative is estimated) by ensuring
#' that the third and fourth derivatives are continuous at each knot, and equal
#' to zero at the first and the last point.
#'
#' @param xk A vector of interpolation points.
#' @param yk A vector of values at each interpolation point.
#' @param s1 The value of the first derivative at the first point.
#' @param an The value of the second derivative at the last point.
#'
#' @return A list with two components:
#' \itemize{
#' \item{\code{sk}} A vector of first derivatives at each point.
#' \item{\code{ak}} A vector of second derivatives at each point.
#' }
#'
#' @importFrom MASS ginv
clamped2_quintic_spline_noderiv <- function(xk, yk, s1, an) {
n <- length(xk)
x1 <- xk[1:(n - 2)]
x2 <- xk[2:(n - 1)]
x3 <- xk[3:(n - 0)]
y1 <- yk[1:(n - 2)]
y2 <- yk[2:(n - 1)]
y3 <- yk[3:(n - 0)]
# Matrix representing the system of equations
A <- matrix(data=0, nrow=2*n, ncol=2*n)
# sk, third derivative
A[cbind(1:(n - 1), 2:n)] <- c(
-24/(xk[2] - xk[1])^2,
24/(x3 - x2)^2
)
A[cbind(1:n, 1:n)] <- c(
-36/(xk[2] - xk[1])^2,
36/(x3 - x2)^2 - 36/(x2 - x1)^2,
-36/(xk[n] - xk[n - 1])^2
)
A[cbind(2:n, 1:(n - 1))] <- c(
-24/(x2 - x1)^2,
-24/(xk[n] - xk[n - 1])^2
)
# ak, third derivative
A[cbind(1:(n - 1), (n + 2):(2*n))] <- c(
3/(xk[2] - xk[1]),
-3/(x3 - x2)
)
A[cbind(1:n, (n + 1):(2*n))] <- c(
-9/(xk[2] - xk[1]),
9/(x2 - x1) + 9/(x3 - x2),
9/(xk[n] - xk[n - 1])
)
A[cbind(2:n, (n + 1):(2*n - 1))] <- c(
-3/(x2 - x1),
-3/(xk[n] - xk[n - 1])
)
# sk, fourth derivative
A[cbind((n + 1):(2*n - 1), 2:n)] <- c(
168/(xk[2] - xk[1])^3,
-168/(x3 - x2)^3
)
A[cbind((n + 1):(2*n), 1:n)] <- c(
192/(xk[2] - xk[1])^3,
-192/(x2 - x1)^3 - 192/(x3 - x2)^3,
-192/(xk[n] - xk[n - 1])^3
)
A[cbind((n + 2):(2*n), 1:(n - 1))] <- c(
-168/(x2 - x1)^3,
-168/(xk[n] - xk[n - 1])^3
)
# ak, fourth derivative
A[cbind((n + 1):(2*n - 1), (n + 2):(2*n))] <- c(
-24/(xk[2] - xk[1])^2,
24/(x3 - x2)^2
)
A[cbind((n + 1):(2*n), (n + 1):(2*n))] <- c(
36/(xk[2] - xk[1])^2,
36/(x2 - x1)^2 - 36/(x3 - x2)^2,
36/(xk[n] - xk[n - 1])^2
)
A[cbind((n + 2):(2*n), (n + 1):(2*n - 1))] <- c(
-24/(x2 - x1)^2,
-24/(xk[n] - xk[n - 1])^2
)
# Clamping
A[n + 1, 1:(2*n)] <- c(1, rep(0, 2*n - 1))
A[2*n, 1:(2*n)] <- c(rep(0, 2*n - 1), 1)
b <- c(
-60*(yk[2] - yk[1])/(xk[2] - xk[1])^3,
(60*y1)/(x2 - x1)^3 - (60/(x2 - x1)^3 + 60/(x3 - x2)^3)*y2 + (60*y3)/(x3 - x2)^3,
-60*(yk[n] - yk[n - 1])/(xk[n] - xk[n - 1])^3,
s1,
(360*y1)/(x2 - x1)^4 - (360/(x2 - x1)^4 - 360/(x3 - x2)^4)*y2 - (360*y3)/(x3 - x2)^4,
an
)
v <- as.vector(ginv(A) %*% b)
# Solve and return solution
return(list(sk=v[1:n], ak=v[(n + 1):(2*n)]))
}
#' @title Estimate a quintic spline known first derivatives with second
#' derivative clamped at the last point
#'
#' @author Thomas Blanchet, Juliette Fournier, Thomas Piketty
#'
#' @description Estimate a natural quintic spline from known values of the
#' function and its derivatives: that is, estimate the second derivative of
#' the spline at each knot to ensure a continuous third derivative and
#' a zero third derivative at the extremities.
#'
#' @param xk A vector of interpolation points.
#' @param yk A vector of values at each interpolation point.
#' @param sk A vector of slopes at each interpolation point.
#' @param an Second derivative at the last knot.
#'
#' @return The vector \code{ak} of second derivatives at each knot.
#'
#' @importFrom MASS ginv
clamped_quintic_spline <- function(xk, yk, sk, an) {
n <- length(xk)
x1 <- xk[1:(n - 2)]
x2 <- xk[2:(n - 1)]
x3 <- xk[3:(n - 0)]
y1 <- yk[1:(n - 2)]
y2 <- yk[2:(n - 1)]
y3 <- yk[3:(n - 0)]
s1 <- sk[1:(n - 2)]
s2 <- sk[2:(n - 1)]
s3 <- sk[3:(n - 0)]
# Tridiagonal matrix representing the system of equations
A <- matrix(data=0, nrow=n, ncol=n)
A[cbind(1:(n - 1), 2:n)] <- c(
-3/(xk[2] - xk[1]),
-3/(x3 - x2)
)
A[cbind(1:n, 1:n)] <- c(
9/(xk[2] - xk[1]),
9/(x2 - x1) + 9/(x3 - x2),
9/(xk[n] - xk[n - 1])
)
A[cbind(2:n, 1:(n - 1))] <- c(
-3/(x2 - x1),
-3/(xk[n] - xk[n - 1])
)
A[n, 1:n] <- c(rep(0, n - 1), 1)
b <- c(
- yk[1]*60/(xk[2] - xk[1])^3
+ yk[2]*60/(xk[2] - xk[1])^3
- sk[1]*36/(xk[2] - xk[1])^2
- sk[2]*24/(xk[2] - xk[1])^2,
+ y1*60/(x2 - x1)^3
- y2*60/(x2 - x1)^3
- y2*60/(x3 - x2)^3
+ y3*60/(x3 - x2)^3
+ s1*24/(x2 - x1)^2
+ s2*36/(x2 - x1)^2
- s2*36/(x3 - x2)^2
- s3*24/(x3 - x2)^2,
an
)
# Solve and return solution
return(as.vector(ginv(A) %*% b))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.