Experiment

Go back to the Contents page.

---

Press Show to reveal the code chunks.

---

Let us set some global options for all code chunks in this document.

# Create a clipboard button
source(here::here("clipboard.R")); clipboard

# Set seed for reproducibility
set.seed(1982) 
# Set global options for all code chunks
knitr::opts_chunk$set(
  # Disable messages printed by R code chunks
  message = FALSE,    
  # Disable warnings printed by R code chunks
  warning = FALSE,    
  # Show R code within code chunks in output
  echo = TRUE,        
  # Include both R code and its results in output
  include = TRUE,     
  # Evaluate R code chunks
  eval = TRUE,       
  # Enable caching of R code chunks for faster rendering
  cache = FALSE,      
  # Align figures in the center of the output
  fig.align = "center",
  # Enable retina display for high-resolution figures
  retina = 2,
  # Show errors in the output instead of stopping rendering
  error = TRUE,
  # Do not collapse code and output into a single block
  collapse = FALSE
)
# Start the figure counter
fig_count <- 0
# Define the captioner function
captioner <- function(caption) {
  fig_count <<- fig_count + 1
  paste0("Figure ", fig_count, ": ", caption)
}

# remotes::install_github("davidbolin/rspde", ref = "devel")
# remotes::install_github("davidbolin/metricgraph", ref = "devel")
library(rSPDE)
library(MetricGraph)
library(grateful)

library(ggplot2)
library(reshape2)
library(plotly)

capture.output(
  knitr::purl(here::here("functionality.Rmd"), output = here::here("functionality.R")),
  file = here::here("purl_log.txt")
)
source(here::here("functionality.R"))

1 The equation

