Go back to the Contents page.
Press Show to reveal the code chunks.
# 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)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}:\; \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.
Let \(\alpha\in(0,2]\) and \(U_h^\tau\) denote the sequence of approximations of the solution to the weak form of problem \(\eqref{eq:maineq}\) at each time step on a mesh indexed by \(h\). The scheme computes \(U_h^{\tau}\subset V_h\), which solves \[\begin{equation} \left\{ \begin{aligned} \langle\delta U_h^{k+1},\phi\rangle + \mathfrak{a}( U_h^{k+1},\phi) &= \langle f^{k+1},\phi\rangle, \quad\forall\phi\in V_h,\quad k=0,\dots, N-1, \\ U^0_h &= P_hu_0, \end{aligned} \right. \end{equation}\] where \(f^{k+1} = \displaystyle\dfrac{1}{\tau}\int_{t_k}^{t_{k+1}}f(t)dt\). The above expression can be equivalently written as \[\begin{equation} \left\{ \begin{aligned} \langle\dfrac{U_h^{k+1} - U_h^{k}}{\tau},\phi\rangle + \mathfrak{a}( U_h^{k+1},\phi) &= \langle f^{k+1},\phi\rangle, \quad\forall\phi\in V_h,\quad k=0,\dots, N-1, \\ U^0_h &= P_hu_0, \end{aligned} \right. \end{equation}\] or \[\begin{equation} \label{system:fully_discrete_scheme} \tag{3} \left\{ \begin{aligned} \langle U_h^{k+1},\phi\rangle + \tau\mathfrak{a}( U_h^{k+1},\phi) &= \langle U_h^{k},\phi\rangle + \tau\langle f^{k+1},\phi\rangle, \quad\forall\phi\in V_h,\quad k=0,\dots, N-1, \\ U^0_h &= P_hu_0. \end{aligned} \right. \end{equation}\]
At each time step \(t_k\), the finite element solution \(U_h^k\in V_h\) to \(\eqref{system:fully_discrete_scheme}\) can be expressed as a linear combination of the basis functions \(\{\psi^i_h\}_{i=1}^{N_h}\) introduced in the Preliminaries page, namely, \[\begin{align} \label{num_sol} \tag{4} U_h^k(s) = \sum_{i=1}^{N_h}u_i^k\psi^i_h(s), \;s\in\Gamma. \end{align}\]
Replacing \(\eqref{num_sol}\) into \(\eqref{system:fully_discrete_scheme}\) yields the following system \[\begin{align*} \sum_{j=1}^{N_h}u_j^{k+1}[(\psi_h^j,\psi_h^i)_{L_2(\Gamma)}+ \tau\mathfrak{a}(\psi_h^j,\psi_h^i)] = \sum_{j=1}^{N_h}u_j^{k}(\psi_h^j,\psi_h^i)_{L_2(\Gamma)}+\tau( f^{k+1},\psi_h^i)_{L_2(\Gamma)},\quad i = 1,\dots, N_h. \end{align*}\] In matrix notation, \[\begin{align} \label{diff_eq_discrete} (\mathbf{C}+\tau \mathbf{L}^{\alpha/2})\mathbf{U}^{k+1} = \mathbf{C}\mathbf{U}^k+\tau \mathbf{F}^{k+1}, \end{align}\] or by introducing the scaling parameter \(\kappa^2>0\), \[\begin{align} (\mathbf{C}+\tau (\kappa^2)^{\alpha/2}(\mathbf{L}/\kappa^2)^{\alpha/2})\mathbf{U}^{k+1} = \mathbf{C}\mathbf{U}^k+\tau \mathbf{F}^{k+1}, \end{align}\] where \(\mathbf{C}\) has entries \(\mathbf{C}_{i,j} = (\psi_h^i,\psi_h^j)_{L_2(\Gamma)}\), \(\mathbf{L}^{\alpha/2}\) has entries \(\mathfrak{a}(\psi_h^i,\psi_h^j)\), \(\mathbf{U}^k\) has components \(u_j^k\), and \(\mathbf{F}^k\) has components \(( f^{k},\psi_h^i)_{L_2(\Gamma)}\). Applying \((\mathbf{L}/\kappa^2)^{-\alpha/2}\) to both sides yields \[\begin{equation} ((\mathbf{L}/\kappa^2)^{-\alpha/2}\mathbf{C}+\tau (\kappa^2)^{\alpha/2}\mathbf{I})\mathbf{U}^{k+1} = (\mathbf{L}/\kappa^2)^{-\alpha/2}(\mathbf{C}\mathbf{U}^k+\tau \mathbf{F}^{k+1}). \end{equation}\] Following Bolin and Kirchner (2020), we approximate \((\mathbf{L}/\kappa^2)^{-\alpha/2}\) by \(\mathbf{P}_\ell^{-\top}\mathbf{P}_r^\top\) to arrive at \[\begin{equation} \label{eq:scheme2} \tag{5} (\mathbf{P}_\ell^{-\top}\mathbf{P}_r^\top \mathbf{C}+\tau(\kappa^2)^{\alpha/2} \mathbf{I})\mathbf{U}^{k+1} = \mathbf{P}_\ell^{-\top}\mathbf{P}_r^\top(\mathbf{C}\mathbf{U}^k+\tau \mathbf{F}^{k+1}). \end{equation}\] where \[\begin{equation} \label{eq:PLPR} \tag{6} \mathbf{P}_r = \prod_{i=1}^m \left(\mathbf{I}-r_{1i}\dfrac{\mathbf{C}^{-1}\mathbf{L}}{\kappa^2}\right)\quad\text{and}\quad \mathbf{P}_\ell = \dfrac{1}{\texttt{factor}}\mathbf{C}\prod_{j=1}^{m+1} \left(\mathbf{I}-r_{2j}\dfrac{\mathbf{C}^{-1}\mathbf{L}}{\kappa^2}\right), \end{equation}\] and \(\texttt{factor} = \dfrac{c_m}{b_{m+1}}\), and \(\{r_{1i}\}_{i=1}^m\) and \(\{r_{2j}\}_{j=1}^{m+1}\) are the roots of \(q_1(x) =\sum_{i=0}^mc_ix^{i}\) and \(q_2(x)=\sum_{j=0}^{m+1}b_jx^{j}\), respectively. The coefficients \(\{c_i\}_{i=0}^m\) and \(\{b_j\}_{j=0}^{m+1}\) are determined as the best rational approximation \(q_1/q_2\) of the function \(x^{\alpha/2-1}\) over the interval \(J_h: = [\kappa^{2}\lambda_{N_h,h}^{-1}, \kappa^{2}\lambda_{1,h}^{-1}]\), where \(\lambda_{1,h}, \lambda_{N_h,h}>0\) are the smallest and the largest eigenvalue of \(L_h\), respectively.
For the sake of clarity, we note that the numerical implementation of
Bolin and Kirchner (2020) actually defines
\(\mathbf{P}_r\) and \(\mathbf{P}_\ell\) as \[\begin{equation}
\label{eq:PLPRbolin}
\tag{7}
\mathbf{P}_r = \prod_{i=1}^m
\left(\mathbf{I}-r_{1i}\dfrac{\mathbf{C}^{-1}\mathbf{L}}{\kappa^2}\right)\quad\text{and}\quad
\mathbf{P}_\ell =
\dfrac{\kappa^{2\beta}}{\texttt{factor}}\mathbf{C}\prod_{j=1}^{m+1}
\left(\mathbf{I}-r_{2j}\dfrac{\mathbf{C}^{-1}\mathbf{L}}{\kappa^2}\right),
\end{equation}\] where \(\beta =
\alpha/2\) and the scaling factor \((\kappa^2)^{\alpha/2}\) or \(\kappa^{2\beta}\) is already incorporated
in \(\mathbf{P}_\ell\), a convention we
adopt in the following. With this under consideration, we can rewrite
\(\eqref{eq:scheme2}\) as \[\begin{equation}
\tag{8}
(\mathbf{P}_r^\top \mathbf{C}+\tau \mathbf{P}_\ell^\top)\mathbf{U}^{k+1}
= \mathbf{P}_r^\top(\mathbf{C}\mathbf{U}^k+\tau \mathbf{F}^{k+1}),
\label{eq:scheme}
\end{equation}\] where
\[\begin{equation}
\mathbf{P}_r^\top = \prod_{i=1}^m
\left(\mathbf{I}-r_{1i}\dfrac{\mathbf{L}\mathbf{C}^{-1}}{\kappa^2}\right)\quad\text{and}\quad
\mathbf{P}_\ell^\top =
\dfrac{\kappa^{2\beta}}{\texttt{factor}}\prod_{j=1}^{m+1}
\left(\mathbf{I}-r_{2j}\dfrac{\mathbf{L}\mathbf{C}^{-1}}{\kappa^2}\right)\cdot
\mathbf{C}
\end{equation}\] since \(\mathbf{L}\) and \(\mathbf{C}^{-1}\) are symmetric and the
factors in the product commute. Replacing these two into \(\eqref{eq:scheme}\) yields \[\begin{equation}
\left(\prod_{i=1}^m
\left(\mathbf{I}-r_{1i}\dfrac{\mathbf{L}\mathbf{C}^{-1}}{\kappa^2}\right)+\dfrac{\tau
\kappa^{2\beta}}{\texttt{factor}}\prod_{j=1}^{m+1}
\left(\mathbf{I}-r_{2j}\dfrac{\mathbf{L}\mathbf{C}^{-1}}{\kappa^2}\right)\right)\mathbf{C}\mathbf{U}^{k+1}
= \prod_{i=1}^m
\left(\mathbf{I}-r_{1i}\dfrac{\mathbf{L}\mathbf{C}^{-1}}{\kappa^2}\right)\cdot
(\mathbf{C}\mathbf{U}^k+\tau \mathbf{F}^{k+1}),
\end{equation}\] that is, \[\begin{equation}
\label{eq:final_scheme}
\tag{9}
\mathbf{U}^{k+1} = \mathbf{C}^{-1}\left(\prod_{i=1}^m
\left(\mathbf{I}-r_{1i}\dfrac{\mathbf{L}\mathbf{C}^{-1}}{\kappa^2}\right)+\dfrac{\tau
\kappa^{2\beta}}{\texttt{factor}}\prod_{j=1}^{m+1}
\left(\mathbf{I}-r_{2j}\dfrac{\mathbf{L}\mathbf{C}^{-1}}{\kappa^2}\right)\right)^{-1}
\prod_{i=1}^m
\left(\mathbf{I}-r_{1i}\dfrac{\mathbf{L}\mathbf{C}^{-1}}{\kappa^2}\right)\cdot
(\mathbf{C}\mathbf{U}^k+\tau \mathbf{F}^{k+1}).
\end{equation}\] Considering the partial fraction decomposition
\[\begin{equation}
\label{eq:partial_fraction}
\tag{10}
\dfrac{\prod_{i=1}^m (1-r_{1i}x)}{\prod_{i=1}^m (1-r_{1i}x)+\dfrac{\tau
\kappa^{2\beta}}{\texttt{factor}} \prod_{j=1}^{m+1}
(1-r_{2j}x)}=\sum_{k=1}^{m+1} a_k(x-p_k)^{-1} + r,
\end{equation}\] scheme \(\eqref{eq:final_scheme}\) can be expressed
as \[\begin{equation}
\label{eq:final_scheme2}
\tag{11}
\mathbf{U}^{k+1} = \mathbf{C}^{-1}\left(\sum_{k=1}^{m+1} a_k\left(
\dfrac{\mathbf{L}\mathbf{C}^{-1}}{\kappa^2}-p_k\mathbf{I}\right)^{-1} +
r\mathbf{I}\right) (\mathbf{C}\mathbf{U}^k+\tau \mathbf{F}^{k+1})
\end{equation}\]
In practice, since the rational function in \(\eqref{eq:partial_fraction}\) is proper, there is no remainder \(r\). Moreover, since \(\left( \dfrac{\mathbf{L}\mathbf{C}^{-1}}{\kappa^2}-p_k\mathbf{I}\right)^{-1} = \mathbf{C}\left( \dfrac{\mathbf{L}}{\kappa^2}-p_k\mathbf{C}\right)^{-1}\), we have that \(\eqref{eq:final_scheme2}\) can be rewritten as
\[\begin{equation} \label{eq:final_scheme3} \tag{12} \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}\]
my.get.roots()For each rational order \(m\)
(1,2,3,4) and smoothness parameter \(\beta\) (= \(\alpha/2\) with \(\alpha\) between 0.5 and 2), function
my.get.roots() (adapted from the rSPDE
package) returns \(\texttt{factor} =
\dfrac{c_m}{b_{m+1}}\), and the roots \(\{r_{1i}\}_{i=1}^m\) and \(\{r_{2j}\}_{j=1}^{m+1}\).
# Function to compute the roots and factor for the rational approximation
my.get.roots <- function(m, # rational order, m = 1, 2, 3, or 4
beta # smoothness parameter, beta = alpha/2 with alpha between 0.5 and 2
) {
m1table <- rSPDE:::m1table
m2table <- rSPDE:::m2table
m3table <- rSPDE:::m3table
m4table <- rSPDE:::m4table
mt <- get(paste0("m", m, "table"))
rb <- rep(0, m + 1)
rc <- rep(0, m)
if(m == 1) {
rc = approx(mt$beta, mt[[paste0("rc")]], beta)$y
} else {
rc = sapply(1:m, function(i) {
approx(mt$beta, mt[[paste0("rc.", i)]], beta)$y
})
}
rb = sapply(1:(m+1), function(i) {
approx(mt$beta, mt[[paste0("rb.", i)]], xout = beta)$y
})
factor = approx(mt$beta, mt$factor, xout = beta)$y
return(list(pl_roots = rb, # roots \{r_{2j}\}_{j=1}^{m+1}
pr_roots = rc, # roots \{r_{1i}\}_{i=1}^m
factor = factor # this is c_m/b_{m+1}
))
}poly.from.roots()Function poly.from.roots() computes the coefficients of
a polynomial from its roots.
compute.partial.fraction.param()Given factor\(=\texttt{factor}
= \dfrac{c_m}{b_{m+1}}\), pr_roots\(=\{r_{1i}\}_{i=1}^m\),
pl_roots\(=\{r_{2j}\}_{j=1}^{m+1}\),
time_step\(=\tau\), and
scaling\(=\kappa^{2\beta}\), function
compute.partial.fraction.param() computes the parameters
for the partial fraction decomposition \(\eqref{eq:partial_fraction}\).
# Function to compute the parameters for the partial fraction decomposition
compute.partial.fraction.param <- function(factor, # c_m/b_{m+1}
pr_roots, # roots \{r_{1i}\}_{i=1}^m
pl_roots, # roots \{r_{2j}\}_{j=1}^{m+1}
time_step, # \tau
scaling # \kappa^{2\beta}
) {
pr_coef <- c(0, poly.from.roots(pr_roots))
pl_coef <- poly.from.roots(pl_roots)
factor_pr_coef <- pr_coef
pr_plus_pl_coef <- factor_pr_coef + ((scaling * time_step)/factor) * pl_coef
res <- gsignal::residue(factor_pr_coef, pr_plus_pl_coef)
return(list(r = res$r, # residues \{a_k\}_{k=1}^{m+1}
p = res$p, # poles \{p_k\}_{k=1}^{m+1}
k = res$k # remainder r
))
}my.fractional.operators.frac()Given the Laplacian matrix L, the smoothness parameter
beta, the mass matrix C (not lumped), the
scaling factor scale.factor\(=\kappa^2\), the rational order
m, and the time step time_step\(=\tau\), function
my.fractional.operators.frac() computes the fractional
operator and returns a list containing the necessary matrices and
parameters for the fractional diffusion equation.
# Function to compute the fractional operator
my.fractional.operators.frac <- function(L, # Laplacian matrix
beta, # smoothness parameter beta
C, # mass matrix (not lumped)
scale.factor, # scaling parameter = kappa^2
m = 1, # rational order, m = 1, 2, 3, or 4
time_step # time step = tau
) {
I <- Matrix::Diagonal(dim(C)[1])
L <- L / scale.factor
if(beta == 1){
L <- L * scale.factor^beta
return(list(C = C, # mass matrix
L = L, # Laplacian matrix scaled
m = m, # rational order
beta = beta, # smoothness parameter
LHS = C + time_step * L # left-hand side of the linear system
))
} else {
scaling <- scale.factor^beta
roots <- my.get.roots(m, beta)
poles_rs_k <- compute.partial.fraction.param(roots$factor, roots$pr_roots, roots$pl_roots, time_step, scaling)
partial_fraction_terms <- list()
for (i in 1:(m+1)) {
# Here is where the terms in the sum in eq 12 are computed
partial_fraction_terms[[i]] <- (L - poles_rs_k$p[i] * C)#/poles_rs_k$r[i]
}
return(list(C = C, # mass matrix
L = L, # Laplacian matrix scaled
m = m, # rational order
beta = beta, # smoothness parameter
partial_fraction_terms = partial_fraction_terms, # partial fraction terms
residues = poles_rs_k$r # residues \{a_k\}_{k=1}^{m+1}
))
}
}my.solver.frac()Given the object returned by
my.fractional.operators.frac() and a vector v,
function my.solver.frac() solves the linear \(\eqref{eq:final_scheme2}\) for the vector
v. If beta = 1, it solves the system directly;
otherwise, it uses the partial fraction decomposition.
# Function to solve the iteration
my.solver.frac <- function(obj, # object returned by my.fractional.operators.frac()
v # vector to be solved for
){
beta <- obj$beta
m <- obj$m
if (beta == 1){
return(solve(obj$LHS, v) # solve the linear system directly for beta = 1
)
} else {
partial_fraction_terms <- obj$partial_fraction_terms
residues <- obj$residues
output <- v*0
for (i in 1:(m+1)) {output <- output + residues[i] * solve(partial_fraction_terms[[i]], v)}
return(output # solve the linear system using the partial fraction decomposition
)
}
}
# my.solver.frac <- function(obj, v) {
# beta <- obj$beta
# m <- obj$m
#
# if (beta == 1) {
# return(solve(obj$LHS, v))
# }
#
# partial_fraction_terms <- obj$partial_fraction_terms
# residues <- obj$residues
#
# # --- PARALLEL PART ---
# results <- parallel::mclapply(
# 1:(m+1),
# function(i) residues[i] * solve(partial_fraction_terms[[i]], v),
# mc.cores = parallel::detectCores()
# )
#
# return(Reduce(`+`, results))
# }solve_fractional_evolution()Given the fractional operator object my_op_frac, a time
step time_step, a sequence of time points
time_seq, an initial value val_at_0, and the
right-hand side matrix RHST, function
solve_fractional_evolution() computes the solution to the
fractional diffusion equation at each time step using scheme \(\eqref{eq:final_scheme2}\).
solve_fractional_evolution <- function(my_op_frac, time_step, time_seq, val_at_0, RHST) {
CC <- my_op_frac$C
SOL <- matrix(NA, nrow = nrow(CC), ncol = length(time_seq))
SOL[, 1] <- val_at_0
for (k in 1:(length(time_seq) - 1)) {
rhs <- CC %*% SOL[, k] + time_step * RHST[, k + 1]
SOL[, k + 1] <- as.matrix(my.solver.frac(my_op_frac, rhs))
}
return(SOL)
}gets.graph.tadpole()Given a mesh size h, function
gets.graph.tadpole() builds a tadpole graph and creates a
mesh.
# Function to build a tadpole graph and create a mesh
gets.graph.tadpole <- function(h){
edge1 <- rbind(c(0,0),c(1,0))
theta <- seq(from=-pi,to=pi,length.out = 10000)
edge2 <- cbind(1+1/pi+cos(theta)/pi,sin(theta)/pi)
edges <- list(edge1, edge2)
graph <- metric_graph$new(edges = edges, verbose = 0)
graph$set_manual_edge_lengths(edge_lengths = c(1,2))
graph$build_mesh(h = h)
return(graph)
}tadpole.eig()Given a mode number k and a tadpole graph
graph, function tadpole.eig() computes the
eigenpairs of the tadpole graph.
# Function to compute the eigenfunctions of the tadpole graph
tadpole.eig <- function(k,graph){
x1 <- c(0,graph$get_edge_lengths()[1]*graph$mesh$PtE[graph$mesh$PtE[,1]==1,2])
x2 <- c(0,graph$get_edge_lengths()[2]*graph$mesh$PtE[graph$mesh$PtE[,1]==2,2])
if(k==0){
f.e1 <- rep(1,length(x1))
f.e2 <- rep(1,length(x2))
f1 = c(f.e1[1],f.e2[1],f.e1[-1], f.e2[-1])
f = list(phi=f1/sqrt(3))
} else {
f.e1 <- -2*sin(pi*k*1/2)*cos(pi*k*x1/2)
f.e2 <- sin(pi*k*x2/2)
f1 = c(f.e1[1],f.e2[1],f.e1[-1], f.e2[-1])
if((k %% 2)==1){
f = list(phi=f1/sqrt(3))
} else {
f.e1 <- (-1)^{k/2}*cos(pi*k*x1/2)
f.e2 <- cos(pi*k*x2/2)
f2 = c(f.e1[1],f.e2[1],f.e1[-1],f.e2[-1])
f <- list(phi=f1,psi=f2/sqrt(3/2))
}
}
return(f)
}gets.eigen.params()Given a finite number of modes N_finite, a scaling
parameter kappa, a smoothness parameter alpha,
and a tadpole graph graph, function
gets.eigen.params() computes EIGENVAL_ALPHA (a
vector with entries \(\lambda_j^{\alpha/2}\)),
EIGENVAL_MINUS_ALPHA (a vector with entries \(\lambda_j^{-\alpha/2}\)), and
EIGENFUN (a matrix with columns \(e_j\) on the mesh of
graph).
# Function to compute the eigenpairs of the tadpole graph
gets.eigen.params <- function(N_finite = 4, kappa = 1, alpha = 0.5, graph){
EIGENVAL <- NULL
EIGENVAL_ALPHA <- NULL
EIGENVAL_MINUS_ALPHA <- NULL
EIGENFUN <- NULL
INDEX <- NULL
for (j in 0:N_finite) {
lambda_j <- kappa^2 + (j*pi/2)^2
lambda_j_alpha_half <- lambda_j^(alpha/2)
lambda_j_minus_alpha_half <- lambda_j^(-alpha/2)
e_j <- tadpole.eig(j,graph)$phi
EIGENVAL <- c(EIGENVAL, lambda_j)
EIGENVAL_ALPHA <- c(EIGENVAL_ALPHA, lambda_j_alpha_half)
EIGENVAL_MINUS_ALPHA <- c(EIGENVAL_MINUS_ALPHA, lambda_j_minus_alpha_half)
EIGENFUN <- cbind(EIGENFUN, e_j)
INDEX <- c(INDEX, j)
if (j>0 && (j %% 2 == 0)) {
lambda_j <- kappa^2 + (j*pi/2)^2
lambda_j_alpha_half <- lambda_j^(alpha/2)
lambda_j_minus_alpha_half <- lambda_j^(-alpha/2)
e_j <- tadpole.eig(j,graph)$psi
EIGENVAL <- c(EIGENVAL, lambda_j)
EIGENVAL_ALPHA <- c(EIGENVAL_ALPHA, lambda_j_alpha_half)
EIGENVAL_MINUS_ALPHA <- c(EIGENVAL_MINUS_ALPHA, lambda_j_minus_alpha_half)
EIGENFUN <- cbind(EIGENFUN, e_j)
INDEX <- c(INDEX, j+0.1)
}
}
return(list(EIGENVAL = EIGENVAL,
EIGENVAL_ALPHA = EIGENVAL_ALPHA,
EIGENVAL_MINUS_ALPHA = EIGENVAL_MINUS_ALPHA,
EIGENFUN = EIGENFUN,
INDEX = INDEX))
}construct_piecewise_projection()Given a matrix projected_U_approx with approximated
values at discrete time points, a sequence of time points
time_seq, and an extended sequence of time points
overkill_time_seq, function
construct_piecewise_projection() constructs a piecewise
constant projection of the approximated values over the extended time
sequence.
# Function to construct a piecewise constant projection of approximated values
construct_piecewise_projection <- function(projected_U_approx, time_seq, overkill_time_seq) {
projected_U_piecewise <- matrix(NA, nrow = nrow(projected_U_approx), ncol = length(overkill_time_seq))
# Assign value at t = 0
projected_U_piecewise[, which(overkill_time_seq == 0)] <- projected_U_approx[, 1]
# Assign values for intervals (t_{k-1}, t_k]
for (k in 2:length(time_seq)) {
idxs <- which(overkill_time_seq > time_seq[k - 1] & overkill_time_seq <= time_seq[k])
projected_U_piecewise[, idxs] <- projected_U_approx[, k]
}
return(projected_U_piecewise)
}compute_guiding_lines()Given a vector x_axis_vector, a matrix
errors with error values, a vector
theoretical_rates with theoretical convergence rates, and a
function line_equation_fun (either
loglog_line_equation or exp_line_equation),
function compute_guiding_lines() computes guiding lines for
convergence plots.
loglog_line_equation <- function(x1, y1, slope) {
b <- log10(y1 / (x1 ^ slope))
function(x) {
(x ^ slope) * (10 ^ b)
}
}
exp_line_equation <- function(x1, y1, slope) {
lnC <- log(y1) - slope * x1
function(x) {
exp(lnC + slope * x)
}
}
compute_guiding_lines <- function(x_axis_vector, errors, theoretical_rates, line_equation_fun) {
guiding_lines <- matrix(NA, nrow = length(x_axis_vector), ncol = length(theoretical_rates))
for (j in seq_along(theoretical_rates)) {
guiding_lines_aux <- matrix(NA, nrow = length(x_axis_vector), ncol = length(x_axis_vector))
for (k in seq_along(x_axis_vector)) {
point_x1 <- x_axis_vector[k]
point_y1 <- errors[k, j]
slope <- theoretical_rates[j]
line <- line_equation_fun(x1 = point_x1, y1 = point_y1, slope = slope)
guiding_lines_aux[, k] <- line(x_axis_vector)
}
guiding_lines[, j] <- rowMeans(guiding_lines_aux)
}
return(guiding_lines)
}Below we present closed-form expressions for \(\displaystyle G_j(t)= \int_0^t e^{\lambda^{\alpha/2}_jr}g(r)dr\) corresponding to some choices of \(g(r)\).
\[\begin{aligned} g(r) &= Ae^{-\lambda^{\alpha/2}_j r} &\implies\quad G_j(t) &= At, \quad A\in\mathbb{R} \\[1.5ex] g(r) &= Ae^{\mu r} &\implies\quad G_j(t) &= A \frac{e^{(\lambda^{\alpha/2}_j+\mu)t} - 1}{\lambda^{\alpha/2}_j + \mu}, \quad -\lambda^{\alpha/2}_j \ne \mu \in \mathbb{R} \\[1.5ex] g(r) &= Ar^n &\implies\quad G_j(t) &= A\frac{(-1)^{n+1} n!}{(\lambda_j^{\alpha/2})^{n+1}} \left(1 - e^{\lambda_j^{\alpha/2} t} \sum_{k=0}^n \frac{(-\lambda_j^{\alpha/2} t)^k}{k!} \right), \quad n=0,1,\dots \\[1.5ex] g(r) &= A\sin(\omega r) &\implies\quad G_j(t) &= A \frac{e^{\lambda_j^{\alpha/2} t} \left( \lambda_j^{\alpha/2} \sin(\omega t) - \omega \cos(\omega t) \right) + \omega}{(\lambda_j^{\alpha/2})^2 + \omega^2}, \quad \omega \in \mathbb{R} \\[1.5ex] g(r) &= A\cos(\theta r) &\implies\quad G_j(t) &= A \frac{e^{\lambda_j^{\alpha/2} t} \left( \lambda_j^{\alpha/2} \cos(\theta t) + \theta \sin(\theta t) \right) - \lambda_j^{\alpha/2}}{(\lambda_j^{\alpha/2})^2 + \theta^2}, \quad \theta \in \mathbb{R} \end{aligned}\]# Functions to compute the exact solution to the fractional diffusion equation
g_linear <- function(r, A, lambda_j_alpha_half) {
return(A * exp(-lambda_j_alpha_half * r))
}
G_linear <- function(t, A) {
return(A * t)
}
g_exp <- function(r, A, mu) {
return(A * exp(mu * r))
}
G_exp <- function(t, A, lambda_j_alpha_half, mu) {
exponent <- lambda_j_alpha_half + mu
return(A * (exp(exponent * t) - 1) / exponent)
}
g_poly <- function(r, A, n) {
return(A * r^n)
}
G_poly <- function(t, A, lambda_j_alpha_half, n) {
t <- as.vector(t)
k_vals <- 0:n
sum_term <- sapply(t, function(tt) {
sum(((-lambda_j_alpha_half * tt)^k_vals) / factorial(k_vals))
})
coeff <- ((-1)^(n + 1)) * factorial(n) / (lambda_j_alpha_half^(n + 1))
return(A * coeff * (1 - exp(lambda_j_alpha_half * t) * sum_term))
}
g_sin <- function(r, A, omega) {
return(A * sin(omega * r))
}
G_sin <- function(t, A, lambda_j_alpha_half, omega) {
denom <- lambda_j_alpha_half^2 + omega^2
numerator <- exp(lambda_j_alpha_half * t) * (lambda_j_alpha_half * sin(omega * t) - omega * cos(omega * t)) + omega
return(A * numerator / denom)
}
g_cos <- function(r, A, theta) {
return(A * cos(theta * r))
}
G_cos <- function(t, A, lambda_j_alpha_half, theta) {
denom <- lambda_j_alpha_half^2 + theta^2
numerator <- exp(lambda_j_alpha_half * t) * (lambda_j_alpha_half * cos(theta * t) + theta * sin(theta * t)) - lambda_j_alpha_half
return(A * numerator / denom)
}gets.mesh.and.FEM.on.rectangle <- function(a, b, n) {
mesh <- fmesher::fm_rcdt_2d(
lattice = fmesher::fm_lattice_2d(
x = seq(0, a, length.out = n + 1),
y = seq(0, b, length.out = n + 1)
), extend = FALSE)
FEM <- fmesher::fm_fem(mesh, order = 1)
return(list(mesh = mesh,
Cl = FEM$c0,
C = FEM$c1,
G = FEM$g1))
}compute_true_eigen_rectangle <- function(a = 1,
b = 1,
loc,
kappa = 1,
alpha = 1,
m_max = 3,
n_max = 3) {
loc_x <- loc[,1]
loc_y <- loc[,2]
N <- (m_max+1)*(n_max+1)
# Prepare lists
EIGENVAL <- numeric(N)
m_vector <- numeric(N)
n_vector <- numeric(N)
EIGENFUN <- matrix(0, nrow = nrow(loc), ncol = N)
i <- 0
# Loop over mode indices
for (m in 0:m_max) {
for (n in 0:n_max) {
lambda_mn <- kappa^2 + pi^2 * ((m^2 / a^2) + (n^2 / b^2))
EIGENVAL[i + 1] <- lambda_mn
# Evaluate eigenfunction on mesh grid
phi_mn <- cos(m*pi*loc_x/a) * cos(n*pi*loc_y/b)
EIGENFUN[, i + 1] <- phi_mn
m_vector[i + 1] <- m
n_vector[i + 1] <- n
i <- i + 1
}
}
# Sort eigenvalues and corresponding eigenfunctions
idx <- order(EIGENVAL)
EIGENVAL <- EIGENVAL[idx]
EIGENVAL_ALPHA <- EIGENVAL^(alpha/2)
EIGENFUN <- EIGENFUN[, idx]
m_vector <- m_vector[idx]
n_vector <- n_vector[idx]
return(list(EIGENVAL = EIGENVAL,
EIGENVAL_ALPHA = EIGENVAL_ALPHA,
EIGENFUN = EIGENFUN,
m_vector = m_vector,
n_vector = n_vector))
}# Function to get globe from h using spline
globe_from_h <- function(h) {
readed_data <- readRDS(
file = here::here("data_files/h_min_max_vs_globe.RDS")
)
globe_vector <- readed_data$globe_vector
h_min_vector <- readed_data$h_min_vector
ord <- order(h_min_vector)
h_sorted <- h_min_vector[ord]
globe_sorted <- globe_vector[ord]
# Create a smooth spline function: globe as function of h_max
spline_func <- splinefun(x = h_sorted,
y = globe_sorted,
method = "monoH.FC")
return(round(spline_func(h)))
}calculate_laplace_beltrami_eigenvalues <- function(kappa = 0, L_max = 5, rot.inv = FALSE) {
eigenvalues <- numeric(0)
for (l in 0:L_max) {
lambda_l <- kappa^2 + l * (l + 1)
if (rot.inv) {
# Only one eigenvalue per l (rotational invariance)
eigenvalues <- c(eigenvalues, lambda_l)
} else {
# Full multiplicity: (m = -l,...,l)
multiplicity <- 2 * l + 1
eigenvalues <- c(eigenvalues, rep(lambda_l, multiplicity))
}
}
return(eigenvalues)
}
compute_true_eigen_sphere <- function(mesh,
kappa,
alpha,
L_max,
rot.inv){
EIGENFUN <- fmesher::fm_raw_basis(
mesh = mesh,
type = "sph.harm",
n = L_max,
rot.inv = rot.inv)
EIGENVAL <- calculate_laplace_beltrami_eigenvalues(
kappa = kappa,
L_max = L_max,
rot.inv = rot.inv)
idx <- order(EIGENVAL)
EIGENVAL <- EIGENVAL[idx]
EIGENVAL_ALPHA <- EIGENVAL^(alpha/2)
EIGENFUN <- EIGENFUN[, idx]
return(list(EIGENVAL = EIGENVAL,
EIGENVAL_ALPHA = EIGENVAL_ALPHA,
EIGENFUN = EIGENFUN))
}plotting.order()Given a vector v and a graph object graph,
function plotting.order() orders the mesh values for
plotting.
global.scene.setter()Given ranges for the x, y, and
z axes, and an optional aspect ratio for the z
axis, function global.scene.setter() sets the scene for 3D
plots so that all plots have the same aspect ratio and camera
position.
# Function to set the scene for 3D plots
global.scene.setter <- function(x_range, y_range, z_range, z_aspectratio = 4) {
return(list(xaxis = list(title = "x", range = x_range),
yaxis = list(title = "y", range = y_range),
zaxis = list(title = "z", range = z_range),
aspectratio = list(x = 2*(1+2/pi),
y = 2*(2/pi),
z = z_aspectratio*(2/pi)),
camera = list(eye = list(x = (1+2/pi)/2,
y = 4,
z = 2),
center = list(x = (1+2/pi)/2,
y = 0,
z = 0))))
}global.scene.setter.rec.and.sphere <- function(x_range, y_range, z_range) {
return(list(xaxis = list(title = "x", range = x_range),
yaxis = list(title = "y", range = y_range),
zaxis = list(title = "z", range = z_range),
aspectratio = list(x = 1,
y = 1,
z = 1),
camera = list(eye = list(x = 2,
y = 2,
z = 2),
center = list(x = 0,
y = 0,
z = 0))))
}
global.scene.setter.rec.and.sphere.for.hat <- function(x_range, y_range, z_range) {
return(list(xaxis = list(title = "x", range = x_range),
yaxis = list(title = "y", range = y_range),
zaxis = list(title = "z", range = z_range),
aspectratio = list(x = 1,
y = 1,
z = 0.4),
camera = list(eye = list(x = 2,
y = 2,
z = 2),
center = list(x = 0,
y = 0,
z = 0))))
}graph.plotter.3d()Given a graph object graph, a sequence of time points
time_seq, and one or more matrices ...
representing function values defined on the mesh of graph
at each time in time_seq, the
graph.plotter.3d() function generates an interactive 3D
visualization of these values over time.
# Function to plot in 3D
graph.plotter.3d.old <- function(graph, time_seq, frame_val_to_display, ...) {
U_list <- list(...)
U_names <- sapply(substitute(list(...))[-1], deparse)
# Spatial coordinates
x <- plotting.order(graph$mesh$V[, 1], graph)
y <- plotting.order(graph$mesh$V[, 2], graph)
weights <- graph$mesh$weights
# Apply plotting.order to each U
U_list <- lapply(U_list, function(U) apply(U, 2, plotting.order, graph = graph))
n_vars <- length(U_list)
# Create plot_data frame with time and position replicated
n_time <- ncol(U_list[[1]])
base_data <- data.frame(
x = rep(x, times = n_time),
y = rep(y, times = n_time),
the_graph = 0,
frame = rep(time_seq, each = length(x))
)
# Add U columns to plot_data
for (i in seq_along(U_list)) {
base_data[[paste0("u", i)]] <- as.vector(U_list[[i]])
}
plot_data <- base_data
# Generate vertical lines
vertical_lines_list <- lapply(seq_along(U_list), function(i) {
do.call(rbind, lapply(time_seq, function(t) {
idx <- which(plot_data$frame == t)
z_vals <- plot_data[[paste0("u", i)]][idx]
data.frame(
x = rep(plot_data$x[idx], each = 3),
y = rep(plot_data$y[idx], each = 3),
z = as.vector(t(cbind(0, z_vals, NA))),
frame = rep(t, each = length(idx) * 3)
)
}))
})
# Set axis ranges
z_range <- range(unlist(U_list))
x_range <- range(x)
y_range <- range(y)
# Create plot
p <- plot_ly(plot_data, frame = ~frame) %>%
add_trace(x = ~x, y = ~y, z = ~the_graph, type = "scatter3d", mode = "lines",
name = "", showlegend = FALSE,
line = list(color = "black", width = 3))
# Add traces for each variable
colors <- rev(viridisLite::viridis(n_vars)) #RColorBrewer::brewer.pal(min(n_vars, 8), "Set1")
for (i in seq_along(U_list)) {
p <- add_trace(p,
x = ~x, y = ~y, z = as.formula(paste0("~u", i)),
type = "scatter3d", mode = "lines", name = U_names[i],
line = list(color = colors[i], width = 3))
}
# Add vertical lines
for (i in seq_along(vertical_lines_list)) {
p <- add_trace(p,
data = vertical_lines_list[[i]],
x = ~x, y = ~y, z = ~z, frame = ~frame,
type = "scatter3d", mode = "lines",
line = list(color = "gray", width = 0.5),
name = "Vertical lines",
showlegend = FALSE)
}
frame_name <- deparse(substitute(frame_val_to_display))
# Layout and animation controls
p <- p %>%
layout(
scene = global.scene.setter(x_range, y_range, z_range),
updatemenus = list(list(type = "buttons", showactive = FALSE,
buttons = list(
list(label = "Play", method = "animate",
args = list(NULL, list(frame = list(duration = 2000 / length(time_seq), redraw = TRUE), fromcurrent = TRUE))),
list(label = "Pause", method = "animate",
args = list(NULL, list(mode = "immediate", frame = list(duration = 0), redraw = FALSE)))
)
)),
title = paste0(frame_name,": ", formatC(frame_val_to_display[1], format = "f", digits = 4))
) %>%
plotly_build()
for (i in seq_along(p$x$frames)) {
p$x$frames[[i]]$layout <- list(title = paste0(frame_name,": ", formatC(frame_val_to_display[i], format = "f", digits = 4)))
}
return(p)
}graph.plotter.3d()Given a graph object graph, a sequence of time points
time_seq, a vector frame_val_to_display with
values to display at each frame, and a list U_list
representing function values defined on the mesh of graph
at each time in time_seq, the
graph.plotter.3d() function generates an interactive 3D
visualization of these functions over time.
graph.plotter.3d <- function(graph, time_seq, frame_val_to_display, U_list) {
U_names <- names(U_list)
# Spatial coordinates
x <- plotting.order(graph$mesh$V[, 1], graph)
y <- plotting.order(graph$mesh$V[, 2], graph)
weights <- graph$mesh$weights
# Apply plotting.order to each U
U_list <- lapply(U_list, function(U) apply(U, 2, plotting.order, graph = graph))
n_vars <- length(U_list)
# Create plot_data frame with time and position replicated
n_time <- ncol(U_list[[1]])
base_data <- data.frame(
x = rep(x, times = n_time),
y = rep(y, times = n_time),
the_graph = 0,
frame = rep(time_seq, each = length(x))
)
# Add U columns to plot_data
for (i in seq_along(U_list)) {
base_data[[paste0("u", i)]] <- as.vector(U_list[[i]])
}
plot_data <- base_data
# Generate vertical lines
vertical_lines_list <- lapply(seq_along(U_list), function(i) {
do.call(rbind, lapply(time_seq, function(t) {
idx <- which(plot_data$frame == t)
z_vals <- plot_data[[paste0("u", i)]][idx]
data.frame(
x = rep(plot_data$x[idx], each = 3),
y = rep(plot_data$y[idx], each = 3),
z = as.vector(t(cbind(0, z_vals, NA))),
frame = rep(t, each = length(idx) * 3)
)
}))
})
# Set axis ranges
z_range <- range(unlist(U_list))
x_range <- range(x)
y_range <- range(y)
# Create plot
p <- plot_ly(plot_data, frame = ~frame) %>%
add_trace(x = ~x, y = ~y, z = ~the_graph, type = "scatter3d", mode = "lines",
name = "", showlegend = FALSE,
line = list(color = "black", width = 3))
if (n_vars == 2) {
colors <- RColorBrewer::brewer.pal(min(n_vars, 8), "Set1")
} else {
colors <- rev(viridisLite::viridis(n_vars))
}
# RColorBrewer::brewer.pal(min(n_vars, 8), "Set1")
for (i in seq_along(U_list)) {
p <- add_trace(p,
x = ~x, y = ~y, z = as.formula(paste0("~u", i)),
type = "scatter3d", mode = "lines", name = U_names[i],
line = list(color = colors[i], width = 3))
}
# Add vertical lines
for (i in seq_along(vertical_lines_list)) {
p <- add_trace(p,
data = vertical_lines_list[[i]],
x = ~x, y = ~y, z = ~z, frame = ~frame,
type = "scatter3d", mode = "lines",
line = list(color = "gray", width = 0.5),
name = "Vertical lines",
showlegend = FALSE)
}
frame_name <- deparse(substitute(frame_val_to_display))
# Layout and animation controls
p <- p %>%
layout(
scene = global.scene.setter(x_range, y_range, z_range),
updatemenus = list(list(type = "buttons", showactive = FALSE,
buttons = list(
list(label = "Play", method = "animate",
args = list(NULL, list(frame = list(duration = 2000 / length(time_seq), redraw = TRUE), fromcurrent = TRUE))),
list(label = "Pause", method = "animate",
args = list(NULL, list(mode = "immediate", frame = list(duration = 0), redraw = FALSE)))
)
)),
title = paste0(frame_name,": ", formatC(frame_val_to_display[1], format = "f", digits = 4))
) %>%
plotly_build()
for (i in seq_along(p$x$frames)) {
p$x$frames[[i]]$layout <- list(title = paste0(frame_name,": ", formatC(frame_val_to_display[i], format = "f", digits = 4)))
}
return(p)
}error.at.each.time.plotter()Given a graph object graph, a matrix U_true
of true values, a matrix U_approx of approximated values, a
sequence of time points time_seq, and a time step
time_step, function
error.at.each.time.plotter() computes the error at each
time step and generates a plot showing the error over time.
# Function to plot the error at each time step
error.at.each.time.plotter <- function(graph, U_true, U_approx, time_seq, time_step) {
weights <- graph$mesh$weights
error_at_each_time <- t(weights) %*% (U_true - U_approx)^2
error <- sqrt(as.double(t(weights) %*% (U_true - U_approx)^2 %*% rep(time_step, ncol(U_true))))
p <- plot_ly() %>%
add_trace(
x = ~time_seq, y = ~error_at_each_time, type = 'scatter', mode = 'lines+markers',
line = list(color = 'blue', width = 2),
marker = list(color = 'blue', size = 4),
name = "",
showlegend = TRUE
) %>%
layout(
title = paste0("Error at Each Time Step (Total error = ", formatC(error, format = "f", digits = 9), ")"),
xaxis = list(title = "t"),
yaxis = list(title = "Error"),
legend = list(x = 0.1, y = 0.9)
)
return(p)
}graph.plotter.3d.comparer()Given a graph object graph, matrices U_true
and U_approx representing true and approximated values, and
a sequence of time points time_seq, function
graph.plotter.3d.comparer() generates a 3D plot comparing
the true and approximated values over time, with color-coded traces for
each time point.
# Function to plot the 3D comparison of U_true and U_approx
graph.plotter.3d.comparer <- function(graph, U_true, U_approx, time_seq) {
x <- graph$mesh$V[, 1]; y <- graph$mesh$V[, 2]
x <- plotting.order(x, graph); y <- plotting.order(y, graph)
U_true <- apply(U_true, 2, plotting.order, graph = graph)
U_approx <- apply(U_approx, 2, plotting.order, graph = graph)
n_times <- length(time_seq)
x_range <- range(x); y_range <- range(y); z_range <- range(c(U_true, U_approx))
# Normalize time_seq
time_normalized <- (time_seq - min(time_seq)) / (max(time_seq) - min(time_seq))
blues <- colorRampPalette(c("lightblue", "blue"))(n_times)
reds <- colorRampPalette(c("mistyrose", "red"))(n_times)
# Accurate colorscales
colorscale_greens <- Map(function(t, col) list(t, col), time_normalized, blues)
colorscale_reds <- Map(function(t, col) list(t, col), time_normalized, reds)
p <- plot_ly()
# Static black graph structure
p <- p %>%
add_trace(x = x, y = y, z = rep(0, length(x)),
type = "scatter3d", mode = "lines",
line = list(color = "black", width = 4),
name = "Graph", showlegend = FALSE)
# U_true traces (green)
for (i in seq_len(n_times)) {
z <- U_true[, i]
p <- add_trace(
p,
type = "scatter3d",
mode = "lines",
x = x, y = y, z = z,
line = list(color = blues[i], width = 4),
showlegend = FALSE,
scene = "scene"
)
}
# U_approx traces (dashed red)
for (i in seq_len(n_times)) {
z <- U_approx[, i]
p <- add_trace(
p,
type = "scatter3d",
mode = "lines",
x = x, y = y, z = z,
line = list(color = reds[i], width = 4, dash = "dot"),
showlegend = FALSE,
scene = "scene"
)
}
# Dummy green colorbar (True) – with ticks
p <- add_trace(
p,
type = "heatmap",
z = matrix(time_seq, nrow = 1),
showscale = TRUE,
colorscale = colorscale_greens,
colorbar = list(
title = list(font = list(size = 12, color = "black"), text = "Time", side = "top"),
len = 0.9,
thickness = 15,
x = 1.02,
xanchor = "left",
y = 0.5,
yanchor = "middle",
tickvals = NULL, # hide tick values
ticktext = NULL,
ticks = "" # also hides tick marks
),
x = matrix(time_seq, nrow = 1),
y = matrix(1, nrow = 1),
hoverinfo = "skip",
opacity = 0
)
# Dummy red colorbar (Approx) – no ticks
p <- add_trace(
p,
type = "heatmap",
z = matrix(time_seq, nrow = 1),
showscale = TRUE,
colorscale = colorscale_reds,
colorbar = list(
title = list(font = list(size = 12, color = "black"), text = ".", side = "top"),
len = 0.9,
thickness = 15,
x = 1.05,
xanchor = "left",
y = 0.5,
yanchor = "middle"
),
x = matrix(time_seq, nrow = 1),
y = matrix(1, nrow = 1),
hoverinfo = "skip",
opacity = 0
)
p <- p %>%
add_trace(x = x, y = y, z = rep(0, length(x)),
type = "scatter3d", mode = "lines",
line = list(color = "black", width = 4),
name = "Graph", showlegend = FALSE)
p <- layout(p,
scene = global.scene.setter(x_range, y_range, z_range),
xaxis = list(visible = FALSE),
yaxis = list(visible = FALSE),
annotations = list(
list(
text = "Exact",
x = 1.045,
y = 0.5,
xref = "paper",
yref = "paper",
showarrow = FALSE,
font = list(size = 12, color = "black"),
textangle = -90
),
list(
text = "Approx",
x = 1.075,
y = 0.5,
xref = "paper",
yref = "paper",
showarrow = FALSE,
font = list(size = 12, color = "black"),
textangle = -90
)
)
)
return(p)
}graph.plotter.3d.single()Given a graph object graph, a matrix U_true
representing true values, and a sequence of time points
time_seq, function graph.plotter.3d.single()
generates a 3D plot of the true values over time, with color-coded
traces for each time point.
# Function to plot a single 3D line for
graph.plotter.3d.single <- function(graph, U_true, time_seq) {
x <- graph$mesh$V[, 1]; y <- graph$mesh$V[, 2]
x <- plotting.order(x, graph); y <- plotting.order(y, graph)
U_true <- apply(U_true, 2, plotting.order, graph = graph)
n_times <- length(time_seq)
x_range <- range(x); y_range <- range(y); z_range <- range(U_true)
z_range[1] <- z_range[1] - 10^-6
viridis_colors <- viridisLite::viridis(100)
# Normalize time_seq
time_normalized <- (time_seq - min(time_seq)) / (max(time_seq) - min(time_seq))
#greens <- colorRampPalette(c("palegreen", "darkgreen"))(n_times)
greens <- colorRampPalette(c(viridis_colors[1], viridis_colors[50], viridis_colors[100]))(n_times)
# Accurate colorscales
colorscale_greens <- Map(function(t, col) list(t, col), time_normalized, greens)
p <- plot_ly()
# Add the 3D lines with fading green color
for (i in seq_len(n_times)) {
z <- U_true[, i]
p <- add_trace(
p,
type = "scatter3d",
mode = "lines",
x = x,
y = y,
z = z,
line = list(color = greens[i], width = 2),
showlegend = FALSE,
scene = "scene"
)
}
p <- p %>%
add_trace(x = x, y = y, z = rep(0, length(x)),
type = "scatter3d", mode = "lines",
line = list(color = "black", width = 5),
name = "Graph", showlegend = FALSE)
# Add dummy heatmap to show colorbar (not part of scene)
p <- add_trace(
p,
type = "heatmap",
z = matrix(time_seq, nrow = 1),
showscale = TRUE,
colorscale = colorscale_greens,
colorbar = list(
title = list(font = list(size = 12, color = "black"), text = "Time", side = "top"),
len = 0.9, # height (0 to 1)
thickness = 15, # width in pixels
x = 1.02, # shift it slightly right of the plot
xanchor = "left",
y = 0.5,
yanchor = "middle"),
x = matrix(time_seq, nrow = 1),
y = matrix(1, nrow = 1),
hoverinfo = "skip",
opacity = 0
)
p <- layout(p,
scene = global.scene.setter(x_range, y_range, z_range),
xaxis = list(visible = FALSE),
yaxis = list(visible = FALSE)
)
return(p)
}error.convergence.plotter()Given a vector x_axis_vector, a vector
alpha_vector of parameter values, a matrix
errors with error values, vectors
theoretical_rates and observed_rates with
convergence rates, a function line_equation_fun (either
loglog_line_equation or exp_line_equation),
and plot titles and labels, function
error.convergence.plotter() generates a convergence plot
showing errors and guiding lines.
# Function to plot the error convergence
error.convergence.plotter <- function(x_axis_vector,
alpha_vector,
errors,
theoretical_rates,
observed_rates,
line_equation_fun,
fig_title,
x_axis_label,
apply_sqrt = FALSE) {
x_vec <- if (apply_sqrt) sqrt(x_axis_vector) else x_axis_vector
guiding_lines <- compute_guiding_lines(x_axis_vector = x_vec,
errors = errors,
theoretical_rates = theoretical_rates,
line_equation_fun = line_equation_fun)
default_colors <- scales::hue_pal()(length(alpha_vector))
plot_lines <- lapply(1:ncol(guiding_lines), function(i) {
geom_line(
data = data.frame(x = x_vec, y = guiding_lines[, i]),
aes(x = x, y = y),
color = default_colors[i],
linetype = "dashed",
show.legend = FALSE
)
})
df <- as.data.frame(cbind(x_vec, errors))
colnames(df) <- c("x_axis_vector", alpha_vector)
df_melted <- melt(df, id.vars = "x_axis_vector", variable.name = "column", value.name = "value")
custom_labels <- paste0(formatC(alpha_vector, format = "f", digits = 2),
" | ",
formatC(theoretical_rates, format = "f", digits = 4),
" | ",
formatC(observed_rates, format = "f", digits = 4))
df_melted$column <- factor(df_melted$column, levels = alpha_vector, labels = custom_labels)
p <- ggplot() +
geom_line(data = df_melted, aes(x = x_axis_vector, y = value, color = column)) +
geom_point(data = df_melted, aes(x = x_axis_vector, y = value, color = column)) +
plot_lines +
labs(
title = fig_title,
x = x_axis_label,
y = expression(Error),
color = " α | theo | obs"
) +
(if (apply_sqrt) {
scale_x_continuous(breaks = x_vec, labels = round(x_axis_vector, 4))
} else {
scale_x_log10(breaks = x_axis_vector, labels = round(x_axis_vector, 4))
}) +
(if (apply_sqrt) {
scale_y_continuous(trans = "log", labels = scales::scientific_format())
} else {
scale_y_log10(labels = scales::scientific_format())
}) +
theme_minimal() +
theme(text = element_text(family = "Palatino"),
legend.position = "bottom",
legend.direction = "vertical",
plot.margin = margin(0, 0, 0, 0),
plot.title = element_text(hjust = 0.5, size = 18, face = "bold"))
return(p)
}graph.plotter.3d.static()Given a graph object graph and a list
z_list of function values defined on the mesh of
graph, function graph.plotter.3d.static()
generates a static 3D plot of these values.
graph.plotter.3d.static <- function(graph, z_list) {
x <- plotting.order(graph$mesh$V[, 1], graph)
y <- plotting.order(graph$mesh$V[, 2], graph)
U_names <- names(z_list)
n_vars <- length(z_list)
z_list <- lapply(z_list, function(z) plotting.order(z, graph))
# Axis ranges
z_range <- range(unlist(z_list))
x_range <- range(x)
y_range <- range(y)
if (n_vars == 2) {
colors <- RColorBrewer::brewer.pal(min(n_vars, 8), "Set1")
} else {
colors <- rev(viridisLite::viridis(n_vars))
}
p <- plot_ly()
for (i in seq_along(z_list)) {
z <- z_list[[i]]
# Main 3D curve
p <- add_trace(
p,
x = x, y = y, z = z,
type = "scatter3d", mode = "lines",
line = list(color = colors[i], width = 3),
name = U_names[i], showlegend = TRUE
)
# Efficient vertical lines: one trace with breaks (NA)
x_vert <- rep(x, each = 3)
y_vert <- rep(y, each = 3)
z_vert <- unlist(lapply(z, function(zj) c(0, zj, NA)))
p <- add_trace(
p,
x = x_vert, y = y_vert, z = z_vert,
type = "scatter3d", mode = "lines",
line = list(color = "gray", width = 0.5),
showlegend = FALSE
)
}
p <- p %>% add_trace(x = x, y = y, z = x*0, type = "scatter3d", mode = "lines",
line = list(color = "black", width = 3),
name = "thegraph", showlegend = FALSE) %>%
layout(scene = global.scene.setter(x_range, y_range, z_range))
return(p)
}graph.plotter.3d.two.meshes.time()Given two graph objects graph_finer and
graph_coarser, sequences of time points
time_seq and frame_val_to_display, and lists
fs_finer and fs_coarser of function values
defined on the meshes of graph_finer and
graph_coarser respectively, function
graph.plotter.3d.two.meshes.time() generates an interactive
3D visualization comparing these functions over time.
graph.plotter.3d.two.meshes.time <- function(graph_finer, graph_coarser,
time_seq, frame_val_to_display,
fs_finer = list(), fs_coarser = list()) {
# Spatial coordinates (ordered for plotting)
x_finer <- plotting.order(graph_finer$mesh$V[, 1], graph_finer)
y_finer <- plotting.order(graph_finer$mesh$V[, 2], graph_finer)
x_coarser <- plotting.order(graph_coarser$mesh$V[, 1], graph_coarser)
y_coarser <- plotting.order(graph_coarser$mesh$V[, 2], graph_coarser)
n_time <- if (length(fs_finer) > 0) ncol(fs_finer[[1]]) else ncol(fs_coarser[[1]])
# Helper: make dataframe from one function
make_df <- function(f_mat, graph, x, y, mesh_name) {
z <- apply(f_mat, 2, plotting.order, graph = graph)
data.frame(
x = rep(x, times = n_time),
y = rep(y, times = n_time),
z = as.vector(z),
frame = rep(time_seq, each = length(x)),
mesh = mesh_name
)
}
# Build data for finer functions
data_finer_list <- lapply(names(fs_finer), function(nm) {
make_df(fs_finer[[nm]], graph_finer, x_finer, y_finer, nm)
})
# Build data for coarser functions
data_coarser_list <- lapply(names(fs_coarser), function(nm) {
make_df(fs_coarser[[nm]], graph_coarser, x_coarser, y_coarser, nm)
})
# Combine
all_data <- c(data_finer_list, data_coarser_list)
# Baseline graph (on finer mesh for consistency)
data_graph <- data.frame(
x = rep(x_finer, times = n_time),
y = rep(y_finer, times = n_time),
z = 0,
frame = rep(time_seq, each = length(x_finer)),
mesh = "Graph"
)
# --------- Vertical lines helper ----------
vertical_lines <- function(x, y, z, frame_vals, mesh_name) {
do.call(rbind, lapply(seq_along(frame_vals), function(i) {
idx <- ((i - 1) * length(x) + 1):(i * length(x))
data.frame(
x = rep(x, each = 3),
y = rep(y, each = 3),
z = as.vector(t(cbind(0, z[idx], NA))),
frame = rep(frame_vals[i], each = length(x) * 3),
mesh = mesh_name
)
}))
}
# --------- Compute vertical lines per mesh using max absolute value ---------
make_vertical_from_list <- function(data_list, x, y, mesh_name) {
if (length(data_list) == 0) return(NULL)
# Reshape each function's z back to matrix: (nodes × time)
z_mats <- lapply(data_list, function(df) {
matrix(df$z, nrow = length(x), ncol = length(time_seq))
})
# Stack into 3D array: (nodes × time × functions)
arr <- array(unlist(z_mats), dim = c(length(x), length(time_seq), length(z_mats)))
# For each node × time, select the entry with largest absolute value (keep sign)
idx <- apply(arr, c(1, 2), function(v) which.max(abs(v)))
z_signed_max <- mapply(function(i, j) arr[i, j, idx[i, j]],
rep(1:length(x), times = length(time_seq)),
rep(1:length(time_seq), each = length(x)))
# Flatten back into long vector
z_signed_max <- as.vector(z_signed_max)
vertical_lines(x, y, z_signed_max, time_seq, mesh_name)
}
vertical_finer <- make_vertical_from_list(data_finer_list, x_finer, y_finer, "finer")
vertical_coarser <- make_vertical_from_list(data_coarser_list, x_coarser, y_coarser, "coarser")
# Compute ranges
all_z <- unlist(lapply(all_data, function(df) df$z))
x_range <- range(c(x_finer, x_coarser))
y_range <- range(c(y_finer, y_coarser))
z_range <- range(all_z)
# --------- Plotly object ----------
p <- plot_ly(frame = ~frame)
# Add traces for finer + coarser (looping automatically with names)
for (df in all_data) {
p <- p %>%
add_trace(data = df,
x = ~x, y = ~y, z = ~z,
type = "scatter3d", mode = "lines",
line = list(width = 3),
name = unique(df$mesh))
}
# Add baseline
p <- p %>%
add_trace(data = data_graph,
x = ~x, y = ~y, z = ~z,
type = "scatter3d", mode = "lines",
line = list(color = "black", width = 2),
name = "Graph", showlegend = FALSE)
# Add verticals (one per mesh, envelope of all functions)
if (!is.null(vertical_finer)) {
p <- p %>%
add_trace(data = vertical_finer,
x = ~x, y = ~y, z = ~z,
type = "scatter3d", mode = "lines",
line = list(color = "gray", width = 0.5),
name = "Vertical finer", showlegend = FALSE)
}
if (!is.null(vertical_coarser)) {
p <- p %>%
add_trace(data = vertical_coarser,
x = ~x, y = ~y, z = ~z,
type = "scatter3d", mode = "lines",
line = list(color = "gray", width = 0.5),
name = "Vertical coarser", showlegend = FALSE)
}
frame_name <- deparse(substitute(frame_val_to_display))
p <- p %>%
layout(
scene = global.scene.setter(x_range, y_range, z_range),
updatemenus = list(list(
type = "buttons", showactive = FALSE,
buttons = list(
list(label = "Play", method = "animate",
args = list(NULL, list(frame = list(duration = 2000 / length(time_seq), redraw = TRUE),
fromcurrent = TRUE))),
list(label = "Pause", method = "animate",
args = list(NULL, list(mode = "immediate",
frame = list(duration = 0), redraw = FALSE)))
)
)),
title = paste0(frame_name, ": ", formatC(frame_val_to_display[1], format = "f", digits = 4))
) %>%
plotly_build()
# Update frame titles
for (i in seq_along(p$x$frames)) {
p$x$frames[[i]]$layout <- list(
title = paste0(frame_name, ": ", formatC(frame_val_to_display[i], format = "f", digits = 4))
)
}
return(p)
}plot_3d_rectangle_scatter()plot_3d_rectangle_scatter <- function(loc, eigvals, eigfuncs, fixed_colorscale = TRUE) {
x <- loc[,1]
y <- loc[,2]
colorscale <- "Viridis"
# Compute global color limits if fixed
if (fixed_colorscale) {
cmin <- min(eigfuncs)
cmax <- max(eigfuncs)
} else {
cmin <- NULL
cmax <- NULL
}
# Frames ----------------------------------------------------------------
frames <- lapply(seq_len(ncol(eigfuncs)), function(i) {
list(
name = as.character(i),
data = list(list(
x = x,
y = y,
z = eigfuncs[,i],
type = "scatter3d",
mode = "markers",
marker = list(
size = 5,
color = eigfuncs[,i],
colorscale = colorscale,
showscale = TRUE,
cmin = cmin,
cmax = cmax
)
))
)
})
# Initial Plot ----------------------------------------------------------
p <- plot_ly(
x = x,
y = y,
z = eigfuncs[,1],
type = "scatter3d",
mode = "markers",
marker = list(
size = 5,
color = eigfuncs[,1],
colorscale = colorscale,
showscale = TRUE,
cmin = cmin,
cmax = cmax
),
frame = "1"
)
p$x$frames <- frames
z_range <- range(eigfuncs)
x_range <- range(x)
y_range <- range(y)
frame_name <- deparse(substitute(eigvals))
# Layout + slider -------------------------------------------------------
p <- p %>% layout(
title = paste0(frame_name, ": ", eigvals[1]),
scene = global.scene.setter.rec.and.sphere(x_range, y_range, z_range),
sliders = list(
list(
active = 0,
currentvalue = list(prefix = "Frame: "),
pad = list(t = 50),
steps = lapply(seq_len(ncol(eigfuncs)), function(i) {
list(
label = as.character(i),
method = "animate",
args = list(
list(as.character(i)),
list(frame = list(duration = 300, redraw = TRUE),
mode = "immediate")
)
)
})
)
),
updatemenus = list(
list(
type = "buttons",
showactive = FALSE,
y = 1,
x = 1.15,
xanchor = "right",
yanchor = "top",
buttons = list(
list(label = "Play",
method = "animate",
args = list(
NULL,
list(frame = list(duration = 300, redraw = TRUE),
fromcurrent = TRUE)
)),
list(label = "Pause",
method = "animate",
args = list(NULL, list(frame = list(duration = 0))))
)
)
)
) %>% plotly_build()
# Update title in each frame --------------------------------------------
for (i in seq_len(ncol(eigfuncs))) {
p$x$frames[[i]]$layout <- list(
title = paste0(frame_name, ": ", eigvals[i])
)
}
return(p)
}
plot_3d_rectangle_onecol_vir <- function(mesh, fvals) {
colorscale = "Viridis"
opacity = 1
x <- mesh$loc[, 1]
y <- mesh$loc[, 2]
tri <- mesh$graph$tv
# fvals is a vector
x_range <- range(x)
y_range <- range(y)
z_range <- range(fvals)
# Static mesh3d plot
p <- plot_ly(
x = x, y = y, z = fvals,
i = tri[,1] - 1,
j = tri[,2] - 1,
k = tri[,3] - 1,
type = "mesh3d",
intensity = fvals,
colorscale = colorscale,
opacity = opacity,
flatshading = TRUE,
text = paste0(
"x: ", sprintf("%.3f", x), "<br>",
"y: ", sprintf("%.3f", y), "<br>",
"f: ", sprintf("%.5f", fvals)
),
hoverinfo = "text"
) %>% layout(
scene = global.scene.setter.rec.and.sphere(x_range, y_range, z_range)
)
return(p)
}
plot_3d_rectangle_onecol <- function(mesh, fvals) {
opacity <- 1
surface_color <- "blue"
x <- mesh$loc[, 1]
y <- mesh$loc[, 2]
tri <- mesh$graph$tv
x_range <- range(x)
y_range <- range(y)
z_range <- range(fvals)
# Static mesh3d plot
p <- plot_ly(
x = x, y = y, z = fvals,
i = tri[,1] - 1,
j = tri[,2] - 1,
k = tri[,3] - 1,
type = "mesh3d",
color = surface_color, # Fixed blue color
opacity = opacity,
flatshading = TRUE,
text = paste0(
"x: ", sprintf("%.3f", x), "<br>",
"y: ", sprintf("%.3f", y), "<br>",
"f: ", sprintf("%.5f", fvals)
),
hoverinfo = "text"
) %>% layout(
scene = global.scene.setter.rec.and.sphere(x_range, y_range, z_range)
)
return(p)
}
plot_3d_square_mesh <- function(mesh) {
x <- mesh$loc[,1]
y <- mesh$loc[,2]
# Flat square (or use mesh$loc[,3] if available)
if (ncol(mesh$loc) >= 3) {
z <- mesh$loc[,3]
} else {
z <- rep(0, length(x))
}
x_range <- range(x)
y_range <- range(y)
z_range <- c(0,1)
tri <- mesh$graph$tv
# ---- Uniform gray color for faces ----
gray_rgb <- rep(list(rgb(180, 180, 180, maxColorValue = 255)), length(x))
white_rgb <- rep(list(rgb(255, 255, 255, alpha = 255, maxColorValue = 255)), length(x))
df3 <- data.frame(x = x,
y = y,
z = z)
# ---- Plot mesh faces ----
p <- plot_ly() %>% add_trace(
x = x,
y = y,
z = z,
i = tri[,1] - 1,
j = tri[,2] - 1,
k = tri[,3] - 1,
type = "mesh3d",
vertexcolor = white_rgb,
flatshading = TRUE,
opacity = 1,
showscale = FALSE
)
# ---- Extract unique edges ----
edges <- rbind(
tri[,1:2],
tri[,2:3],
tri[,c(3,1)]
)
edges <- t(apply(edges, 1, sort))
edges <- unique(edges)
# ---- Efficient edges as one trace using NA breaks ----
x_edges <- as.vector(t(cbind(x[edges[,1]], x[edges[,2]], NA)))
y_edges <- as.vector(t(cbind(y[edges[,1]], y[edges[,2]], NA)))
z_edges <- as.vector(t(cbind(z[edges[,1]], z[edges[,2]], NA)))
p <- add_trace(
p,
x = x_edges, y = y_edges, z = z_edges,
type = "scatter3d",
mode = "lines",
line = list(color = "black", width = 3),
showlegend = FALSE,
hoverinfo = "none"
) %>%
add_trace(data = df3, x = ~x, y = ~y, z = ~z, mode = "markers", type = "scatter3d",
marker = list(size = 4, color = "gray", symbol = 104)) %>%
layout(
scene = global.scene.setter.rec.and.sphere(x_range, y_range, z_range)
)
return(p)
}
plot_3d_square_mesh_with_hat <- function(mesh,
hat_nodes = NULL,
hat_height = 1,
hat_color = NULL,
hat_alpha = 0.6) {
x <- mesh$loc[,1]
y <- mesh$loc[,2]
n <- length(x)
z0 <- rep(0, n)
tri <- mesh$graph$tv
x_range <- range(x)
y_range <- range(y)
z_range <- c(0, hat_height)
# ---- Extract unique edges once ----
edges <- rbind(tri[,1:2], tri[,2:3], tri[,c(3,1)])
edges <- t(apply(edges, 1, sort))
edges <- unique(edges)
# Base mesh color
white_rgb <- rep(list(rgb(255, 255, 255, alpha = 255, maxColorValue = 255)), n)
df3 <- data.frame(x = x,
y = y,
z = rep(0, length(x)))
# ---- Base mesh ----
p <- plot_ly() %>% add_trace(
x = x, y = y, z = z0,
i = tri[,1] - 1, j = tri[,2] - 1, k = tri[,3] - 1,
type = "mesh3d",
vertexcolor = white_rgb,
flatshading = TRUE,
opacity = 1,
name = "Mesh"
)
# ---- Base mesh edges ----
x_edges <- as.vector(t(cbind(x[edges[,1]], x[edges[,2]], NA)))
y_edges <- as.vector(t(cbind(y[edges[,1]], y[edges[,2]], NA)))
z_edges <- as.vector(t(cbind(z0[edges[,1]], z0[edges[,2]], NA)))
# ---- Hat functions ----
if (!is.null(hat_nodes)) {
if (is.null(hat_color)) {
hat_color <- grDevices::rainbow(length(hat_nodes))
}
if (length(hat_color) != length(hat_nodes)) {
stop("length(hat_color) must equal length(hat_nodes)")
}
for (k in seq_along(hat_nodes)) {
i0 <- hat_nodes[k]
phi <- rep(0, n)
phi[i0] <- 1
z_hat <- hat_height * phi
hat_rgb <- rep(list(adjustcolor(hat_color[k], alpha.f = hat_alpha)), n)
# ---- Hat surface ----
p <- add_trace(p,
x = x, y = y, z = z_hat,
i = tri[,1] - 1, j = tri[,2] - 1, k = tri[,3] - 1,
type = "mesh3d",
vertexcolor = hat_rgb,
flatshading = TRUE,
opacity = 1,
name = paste0("phi_", i0),
showlegend = FALSE
)
# ---- Hat edges (wireframe) ----
x_hat_edges <- as.vector(t(cbind(x[edges[,1]], x[edges[,2]], NA)))
y_hat_edges <- as.vector(t(cbind(y[edges[,1]], y[edges[,2]], NA)))
z_hat_edges <- as.vector(t(cbind(z_hat[edges[,1]], z_hat[edges[,2]], NA)))
p <- add_trace(p,
x = x_hat_edges, y = y_hat_edges, z = z_hat_edges,
type = "scatter3d", mode = "lines",
line = list(color = hat_color[k], width = 2),
showlegend = FALSE
)
# # Mark hat center node
# p <- add_trace(p,
# x = x[i0], y = y[i0], z = z_hat[i0],
# type = "scatter3d", mode = "markers",
# marker = list(size = 8, color = hat_color[k]),
# name = paste0("node_", i0)
# )
}
}
p <- add_trace(p,
x = x_edges, y = y_edges, z = z_edges,
type = "scatter3d", mode = "lines",
line = list(color = "black", width = 4),
showlegend = FALSE
)
p <- p %>%
add_trace(data = df3, x = ~x, y = ~y, z = ~z, mode = "markers", type = "scatter3d",
marker = list(size = 4, color = "gray", symbol = 104), showlegend = FALSE) %>%
layout(
scene = global.scene.setter.rec.and.sphere.for.hat(x_range, y_range, z_range)
)
return(p)
}plot_3d_sphere_scatter()plot_3d_sphere_scatter <- function(mesh,
eigvals,
eigfuncs,
fixed_colorscale = TRUE) {
colorscale = "Viridis"
x <- mesh$loc[,1]
y <- mesh$loc[,2]
z <- mesh$loc[,3]
# Compute global color limits if fixed
if (fixed_colorscale) {
cmin <- min(eigfuncs)
cmax <- max(eigfuncs)
} else {
cmin <- NULL
cmax <- NULL
}
# Create frames
frames <- lapply(seq_len(ncol(eigfuncs)), function(i) {
fvals <- eigfuncs[,i]
list(
name = as.character(i),
data = list(list(
x = x,
y = y,
z = z,
type = "scatter3d",
mode = "markers",
marker = list(
size = 5,
color = fvals,
colorscale = colorscale,
showscale = TRUE,
cmin = cmin,
cmax = cmax
),
text = paste0(
"x: ", sprintf("%.3f", x), "<br>",
"y: ", sprintf("%.3f", y), "<br>",
"z: ", sprintf("%.3f", z), "<br>",
"f: ", sprintf("%.5f", fvals)
),
hoverinfo = "text"
))
)
})
# Initial plot (frame 1)
fvals0 <- eigfuncs[,1]
p <- plot_ly(
x = x, y = y, z = z,
type = "scatter3d",
mode = "markers",
marker = list(
size = 5,
color = fvals0,
colorscale = colorscale,
showscale = TRUE,
cmin = cmin,
cmax = cmax
),
text = paste0(
"x: ", sprintf("%.3f", x), "<br>",
"y: ", sprintf("%.3f", y), "<br>",
"z: ", sprintf("%.3f", z), "<br>",
"f: ", sprintf("%.5f", fvals0)
),
hoverinfo = "text",
frame = "1"
)
frame_name <- deparse(substitute(eigvals))
p$x$frames <- frames
x_range <- range(x)
y_range <- range(y)
z_range <- range(z)
# Slider + play/pause
p <- p %>% layout(
title = paste0(frame_name, ": ", eigvals[1]),
scene = global.scene.setter.rec.and.sphere(x_range, y_range, z_range),
sliders = list(
list(
active = 0,
currentvalue = list(prefix = "Mode: "),
pad = list(t = 50),
steps = lapply(seq_len(ncol(eigfuncs)), function(i) {
list(
label = as.character(i),
method = "animate",
args = list(list(as.character(i)),
list(mode = "immediate",
frame = list(duration = 300, redraw = TRUE),
transition = list(duration = 0)))
)
})
)
),
updatemenus = list(
list(
type = "buttons",
showactive = FALSE,
y = 1,
x = 1.15,
xanchor = "right",
yanchor = "top",
buttons = list(
list(
label = "Play",
method = "animate",
args = list(NULL,
list(frame = list(duration = 300, redraw = TRUE),
fromcurrent = TRUE,
mode = "immediate"))
),
list(
label = "Pause",
method = "animate",
args = list(NULL,
list(frame = list(duration = 0, redraw = FALSE),
mode = "immediate"))
)
)
)
)
) %>% plotly_build()
# Update title per frame
for (i in seq_len(ncol(eigfuncs))) {
p$x$frames[[i]]$layout <- list(
title = paste0(frame_name, ": ", eigvals[i])
)
}
return(p)
}
plot_3d_slider_sphere <- function(mesh, eigvals, eigfuncs, fixed_colorscale = TRUE) {
colorscale = "Viridis"
opacity = 1
x <- mesh$loc[, 1]
y <- mesh$loc[, 2]
z <- mesh$loc[, 3]
tri <- mesh$graph$tv
# Global intensity limits if fixed
if (fixed_colorscale) {
cmin <- min(eigfuncs)
cmax <- max(eigfuncs)
} else {
cmin <- NULL
cmax <- NULL
}
# Create frames
frames <- lapply(seq_len(ncol(eigfuncs)), function(i) {
fvals <- eigfuncs[,i]
list(
name = as.character(i),
data = list(list(
x = x,
y = y,
z = z,
i = tri[,1] - 1,
j = tri[,2] - 1,
k = tri[,3] - 1,
type = "mesh3d",
intensity = fvals,
colorscale = colorscale,
opacity = opacity,
flatshading = TRUE,
cmin = cmin,
cmax = cmax,
text = paste0(
"x: ", sprintf("%.3f", x), "<br>",
"y: ", sprintf("%.3f", y), "<br>",
"z: ", sprintf("%.3f", z), "<br>",
"f: ", sprintf("%.5f", fvals)
),
hoverinfo = "text"
))
)
})
# Initial plot
fvals0 <- eigfuncs[,1]
p <- plot_ly(
x = x, y = y, z = z,
i = tri[,1] - 1,
j = tri[,2] - 1,
k = tri[,3] - 1,
type = "mesh3d",
intensity = fvals0,
colorscale = colorscale,
opacity = opacity,
flatshading = TRUE,
cmin = cmin,
cmax = cmax,
text = paste0(
"x: ", sprintf("%.3f", x), "<br>",
"y: ", sprintf("%.3f", y), "<br>",
"z: ", sprintf("%.3f", z), "<br>",
"f: ", sprintf("%.5f", fvals0)
),
colorbar = list(
thickness = 15,
len = 0.5, # fraction of the height of the plot
y = 0.5, # center vertically (0 = bottom, 1 = top)
yanchor = "middle",
title = ""
),
hoverinfo = "text",
frame = "1"
)
frame_name <- deparse(substitute(eigvals))
p$x$frames <- frames
# Layout + slider + buttons
p <- p %>% layout(
title = paste0(frame_name, ": ", eigvals[1]),
scene = global.scene.setter.rec.and.sphere(range(x), range(y), range(z)),
sliders = list(
list(
active = 0,
currentvalue = list(prefix = "Frame: "),
pad = list(t = 50),
steps = lapply(seq_len(ncol(eigfuncs)), function(i) {
list(
label = as.character(i),
method = "animate",
args = list(list(as.character(i)),
list(mode = "immediate",
frame = list(duration = 300, redraw = TRUE),
transition = list(duration = 0)))
)
})
)
),
updatemenus = list(
list(
type = "buttons",
showactive = FALSE,
y = 1,
x = 1.15,
xanchor = "right",
yanchor = "top",
buttons = list(
list(label = "Play",
method = "animate",
args = list(NULL, list(frame = list(duration = 300, redraw = TRUE),
fromcurrent = TRUE, mode = "immediate"))),
list(label = "Pause",
method = "animate",
args = list(NULL, list(frame = list(duration = 0, redraw = FALSE),
mode = "immediate")))
)
)
)
) %>% plotly_build()
# Add titles for frames
for (i in seq_len(ncol(eigfuncs))) {
p$x$frames[[i]]$layout <- list(
title = paste0(frame_name, ": ", eigvals[i])
)
}
return(p)
}
plot_3d_sphere_scatter_onecol <- function(mesh, eigfuncs) {
colorscale = "Viridis"
x <- mesh$loc[,1]
y <- mesh$loc[,2]
z <- mesh$loc[,3]
# eigfuncs is a vector
fvals <- eigfuncs
x_range <- range(x)
y_range <- range(y)
z_range <- range(z)
# Build static 3D scatter plot
p <- plot_ly(
x = x, y = y, z = z,
type = "scatter3d",
mode = "markers",
marker = list(
size = 5,
color = fvals,
colorscale = colorscale,
showscale = TRUE
),
text = paste0(
"x: ", sprintf("%.3f", x), "<br>",
"y: ", sprintf("%.3f", y), "<br>",
"z: ", sprintf("%.3f", z), "<br>",
"f: ", sprintf("%.5f", fvals)
),
hoverinfo = "text"
) %>% layout(
scene = global.scene.setter.rec.and.sphere(x_range, y_range, z_range)
)
return(p)
}
# plot_3d_sphere_surface_onecol <- function(mesh, fvals) {
#
# colorscale = "Viridis"
# opacity = 1
#
#
# x <- mesh$loc[, 1]
# y <- mesh$loc[, 2]
# z <- mesh$loc[, 3]
# tri <- mesh$graph$tv
#
# # fvals is a vector
#
# x_range <- range(x)
# y_range <- range(y)
# z_range <- range(z)
#
# # Static mesh3d plot
# p <- plot_ly(
# x = x, y = y, z = z,
# i = tri[,1] - 1,
# j = tri[,2] - 1,
# k = tri[,3] - 1,
# type = "mesh3d",
# intensity = fvals,
# colorscale = colorscale,
# opacity = opacity,
# flatshading = TRUE,
# text = paste0(
# "x: ", sprintf("%.3f", x), "<br>",
# "y: ", sprintf("%.3f", y), "<br>",
# "z: ", sprintf("%.3f", z), "<br>",
# "f: ", sprintf("%.5f", fvals)
# ),
# hoverinfo = "text"
# ) %>% layout(
# scene = global.scene.setter.rec.and.sphere(x_range, y_range, z_range)
# )
#
# return(p)
# }
plot_3d_sphere_surface_onecol <- function(mesh, fvals) {
colorscale = "Viridis"
opacity = 1
x <- mesh$loc[, 1]
y <- mesh$loc[, 2]
z <- mesh$loc[, 3]
tri <- mesh$graph$tv
x_range <- range(x)
y_range <- range(y)
z_range <- range(z)
# Static mesh3d plot
p <- plot_ly(
x = x, y = y, z = z,
i = tri[,1] - 1,
j = tri[,2] - 1,
k = tri[,3] - 1,
type = "mesh3d",
intensity = fvals,
colorscale = colorscale,
opacity = opacity,
flatshading = TRUE,
text = paste0(
"x: ", sprintf("%.3f", x), "<br>",
"y: ", sprintf("%.3f", y), "<br>",
"z: ", sprintf("%.3f", z), "<br>",
"f: ", sprintf("%.5f", fvals)
),
hoverinfo = "text",
colorbar = list(
thickness = 15,
len = 0.5, # fraction of the height of the plot
y = 0.5, # center vertically (0 = bottom, 1 = top)
yanchor = "middle",
title = ""
)
) %>% layout(
scene = global.scene.setter.rec.and.sphere(x_range, y_range, z_range)
)
return(p)
}
plot_3d_mesh_edges_faces_old <- function(mesh) {
x <- mesh$loc[,1]
y <- mesh$loc[,2]
z <- mesh$loc[,3]
tri <- mesh$graph$tv
# ---- Uniform gray color for faces (vertexcolor) ----
gray_rgb <- rep(list(rgb(180, 180, 180, maxColorValue = 255)), length(x))
p <- plot_ly() %>% add_trace(
x = x,
y = y,
z = z,
i = tri[,1] - 1,
j = tri[,2] - 1,
k = tri[,3] - 1,
type = "mesh3d",
vertexcolor = gray_rgb, # <- the correct way
flatshading = TRUE,
opacity = 1,
showscale = FALSE
)
# ---- Extract unique edges ----
edges <- rbind(
tri[,1:2],
tri[,2:3],
tri[,c(3,1)]
)
edges <- t(apply(edges, 1, sort))
edges <- unique(edges)
# ---- Add edges as blue lines ----
for (e in 1:nrow(edges)) {
i1 <- edges[e,1]
i2 <- edges[e,2]
p <- p %>% add_trace(
x = c(x[i1], x[i2]),
y = c(y[i1], y[i2]),
z = c(z[i1], z[i2]),
type = "scatter3d",
mode = "lines",
line = list(color = "blue", width = 3),
hoverinfo = "none",
showlegend = FALSE
)
}
return(p)
}
plot_3d_mesh_edges_faces <- function(mesh) {
x <- mesh$loc[,1]
y <- mesh$loc[,2]
z <- mesh$loc[,3]
tri <- mesh$graph$tv
x_range <- range(x)
y_range <- range(y)
z_range <- range(z)
# ---- Uniform gray color for faces ----
gray_rgb <- rep(list(rgb(180, 180, 180, maxColorValue = 255)), length(x))
white_rgb <- rep(list(rgb(255, 255, 255, alpha = 255, maxColorValue = 255)), length(x))
df3 <- data.frame(x = x,
y = y,
z = z)
# ---- Plot mesh faces ----
p <- plot_ly() %>% add_trace(
x = x,
y = y,
z = z,
i = tri[,1] - 1,
j = tri[,2] - 1,
k = tri[,3] - 1,
type = "mesh3d",
vertexcolor = white_rgb,
flatshading = TRUE,
opacity = 1,
showscale = FALSE
)
# ---- Extract unique edges ----
edges <- rbind(
tri[,1:2],
tri[,2:3],
tri[,c(3,1)]
)
edges <- t(apply(edges, 1, sort))
edges <- unique(edges)
# ---- Efficient edges as a single trace using NA breaks ----
x_edges <- as.vector(t(cbind(x[edges[,1]], x[edges[,2]], NA)))
y_edges <- as.vector(t(cbind(y[edges[,1]], y[edges[,2]], NA)))
z_edges <- as.vector(t(cbind(z[edges[,1]], z[edges[,2]], NA)))
p <- add_trace(
p,
x = x_edges, y = y_edges, z = z_edges,
type = "scatter3d",
mode = "lines",
line = list(color = "black", width = 3),
showlegend = FALSE,
hoverinfo = "none"
) %>%
add_trace(data = df3, x = ~x, y = ~y, z = ~z, mode = "markers", type = "scatter3d",
marker = list(size = 4, color = "gray", symbol = 104)) %>%
layout(
scene = global.scene.setter.rec.and.sphere(x_range, y_range, z_range)
)
return(p)
}
plot_3d_mesh_edges_faces_with_hat <- function(mesh,
hat_nodes = NULL,
hat_height = 0.1,
hat_color = NULL,
hat_alpha = 0.6) {
x <- mesh$loc[,1]
y <- mesh$loc[,2]
z <- mesh$loc[,3]
n <- length(x)
tri <- mesh$graph$tv
x_range <- range(x)*2
y_range <- range(y)*2
z_range <- range(z)*2
# ---- Extract unique edges once ----
edges <- rbind(tri[,1:2], tri[,2:3], tri[,c(3,1)])
edges <- t(apply(edges, 1, sort))
edges <- unique(edges)
# Base mesh colors
white_rgb <- rep(list(rgb(255, 255, 255, alpha = 255, maxColorValue = 255)), n)
df3 <- data.frame(x = x, y = y, z = z)
# ---- Base spherical mesh ----
p <- plot_ly() %>% add_trace(
x = x, y = y, z = z,
i = tri[,1] - 1, j = tri[,2] - 1, k = tri[,3] - 1,
type = "mesh3d",
vertexcolor = white_rgb,
flatshading = TRUE,
opacity = 1,
showscale = FALSE,
name = "Sphere mesh"
)
# ---- Base mesh edges (single efficient trace) ----
x_edges <- as.vector(t(cbind(x[edges[,1]], x[edges[,2]], NA)))
y_edges <- as.vector(t(cbind(y[edges[,1]], y[edges[,2]], NA)))
z_edges <- as.vector(t(cbind(z[edges[,1]], z[edges[,2]], NA)))
# ---- Hat functions on sphere ----
if (!is.null(hat_nodes)) {
# Validate indices
if (any(hat_nodes < 1 | hat_nodes > n)) {
stop("hat_nodes contains invalid node indices")
}
# Color recycling
if (is.null(hat_color)) {
hat_color <- grDevices::rainbow(length(hat_nodes))
}
hat_color <- rep(hat_color, length.out = length(hat_nodes))
for (k in seq_along(hat_nodes)) {
i0 <- hat_nodes[k]
# Nodal hat basis
phi <- rep(0, n)
phi[i0] <- 1
# Lift radially outward (important for sphere!)
r <- sqrt(x^2 + y^2 + z^2)
nx <- x / r
ny <- y / r
nz <- z / r
x_hat <- x + hat_height * phi * nx
y_hat <- y + hat_height * phi * ny
z_hat <- z + hat_height * phi * nz
# Hat surface color
hat_rgb <- rep(list(adjustcolor(hat_color[k], alpha.f = hat_alpha)), n)
# ---- Hat surface ----
p <- add_trace(p,
x = x_hat, y = y_hat, z = z_hat,
i = tri[,1] - 1, j = tri[,2] - 1, k = tri[,3] - 1,
type = "mesh3d",
vertexcolor = hat_rgb,
flatshading = TRUE,
opacity = 1,
name = paste0("phi_", i0),
showlegend = FALSE
)
# ---- Hat wireframe edges ----
x_hat_edges <- as.vector(t(cbind(x_hat[edges[,1]], x_hat[edges[,2]], NA)))
y_hat_edges <- as.vector(t(cbind(y_hat[edges[,1]], y_hat[edges[,2]], NA)))
z_hat_edges <- as.vector(t(cbind(z_hat[edges[,1]], z_hat[edges[,2]], NA)))
p <- add_trace(p,
x = x_hat_edges, y = y_hat_edges, z = z_hat_edges,
type = "scatter3d", mode = "lines",
line = list(color = hat_color[k], width = 2),
showlegend = FALSE,
hoverinfo = "none"
)
}
}
p <- add_trace(p,
x = x_edges, y = y_edges, z = z_edges,
type = "scatter3d", mode = "lines",
line = list(color = "black", width = 3),
showlegend = FALSE,
hoverinfo = "none"
)
# ---- Nodes ----
p <- add_trace(p,
data = df3,
x = ~x, y = ~y, z = ~z,
mode = "markers",
type = "scatter3d",
marker = list(size = 4, color = "gray", symbol = 104),
showlegend = FALSE
)
# ---- Scene ----
p <- p %>% layout(
scene = global.scene.setter.rec.and.sphere(x_range, y_range, z_range)
)
return(p)
}graph.plotter.3d.onecol <- function(graph, vec) {
# Coordinates on the mesh
x <- plotting.order(graph$mesh$V[, 1], graph)
y <- plotting.order(graph$mesh$V[, 2], graph)
z <- plotting.order(vec, graph) # vector to plot
# Axis ranges
x_range <- range(x)
y_range <- range(y)
z_range <- range(z)
z_range[1] <- z_range[1] - 10^-4
# Vertical lines (vectorized)
n <- length(x)
x_vert <- rep(x, each = 3)
y_vert <- rep(y, each = 3)
z_vert <- as.vector(t(cbind(0, z, NA)))
# Create plot
p <- plot_ly() %>%
# Main 3D curve over the graph
add_trace(
x = x, y = y, z = z,
type = "scatter3d", mode = "lines",
line = list(color = "blue", width = 3),
showlegend = FALSE
) %>%
# Graph base
add_trace(
x = x, y = y, z = z*0,
type = "scatter3d", mode = "lines",
line = list(color = "black", width = 3),
showlegend = FALSE
) %>%
# Vertical lines from base to curve
add_trace(
x = x_vert, y = y_vert, z = z_vert,
type = "scatter3d", mode = "lines",
line = list(color = "gray", width = 0.5),
showlegend = FALSE
) %>%
layout(
scene = list(xaxis = list(title = "x", range = x_range),
yaxis = list(title = "y", range = y_range),
zaxis = list(
title = "z",
range = z_range,
exponentformat = "e",
tickformat = NULL#".4f" # prevents SI prefixes like μ, m, k
),
aspectratio = list(x = 2*(1+2/pi),
y = 2*(2/pi),
z = 1*(2/pi)),
camera = list(eye = list(x = 5,
y = 3,
z = 3.5),
center = list(x = (1+2/pi)/2,
y = 0,
z = 0)))
)
return(p)
}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).