# Rational approximations of fractional operators

Source:`R/fractional.operators.R`

`fractional.operators.Rd`

`fractional.operators`

is used for computing an approximation,
which can be used for inference and simulation, of the fractional SPDE
$$L^\beta (\tau u(s)) = W.$$
Here \(L\) is a differential operator, \(\beta>0\) is
the fractional power, \(\tau\) is a positive scalar or vector that
scales the variance of the solution \(u\), and \(W\) is white noise.

## Arguments

- L
A finite element discretization of the operator \(L\).

- beta
The positive fractional power.

- C
The mass matrix of the finite element discretization.

- scale.factor
A constant \(c\) is a lower bound for the the smallest eigenvalue of the non-discretized operator \(L\).

- m
The order of the rational approximation, which needs to be a positive integer. The default value is 1. Higer values gives a more accurate approximation, which are more computationally expensive to use for inference. Currently, the largest value of m that is implemented is 4.

- tau
The constant or vector that scales the variance of the solution. The default value is 1.

## Value

`fractional.operators`

returns an object of class "rSPDEobj".
This object contains the following quantities:

- Pl
The operator \(P_l\).

- Pr
The operator \(P_r\).

- C
The mass lumped mass matrix.

- Ci
The inverse of

`C`

.- m
The order of the rational approximation.

- beta
The fractional power.

- type
String indicating the type of approximation.

- Q
The matrix

`t(Pl) %*% solve(C,Pl)`

.- type
String indicating the type of approximation.

- Pl.factors
List with elements that can be used to assemble \(P_l\).

- Pr.factors
List with elements that can be used to assemble \(P_r\).

## Details

The approximation is based on a rational approximation of
the fractional operator, resulting in an
approximate model on the form $$P_l u(s) = P_r W,$$
where \(P_j = p_j(L)\) are non-fractional operators defined in terms of
polynomials \(p_j\) for \(j=l,r\). The order of \(p_r\) is given by
`m`

and the order of \(p_l\) is \(m + m_\beta\)
where \(m_\beta\) is the integer part of \(\beta\) if \(\beta>1\) and
\(m_\beta = 1\) otherwise.

The discrete approximation can be written as \(u = P_r x\) where
\(x \sim N(0,Q^{-1})\) and \(Q = P_l^T C^{-1} P_l\).
Note that the matrices \(P_r\) and \(Q\) may be be ill-conditioned
for \(m>1\). In this case, the methods in `operator.operations()`

should be used for operations involving the matrices, since these methods
are more numerically stable.

## Examples

```
# Compute rational approximation of a Gaussian process with a
# Matern covariance function on R
kappa <- 10
sigma <- 1
nu <- 0.8
# create mass and stiffness matrices for a FEM discretization
x <- seq(from = 0, to = 1, length.out = 101)
fem <- rSPDE.fem1d(x)
# compute rational approximation of covariance function at 0.5
tau <- sqrt(gamma(nu) / (sigma^2 * kappa^(2 * nu) *
(4 * pi)^(1 / 2) * gamma(nu + 1 / 2)))
op <- fractional.operators(
L = fem$G + kappa^2 * fem$C, beta = (nu + 1 / 2) / 2,
C = fem$C, scale.factor = kappa^2, tau = tau
)
v <- t(rSPDE.A1d(x, 0.5))
c.approx <- Sigma.mult(op, v)
# plot the result and compare with the true Matern covariance
plot(x, matern.covariance(abs(x - 0.5), kappa, nu, sigma),
type = "l", ylab = "C(h)",
xlab = "h", main = "Matern covariance and rational approximations"
)
lines(x, c.approx, col = 2)
```