We analyze and numerically approximate solutions to fractional diffusion equations on metric graphs of the form \[\begin{equation} \label{eq:maineq} \tag{1} \left\{ \begin{aligned} \partial_t u(s,t) + (\kappa^2 - \Delta_\Gamma)^{\alpha/2} u(s,t) &= f(s,t), && \quad (s,t) \in \Gamma \times (0, T), \\ u(s,0) &= u_0(s), && \quad s \in \Gamma, \end{aligned} \right. \end{equation}\] with \(u(\cdot,t)\) satisfying the Kirchhoff vertex conditions \[\begin{equation} \label{eq:Kcond} \tag{2} \mathcal{K} = \left\{\phi\in C(\Gamma)\;\middle|\; \forall v\in \mathcal{V}:\; \textstyle\sum_{e\in\mathcal{E}_v}\partial_e \phi(v)=0 \right\}. \end{equation}\] Here \(\Gamma = (\mathcal{V},\mathcal{E})\) is a metric graph, \(\kappa>0\), \(\alpha\in(0,2]\) determines the smoothness of \(u(\cdot,t)\), \(\Delta_{\Gamma}\) is the so-called Kirchhoff–Laplacian, and \(f:\Gamma\times (0,T)\to\mathbb{R}\) and \(u_0: \Gamma \to \mathbb{R}\) are fixed functions, called right-hand side and initial condition, respectively.

2 Eigenfunction-based construction of an exact solution

The exact solution to \(\eqref{eq:maineq}\) can be expressed as \[\begin{equation} \label{eq:sol_reprentation} \tag{3} u(s,t) = \displaystyle\sum_{j\in\mathbb{N}}e^{-\lambda^{\alpha/2}_jt}\left(u_0, e_j\right)_{L_2(\Gamma)}e_j(s) + \int_0^t \displaystyle\sum_{j\in\mathbb{N}}e^{-\lambda^{\alpha/2}_j(t-r)}\left(f(\cdot, r), e_j\right)_{L_2(\Gamma)}e_j(s)dr. \end{equation}\]

To construct an exact solution, we expand both the initial condition and the right-hand side in terms of the eigenfunctions. We choose coefficients \(\{x_j\}_{j=1}^{\infty}\) and \(\{y_j\}_{j=1}^{\infty}\) and a scalar function \(g(t)\) and set \[\begin{align} \tag{4} u_0(s) = \sum_{j=1}^{\infty} x_j e_j(s)\quad \text{ and }\quad f(s,t) = g(t) \sum_{j=1}^{\infty} y_j e_j(s). \end{align}\] Substituting these expressions into the solution formula \(\eqref{eq:sol_reprentation}\), we obtain \[\begin{align} \label{sollll} \tag{5} u(s,t) = \sum_{j=1}^{\infty}(x_j+y_j G_j(t))e^{-\lambda^{\alpha/2}_jt}e_j(s),\quad G_j(t)= \int_0^t e^{\lambda^{\alpha/2}_jr}g(r)dr, \end{align}\] where the integral can be evaluated analytically for some choices of \(g(t)\), as shown in the Functionality page. The following list shows how the above expressions are implemented in \(\texttt{R}\). \[\begin{aligned} \text{expression} \iff&\quad \text{In } \texttt{R}\\ \{x_j\}_{j=1}^{\infty} \iff&\quad \texttt{coeff_U_0}\\ \{y_j\}_{j=1}^{\infty} \iff&\quad \texttt{coeff_FF}\\ [e_1 e_2 \dots e_{N_f}] \iff&\quad \texttt{EIGENFUN}\\ u_0(s) \iff&\quad \texttt{U_0 <- EIGENFUN %*% coeff_U_0}\\ f(s,t) \iff&\quad \texttt{FF_true <- EIGENFUN %*%}\\ &\quad\texttt{ outer(1:length(coeff_FF),}\\ &\qquad\qquad\texttt{ 1:length(time_seq),}\\ &\qquad\qquad\texttt{ function(i, j) coeff_FF[i] * g_sin(r = time_seq[j], A = AAA, omega = OMEGA))}\\ u(s,t) \iff&\quad \texttt{U_true <- EIGENFUN %*%}\\ &\quad\texttt{ outer(1:length(coeff_U_0),}\\ &\qquad\qquad\texttt{ 1:length(time_seq),}\\ &\qquad\qquad\texttt{ function(i, j) (coeff_U_0[i] + coeff_FF[i] * G_sin(t = time_seq[j], A = AAA, lambda_j_alpha_half = EIGENVAL_ALPHA[i], omega = OMEGA)) *}\\ &\qquad\qquad\qquad\qquad\qquad\qquad\quad\texttt{exp(-EIGENVAL_ALPHA[i] * time_seq[j]))} \end{aligned}\]

Let us check that \(\eqref{sollll}\) is a solution of \(\eqref{eq:maineq}\). At \(t=0\), we get \[\begin{align*} u(s,0) = \sum_{j=1}^{\infty}(x_j+y_j G_j(0))e^{-\lambda^{\alpha/2}_j0}e_j(s) = \sum_{j=1}^{\infty}x_je_j(s) =u_0(s), \end{align*}\] while for \((s,t) \in \Gamma \times (0, T)\), we obtain \[\begin{align*} \partial_t u(s,t) = \sum_{j=1}^{\infty}y_j G'_j(t)e^{-\lambda^{\alpha/2}_jt}e_j(s) - \sum_{j=1}^{\infty}(x_j+y_j G_j(t))e^{-\lambda^{\alpha/2}_jt}\lambda^{\alpha/2}_je_j(s) \end{align*}\] and \[\begin{align*} (\kappa^2 - \Delta_\Gamma)^{\alpha/2} u(s,t) &= \sum_{j=1}^{\infty}(x_j+y_j G_j(t))e^{-\lambda^{\alpha/2}_jt}(\kappa^2 - \Delta_\Gamma)^{\alpha/2}e_j(s)\\ &= \sum_{j=1}^{\infty}(x_j+y_j G_j(t))e^{-\lambda^{\alpha/2}_jt}\lambda^{\alpha/2}_je_j(s). \end{align*}\] Adding these two and using the fact that \(G'_j(t) = e^{\lambda^{\alpha/2}_jt}g(t)\) yield \[\begin{align*} \partial_t u(s,t) + (\kappa^2 - \Delta_\Gamma)^{\alpha/2} u(s,t)&= g(t)\sum_{j=1}^{\infty}y_je^{\lambda^{\alpha/2}_jt} e^{-\lambda^{\alpha/2}_jt}e_j(s) = f(s,t). \end{align*}\] This shows that \(\eqref{sollll}\) is indeed a solution of \(\eqref{eq:maineq}\).

From the Functionality page, we know that the solution to \(\eqref{eq:maineq}\) can be approximated by a numerical scheme of the form \[\begin{equation} \label{eq:final_scheme2} \tag{6} \mathbf{U}^{k+1} = \mathbf{R}(\mathbf{C}\mathbf{U}^k+\tau \mathbf{F}^{k+1}),\quad \mathbf{R} = \sum_{k=1}^{m+1} a_k\left( \dfrac{\mathbf{L}}{\kappa^2}-p_k\mathbf{C}\right)^{-1} \end{equation}\] Let \(\mathbf{F}\) be the matrix with columns \(\mathbf{F}^k\). This matrix has entries \(\mathbf{F}_{ik}\) given by \[\begin{align} \label{eq:innerprod} \tag{7} (f^{k},\psi^i_h)_{L_2(\Gamma)} = (f(\cdot,t_k),\psi^i_h)_{L_2(\Gamma)} = \sum_{j=1}^{N_f} y_j g(t_k) (e_j,\psi^i_h)_{L_2(\Gamma)}. \end{align}\] To approximate the inner products \((e_j,\psi^i_h)_{L_2(\Gamma)}\), we use the quadrature formula \(\boldsymbol{\psi}_i^\top\mathbf{C}^{\star}\mathbf{e}_j\), where \(\mathbf{e}_j=[e_j(s_1)\dots e_j(s_{N_{h_{\star}}})]^\top\) and \(\boldsymbol{\psi}_i=[\psi^i_h(s_1)\dots \psi^i_h(s_{N_{h_{\star}}})]^\top\) are vectors of function evaluations on a fine spatial mesh with mesh size \(h_{\star}\) and nodes \(\{s_\ell\}_{\ell=1}^{N_{h_{\star}}}\). Matrix \(\mathbf{C}^{\star}\) contains the corresponding quadrature weights, with entries \(\mathbf{C}^{\star}_{nm} = (\psi_{h_\star}^n,\psi_{h_\star}^m)_{L_2(\Gamma)}\) for \(n,m = 1,\dots, N_{h_\star}\). We emphasize that two spatial meshes are involved in this construction. The basis functions \(\{\psi^i_h\}_{i=1}^{N_h}\) are defined on a coarse spatial mesh with mesh size \(h\), while the quadrature is carried out over a finer mesh with associated basis functions \(\{\psi^\ell_{h_\star}\}_{\ell=1}^{N_{h_\star}}\) and mesh size \(h_\star\). If \(\mathbf{N}\) denotes the matrix with entries \(\mathbf{N}_{jk} = y_j g(t_k)\), then \(\mathbf{F}\) can be approximated as \([\boldsymbol{\psi}_1\dots \boldsymbol{\psi}_{N_h}]^\top \mathbf{C}^{\star} [\mathbf{e}_1\dots \mathbf{e}_{N_f}] \mathbf{N}\). In \(\texttt{R}\), \(\mathbf{F}\) and \(\mathbf{U}^{k+1}\) are implemented as \[\begin{aligned} \text{expression} \iff&\quad \text{In } \texttt{R}\\ [\boldsymbol{\psi}_1\dots \boldsymbol{\psi}_{N_h}]^\top \iff&\quad \texttt{t(graph\$fem_basis(overkill_graph\$get_mesh_locations()))}\\ \mathbf{C}^{\star} \iff&\quad \texttt{overkill_graph\$mesh\$C}\\ [\mathbf{e}_1\dots \mathbf{e}_{N_f}] \iff&\quad \texttt{overkill_EIGENFUN}\\ \mathbf{N} \iff&\quad \texttt{outer(1:length(coeff_FF),}\\ &\qquad\quad\;\;\texttt{ 1:length(time_seq),}\\ &\qquad\quad\;\;\texttt{ function(i, j) coeff_FF[i] * g_sin(r = time_seq[j], A = AAA, omega = OMEGA))}\\ \mathbf{F}\iff&\quad \texttt{FF_innerprod <- t(graph\$fem_basis(overkill_graph\$get_mesh_locations())) %*%}\\ &\quad\texttt{overkill_graph\$mesh\$C %*%}\\ &\quad\texttt{overkill_EIGENFUN %*%}\\ &\quad\texttt{outer(1:length(coeff_FF),}\\ &\qquad\qquad\texttt{1:length(time_seq),}\\ &\qquad\qquad\texttt{function(i, j) coeff_FF[i] * g_sin(r = time_seq[j], A = AAA, omega = OMEGA))}\\ \mathbf{U}^{k+1}\iff&\quad \texttt{U_approx[, k + 1] <- as.matrix(my.solver.frac(my_op_frac, my_op_frac\$C %*% U_approx[, k] + time_step * FF_innerprod[, k + 1]))} \end{aligned}\]

T_final <- 2
time_step <- 0.05
h <- 0.05
kappa <- 4
alpha <- 1.8 
m = 4
beta <- alpha/2
N_finite = 4 # choose even
adjusted_N_finite <- N_finite + N_finite/2 + 1
# Coefficients for u_0 and f
coeff_U_0 <- 50*(1:adjusted_N_finite)^-1
coeff_U_0[-5] <- 0
coeff_FF <- rep(0, adjusted_N_finite)
coeff_FF[7] <- 10


time_seq <- seq(0, T_final, by = time_step)
graph <- gets.graph.tadpole(h = h)
graph$compute_fem()
G <- graph$mesh$G
C <- graph$mesh$C
L <- kappa^2*C + G
I <- Matrix::Diagonal(nrow(C))

eigen_params <- gets.eigen.params(N_finite = N_finite, kappa = kappa, alpha = alpha, graph = graph)
EIGENVAL_ALPHA <- eigen_params$EIGENVAL_ALPHA
EIGENFUN <- eigen_params$EIGENFUN

AAA = 1
OMEGA = pi

U_0 <- EIGENFUN %*% coeff_U_0

U_true <- EIGENFUN %*% 
  outer(1:length(coeff_U_0), 
        1:length(time_seq), 
        function(i, j) (coeff_U_0[i] + coeff_FF[i] * G_sin(t = time_seq[j], A = AAA, lambda_j_alpha_half = EIGENVAL_ALPHA[i], omega = OMEGA) ) * exp(-EIGENVAL_ALPHA[i] * time_seq[j]))

overkill_graph <- gets.graph.tadpole(h = 0.0001)
overkill_graph$compute_fem()
overkill_EIGENFUN <- gets.eigen.params(N_finite = N_finite, kappa = kappa, alpha = alpha, graph = overkill_graph)$EIGENFUN


FF_innerprod <- t(graph$fem_basis(overkill_graph$get_mesh_locations())) %*%
  overkill_graph$mesh$C %*%
  overkill_EIGENFUN %*%
  outer(1:length(coeff_FF), 
        1:length(time_seq), 
        function(i, j) coeff_FF[i] * g_sin(r = time_seq[j], A = AAA, omega = OMEGA))

FF_true <- EIGENFUN %*% outer(1:length(coeff_FF), 
                         1:length(time_seq), 
                         function(i, j) coeff_FF[i] * g_sin(r = time_seq[j], A = AAA, omega = OMEGA))

# Numerical solution
my_op_frac <- my.fractional.operators.frac(L, beta, C, scale.factor = kappa^2, m = m, time_step)


U_approx <- solve_fractional_evolution(my_op_frac, time_step, time_seq, val_at_0 = U_0, RHST = FF_innerprod)
U_approx2 <- U_approx
XX <- U_true - U_approx

error <- sqrt(as.double(time_step * sum(XX * (C %*% XX))))

The following is the error reported in the paper for the considered experiment.

print(paste0("Total error = ", formatC(error, format = "f", digits = 9)))

## [1] "Total error = 0.469623220"

Because of GitHub storage limits, we project the solution (and all other visualizations) to a coarser mesh as follows.

aux_graph <- gets.graph.tadpole(h = 0.01)
A_proj <- graph$fem_basis(aux_graph$get_mesh_locations())
FF_true <- A_proj %*% FF_true
U_true <- A_proj %*% U_true
U_approx <- A_proj %*% U_approx

start_idx <- which.min(abs(time_seq - 0.5))
end_idx <- which.min(abs(time_seq - 1.5))
idx <- seq(start_idx, end_idx, by = 4)

idx_for1 <- which.min(abs(time_seq - 0.5))

exp1_fig <- graph.plotter.3d.single(aux_graph, FF_true[, idx], time_seq[idx])
exp1_fig

Figure 1: Time evolution of the right-hand side function \(f\).

save(exp1_fig, file = here::here("data_files/exp1_fig.RData"))

diff <- U_true - U_approx
abs_error <- abs(diff)[, idx]
exp1_error <- graph.plotter.3d.single(aux_graph, abs_error, time_seq[idx])
exp1_error

Figure 2: Time evolution of the absolute difference between the exact and approximate solution.

save(exp1_error, file = here::here("data_files/exp1_error.RData"))

error_tadpole <- graph.plotter.3d.onecol(graph, abs(XX[,idx_for1]))
save(error_tadpole, file = here::here("data_files/exp1_error_tadpole.RData"))
error_tadpole

Figure 3: Error at some time point.

f_tadpole <- graph.plotter.3d.onecol(graph, U_approx2[,idx_for1])
save(f_tadpole, file = here::here("data_files/exp1_f_tadpole.RData"))
f_tadpole

Figure 4: Error at some time point.

idx2 <- seq(start_idx, end_idx, by = 10)
exp1_com <- graph.plotter.3d.comparer(aux_graph, U_true[,idx2], U_approx[,idx2], time_seq[idx2])
exp1_com

Figure 5: Comparison of the exact and approximate solution at selected time points.

save(exp1_com, file = here::here("data_files/exp1_com.RData"))

3 References

grateful::cite_packages(output = "paragraph", out.dir = ".")

We used R version 4.5.2 (R Core Team 2025) and the following R packages: akima v. 0.6.3.6 (Akima and Gebhardt 2025), expm v. 1.0.0 (Maechler, Dutang, and Goulet 2024), fmesher v. 0.5.0 (Lindgren 2025), gsignal v. 0.3.7 (Van Boxtel, G.J.M., et al. 2021), here v. 1.0.1 (Müller 2020), htmltools v. 0.5.8.1 (Cheng et al. 2024), INLA v. 25.11.22 (Rue, Martino, and Chopin 2009; Lindgren, Rue, and Lindström 2011; Martins et al. 2013; Lindgren and Rue 2015; De Coninck et al. 2016; Rue et al. 2017; Verbosio et al. 2017; Bakka et al. 2018; Kourounis, Fuchs, and Schenk 2018), inlabru v. 2.13.0 (Yuan et al. 2017; Bachl et al. 2019), knitr v. 1.50 (Xie 2014, 2015, 2025a), Matrix v. 1.7.3 (Bates, Maechler, and Jagan 2025), MetricGraph v. 1.5.0.9000 (Bolin, Simas, and Wallin 2023a, 2023b, 2024, 2025; Bolin et al. 2024), neuralnet v. 1.44.2 (Fritsch, Guenther, and Wright 2019), orthopolynom v. 1.0.6.1 (Novomestky 2022), patchwork v. 1.3.1 (Pedersen 2025), pbmcapply v. 1.5.1 (Kuang, Kong, and Napolitano 2022), plotly v. 4.11.0 (Sievert 2020), posterdown v. 1.0 (Thorne 2019), pracma v. 2.4.4 (Borchers 2023), qrcode v. 0.3.0 (Onkelinx and Teh 2024), RColorBrewer v. 1.1.3 (Neuwirth 2022), RefManageR v. 1.4.0 (McLean 2014, 2017), renv v. 1.1.5 (Ushey and Wickham 2025), reshape2 v. 1.4.4 (Wickham 2007), reticulate v. 1.44.1 (Ushey, Allaire, and Tang 2025), rmarkdown v. 2.30 (Xie, Allaire, and Grolemund 2018; Xie, Dervieux, and Riederer 2020; Allaire et al. 2025), rSPDE v. 2.5.1.9000 (Bolin and Kirchner 2020; Bolin and Simas 2023; Bolin, Simas, and Xiong 2024), RSpectra v. 0.16.2 (Qiu and Mei 2024), scales v. 1.4.0 (Wickham, Pedersen, and Seidel 2025), slackr v. 3.4.0 (Kaye et al. 2025), tidyverse v. 2.0.0 (Wickham et al. 2019), viridisLite v. 0.4.2 (Garnier et al. 2023), xaringan v. 0.31 (Xie 2025b), xaringanExtra v. 0.8.0 (Aden-Buie and Warkentin 2024), xaringanthemer v. 0.4.4 (Aden-Buie 2025).

Aden-Buie, Garrick. 2025. xaringanthemer: Custom “xaringan” CSS Themes. https://doi.org/10.32614/CRAN.package.xaringanthemer.

Aden-Buie, Garrick, and Matthew T. Warkentin. 2024. xaringanExtra: Extras and Extensions for “xaringan” Slides. https://doi.org/10.32614/CRAN.package.xaringanExtra.

Akima, Hiroshi, and Albrecht Gebhardt. 2025. akima: Interpolation of Irregularly and Regularly Spaced Data. https://doi.org/10.32614/CRAN.package.akima.

Allaire, JJ, Yihui Xie, Christophe Dervieux, Jonathan McPherson, Javier Luraschi, Kevin Ushey, Aron Atkins, et al. 2025. rmarkdown: Dynamic Documents for r. https://github.com/rstudio/rmarkdown.

Bachl, Fabian E., Finn Lindgren, David L. Borchers, and Janine B. Illian. 2019. “inlabru: An R Package for Bayesian Spatial Modelling from Ecological Survey Data.” Methods in Ecology and Evolution 10: 760–66. https://doi.org/10.1111/2041-210X.13168.

Bakka, Haakon, Håvard Rue, Geir-Arne Fuglstad, Andrea I. Riebler, David Bolin, Janine Illian, Elias Krainski, Daniel P. Simpson, and Finn K. Lindgren. 2018. “Spatial Modelling with INLA: A Review.” WIRES (Invited Extended Review) xx (Feb): xx–. http://arxiv.org/abs/1802.06350.

Bates, Douglas, Martin Maechler, and Mikael Jagan. 2025. Matrix: Sparse and Dense Matrix Classes and Methods. https://doi.org/10.32614/CRAN.package.Matrix.

Bolin, David, and Kristin Kirchner. 2020. “The Rational SPDE Approach for Gaussian Random Fields with General Smoothness.” Journal of Computational and Graphical Statistics 29 (2): 274–85. https://doi.org/10.1080/10618600.2019.1665537.

Bolin, David, Mihály Kovács, Vivek Kumar, and Alexandre B. Simas. 2024. “Regularity and Numerical Approximation of Fractional Elliptic Differential Equations on Compact Metric Graphs.” Mathematics of Computation 93 (349): 2439–72. https://doi.org/10.1090/mcom/3929.

Bolin, David, and Alexandre B. Simas. 2023. rSPDE: Rational Approximations of Fractional Stochastic Partial Differential Equations. https://CRAN.R-project.org/package=rSPDE.

Bolin, David, Alexandre B. Simas, and Jonas Wallin. 2023a. MetricGraph: Random Fields on Metric Graphs. https://CRAN.R-project.org/package=MetricGraph.

———. 2023b. “Statistical Inference for Gaussian Whittle-Matérn Fields on Metric Graphs.” arXiv Preprint arXiv:2304.10372. https://doi.org/10.48550/arXiv.2304.10372.

———. 2024. “Gaussian Whittle-Matérn Fields on Metric Graphs.” Bernoulli 30 (2): 1611–39. https://doi.org/10.3150/23-BEJ1647.

———. 2025. “Markov Properties of Gaussian Random Fields on Compact Metric Graphs.” Bernoulli. https://doi.org/10.48550/arXiv.2304.03190.

Bolin, David, Alexandre B. Simas, and Zhen Xiong. 2024. “Covariance-Based Rational Approximations of Fractional SPDEs for Computationally Efficient Bayesian Inference.” Journal of Computational and Graphical Statistics 33 (1): 64–74. https://doi.org/10.1080/10618600.2023.2231051.

Borchers, Hans W. 2023. pracma: Practical Numerical Math Functions. https://doi.org/10.32614/CRAN.package.pracma.

Cheng, Joe, Carson Sievert, Barret Schloerke, Winston Chang, Yihui Xie, and Jeff Allen. 2024. htmltools: Tools for HTML. https://github.com/rstudio/htmltools.

De Coninck, Arne, Bernard De Baets, Drosos Kourounis, Fabio Verbosio, Olaf Schenk, Steven Maenhout, and Jan Fostier. 2016. “Needles: Toward Large-Scale Genomic Prediction with Marker-by-Environment Interaction.” Genetics 203 (1): 543–55. https://doi.org/10.1534/genetics.115.179887.

Fritsch, Stefan, Frauke Guenther, and Marvin N. Wright. 2019. neuralnet: Training of Neural Networks. https://doi.org/10.32614/CRAN.package.neuralnet.

Garnier, Simon, Ross, Noam, Rudis, Robert, Camargo, et al. 2023. viridis(Lite) - Colorblind-Friendly Color Maps for r. https://doi.org/10.5281/zenodo.4678327.

Kaye, Matt, Bob Rudis, Andrie de Vries, and Jonathan Sidi. 2025. slackr: Send Messages, Images, r Objects and Files to “Slack” Channels/Users. https://github.com/mrkaye97/slackr.

Kourounis, D., A. Fuchs, and O. Schenk. 2018. “Towards the Next Generation of Multiperiod Optimal Power Flow Solvers.” IEEE Transactions on Power Systems PP (99): 1–10. https://doi.org/10.1109/TPWRS.2017.2789187.

Kuang, Kevin, Quyu Kong, and Francesco Napolitano. 2022. pbmcapply: Tracking the Progress of Mc*pply with Progress Bar. https://doi.org/10.32614/CRAN.package.pbmcapply.

Lindgren, Finn. 2025. fmesher: Triangle Meshes and Related Geometry Tools. https://github.com/inlabru-org/fmesher.

Lindgren, Finn, and Håvard Rue. 2015. “Bayesian Spatial Modelling with R-INLA.” Journal of Statistical Software 63 (19): 1–25. http://www.jstatsoft.org/v63/i19/.

Lindgren, Finn, Håvard Rue, and Johan Lindström. 2011. “An Explicit Link Between Gaussian Fields and Gaussian Markov Random Fields: The Stochastic Partial Differential Equation Approach (with Discussion).” Journal of the Royal Statistical Society B 73 (4): 423–98.

Maechler, Martin, Christophe Dutang, and Vincent Goulet. 2024. expm: Matrix Exponential, Log, “etc”. https://doi.org/10.32614/CRAN.package.expm.

Martins, Thiago G., Daniel Simpson, Finn Lindgren, and Håvard Rue. 2013. “Bayesian Computing with INLA: New Features.” Computational Statistics and Data Analysis 67: 68–83.

McLean, Mathew William. 2014. Straightforward Bibliography Management in r Using the RefManager Package. https://arxiv.org/abs/1403.2036.

———. 2017. “RefManageR: Import and Manage BibTeX and BibLaTeX References in r.” The Journal of Open Source Software. https://doi.org/10.21105/joss.00338.

Müller, Kirill. 2020. here: A Simpler Way to Find Your Files. https://doi.org/10.32614/CRAN.package.here.

Neuwirth, Erich. 2022. RColorBrewer: ColorBrewer Palettes.

Novomestky, Frederick. 2022. orthopolynom: Collection of Functions for Orthogonal and Orthonormal Polynomials. https://doi.org/10.32614/CRAN.package.orthopolynom.

Onkelinx, Thierry, and Victor Teh. 2024. qrcode: Generate QRcodes with r. Version 0.3.0. https://doi.org/10.5281/zenodo.5040088.

Pedersen, Thomas Lin. 2025. patchwork: The Composer of Plots. https://doi.org/10.32614/CRAN.package.patchwork.

Qiu, Yixuan, and Jiali Mei. 2024. RSpectra: Solvers for Large-Scale Eigenvalue and SVD Problems. https://doi.org/10.32614/CRAN.package.RSpectra.

R Core Team. 2025. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. https://www.R-project.org/.

Rue, Håvard, Sara Martino, and Nicholas Chopin. 2009. “Approximate Bayesian Inference for Latent Gaussian Models Using Integrated Nested Laplace Approximations (with Discussion).” Journal of the Royal Statistical Society B 71: 319–92.

Rue, Håvard, Andrea I. Riebler, Sigrunn H. Sørbye, Janine B. Illian, Daniel P. Simpson, and Finn K. Lindgren. 2017. “Bayesian Computing with INLA: A Review.” Annual Reviews of Statistics and Its Applications 4 (March): 395–421. http://arxiv.org/abs/1604.00860.

Sievert, Carson. 2020. Interactive Web-Based Data Visualization with r, Plotly, and Shiny. Chapman; Hall/CRC. https://plotly-r.com.

Thorne, W. Brent. 2019. posterdown: An r Package Built to Generate Reproducible Conference Posters for the Academic and Professional World Where Powerpoint and Pages Just Won’t Cut It. https://github.com/brentthorne/posterdown.

Ushey, Kevin, JJ Allaire, and Yuan Tang. 2025. reticulate: Interface to “Python”. https://doi.org/10.32614/CRAN.package.reticulate.

Ushey, Kevin, and Hadley Wickham. 2025. renv: Project Environments. https://rstudio.github.io/renv/.

Van Boxtel, G.J.M., et al. 2021. gsignal: Signal Processing. https://github.com/gjmvanboxtel/gsignal.

Verbosio, Fabio, Arne De Coninck, Drosos Kourounis, and Olaf Schenk. 2017. “Enhancing the Scalability of Selected Inversion Factorization Algorithms in Genomic Prediction.” Journal of Computational Science 22 (Supplement C): 99–108. https://doi.org/10.1016/j.jocs.2017.08.013.

Wickham, Hadley. 2007. “Reshaping Data with the reshape Package.” Journal of Statistical Software 21 (12): 1–20. http://www.jstatsoft.org/v21/i12/.

Wickham, Hadley, Mara Averick, Jennifer Bryan, Winston Chang, Lucy D’Agostino McGowan, Romain François, Garrett Grolemund, et al. 2019. “Welcome to the tidyverse.” Journal of Open Source Software 4 (43): 1686. https://doi.org/10.21105/joss.01686.

Wickham, Hadley, Thomas Lin Pedersen, and Dana Seidel. 2025. scales: Scale Functions for Visualization. https://scales.r-lib.org.

Xie, Yihui. 2014. “knitr: A Comprehensive Tool for Reproducible Research in R.” In Implementing Reproducible Computational Research, edited by Victoria Stodden, Friedrich Leisch, and Roger D. Peng. Chapman; Hall/CRC.

———. 2015. Dynamic Documents with R and Knitr. 2nd ed. Boca Raton, Florida: Chapman; Hall/CRC. https://yihui.org/knitr/.

———. 2025a. knitr: A General-Purpose Package for Dynamic Report Generation in R. https://yihui.org/knitr/.

———. 2025b. xaringan: Presentation Ninja. https://doi.org/10.32614/CRAN.package.xaringan.

Xie, Yihui, J. J. Allaire, and Garrett Grolemund. 2018. R Markdown: The Definitive Guide. Boca Raton, Florida: Chapman; Hall/CRC. https://bookdown.org/yihui/rmarkdown.

Xie, Yihui, Christophe Dervieux, and Emily Riederer. 2020. R Markdown Cookbook. Boca Raton, Florida: Chapman; Hall/CRC. https://bookdown.org/yihui/rmarkdown-cookbook.

Yuan, Yuan, Bachl, Fabian E., Lindgren, Finn, Borchers, et al. 2017. “Point Process Models for Spatio-Temporal Distance Sampling Data from a Large-Scale Survey of Blue Whales.” Ann. Appl. Stat. 11 (4): 2270–97. https://doi.org/10.1214/17-AOAS1078.