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)
}capture.output(
knitr::purl(here::here("functionality.Rmd"), output = here::here("functionality.R")),
file = here::here("purl_log.txt")
)
source(here::here("functionality.R"))Let each edge \(e\in\mathcal{E}\) be subdivided into \(n_{e}\geq 2\) regular segments of length \(h_{e}\), and be delimited by the nodes \(0 = x_0^{e},x_1^{e},\dots,x_{n_{e}-1}^{e}, x_{n_{e}}^{e} = \ell_{e}\). For each \(j = 1,\dots,n_{e}-1\), we consider the following standard hat basis functions \[\begin{equation*} \varphi_j^{e}(x)=\begin{cases} 1-\dfrac{|x_j^{e}-x|}{h_{e}},&\text{ if }x_{j-1}^{e}\leq x\leq x_{j+1}^{e},\\ 0,&\text{ otherwise}. \end{cases} \end{equation*}\] For each \(e\in\mathcal{E}\), the set of hat functions \(\left\{\varphi_1^{e},\dots,\varphi_{n_{e}-1}^{e}\right\}\) is a basis for the space \[\begin{equation*} V_{h_{e}} = \left\{w\in H_0^1(e)\;\Big|\;\forall j = 0,1,\dots,n_{e}-1:w|_{[x_j^{e}, x_{j+1}^{e}]}\in\mathbb{P}^1\right\}, \end{equation*}\] where \(\mathbb{P}^1\) is the space of linear functions on \([0,\ell_{e}]\). For each vertex \(v\in\mathcal{V}\), we define \[\begin{equation*} \mathcal{N}_v = \left\{\bigcup_{e\in\left\{e\in\mathcal{E}_v: v = x_0^e\right\}}[v,x_1^e]\right\}\bigcup\left\{\bigcup_{e\in\left\{e\in\mathcal{E}_v: v = x^e_{n_e}\right\}}[x^e_{n_e-1},v]\right\}, \end{equation*}\] which is a star-shaped set with center at \(v\) and rays made of the segments contiguous to \(v\). On \(\mathcal{N}_v\), we define the hat functions as \[\begin{equation*} \phi_v(x)=\begin{cases} 1-\dfrac{|x_v^{e}-x|}{h_{e}},&\text{ if }x\in\mathcal{N}_v\cap e \text{ and }e\in\mathcal{E}_v,\\ 0,&\text{ otherwise}, \end{cases} \end{equation*}\] where \(x_v^e\) is either \(x_0^e\) or \(x_{n_e}^e\) depending on the edge direction and its parameterization. See Arioli and Benzi (2018) for more details. Figure 1 shows an illustration of the basis function system \(\{\varphi_j^e, \phi_v\}\) on the tadpole graph (in black). Standard hat functions associated with internal edge nodes are shown in blue, while special vertex-centered functions are highlighted in red.
For an additional illustration, go to the Basis page.
graph <- gets.graph.tadpole(h = 1/4)
graph_to_get_loc <- gets.graph.tadpole(h = 1/40)
loc <- graph_to_get_loc$get_mesh_locations()
A <- as.matrix(graph$fem_basis(loc))
A <- apply(A, 2, function(x) plotting.order(x, graph_to_get_loc))
A <- rbind(A, rep(NA, ncol(A))) # Add a row of NAs for the plotting
x <- graph_to_get_loc$mesh$V[, 1]
y <- graph_to_get_loc$mesh$V[, 2]
x <- c(plotting.order(x, graph_to_get_loc), NA)
y <- c(plotting.order(y, graph_to_get_loc), NA)
x_range <- range(x, na.rm = TRUE)
y_range <- range(y, na.rm = TRUE)
z_range <- c(0,1)
df3 <- data.frame(x = graph$mesh$V[, 1],
y = graph$mesh$V[, 2],
z = rep(0, length(graph$mesh$V[, 1])))
# Start plot
p_b1 <- plot_ly() |>
add_trace(x = rep(x, times = graph$nV),
y = rep(y, times = graph$nV),
z = as.vector(A[, 1:graph$nV]),
type = "scatter3d",
mode = "lines",
line = list(color = "red", width = 3),
showlegend = FALSE) |>
add_trace(x = rep(x, times = ncol(A) - graph$nV),
y = rep(y, times = ncol(A) - graph$nV),
z = as.vector(A[, (graph$nV+1):ncol(A)]),
type = "scatter3d",
mode = "lines",
line = list(color = "blue", width = 3),
showlegend = FALSE) |>
add_trace(x = rep(x, each = 3),
y = rep(y, each = 3),
z = unlist(lapply(apply(A, 1, max, na.rm = TRUE), function(zj) c(0, zj, NA))),
type = "scatter3d",
mode = "lines",
line = list(color = "gray", width = 0.5),
showlegend = FALSE) |>
add_trace(x = x,
y = y,
z = x*0,
type = "scatter3d",
mode = "lines",
line = list(color = "black", width = 3),
showlegend = FALSE) |>
add_trace(x = rep(x, times = graph$nV),
y = rep(y, times = graph$nV),
z = c(replace(rep(NA, nrow(A)), 1:11, 0),
replace(rep(NA, nrow(A)), c(31:51, 111:121), 0)),
type = "scatter3d",
mode = "lines",
line = list(color = "green", width = 3),
showlegend = FALSE) |>
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(x_range, y_range, z_range, z_aspectratio = 1))
layout(scene = 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 = 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))))
p_b1Figure 1: Illustration of the system of basis functions \(\{\varphi_j^e, \phi_v\}\) on the tadpole graph. Notice that the sets \(\mathcal{N}_{v}\) are depicted in green and their corresponding basis functions are shown in red.
Having introduced the system of basis functions \(\{\varphi_j^e, \phi_v\}\), which will be referred to as \(\{\psi_j\}_{j=1}^{N_h}\), we can now define the finite element space \(V_h\subset H^1(\Gamma)\) as \(V_h = \left(\bigoplus_{e\in\mathcal{E}} V_{h_e}\right)\bigoplus V_v\), where \(V_v = \text{span}\left(\{\phi_v:v\in\mathcal{V}\}\right)\) and \(\dim\left(V_h\right)\) is given by \(N_h = |\mathcal{V}| + \sum_{e\in\mathcal{E}}n_e\).
Let \(\Gamma_T = (\mathcal{V},\mathcal{E})\) characterize the tadpole graph with \(\mathcal{V}= \{v_1,v_2\}\) and \(\mathcal{E}= \{e_1,e_2\}\). The left edge \(e_1\) has length 1 and the circular edge \(e_2\) has length 2. A point on \(e_1\) is parameterized via \(s=\left(e_1, t\right)\) for \(t \in[0,1]\) and a point on \(e_2\) via \(s=\left(e_2, t\right)\) for \(t\in[0,2]\). One can verify that \(-\Delta_\Gamma\) has eigenvalues \(0,\left\{(i \pi / 2)^2\right\}_{i \in \mathbb{N}}\) and \(\left\{(i \pi / 2)^2\right\}_{2 i \in \mathbb{N}}\) with corresponding eigenfunctions \(\phi_0\), \(\left\{\phi_i\right\}_{i \in \mathbb{N}}\), and \(\left\{\psi_i\right\}_{2 i \in \mathbb{N}}\) given by \(\phi_0(s)=1 / \sqrt{3}\) and \[\begin{equation*} \phi_i(s)=C_{\phi, i}\begin{cases} -2 \sin \left(\dfrac{i\pi}{2}\right) \cos \left(\dfrac{i \pi s}{2}\right), & s \in e_1, \\ \sin \left(\dfrac{i \pi s}{2}\right), & s \in e_2, \end{cases}, \quad \psi_i(s)=\dfrac{\sqrt{3}}{\sqrt{2}} \begin{cases} (-1)^{i / 2} \cos \left(\dfrac{i \pi s}{2}\right), & s \in e_1, \\ \cos \left(\dfrac{i \pi s}{2}\right), & s \in e_2, \end{cases}, \end{equation*}\] where \(C_{\phi, i}=1\) if \(i\) is even and \(C_{\phi, i}=1 / \sqrt{3}\) otherwise. Moreover, these functions form an orthonormal basis for \(L_2(\Gamma_T)\).
graph <- gets.graph.tadpole(h = 0.02)
eigenpairs <- gets.eigen.params(N_finite = 10, kappa = 1, alpha = 0.5, graph)
EIGENFUN <- eigenpairs$EIGENFUN
EIGENVAL <- eigenpairs$EIGENVAL
INDEX <- eigenpairs$INDEXFigure 2: Illustration of the eigenfunctions on the tadpole graph.
The piecewise constant projection \(U_{\text{piecewise}}(t^*_\ell)\) of the approximated values \(U_{\text{approx}}(t_k)\), defined on a coarse time grid \(\{t_k\}_{k=0}^{N}\), onto a finer grid \(\{t^*_\ell\}_{\ell=0}^{M}\), is given by
\[\begin{equation} U_{\text{piecewise}}(t^*_\ell) = \begin{cases} U_{\text{approx}}(t_0), & \text{if } t^*_\ell = 0 \\\\ U_{\text{approx}}(t_k), & \text{if } t^*_\ell \in (t_{k-1}, t_k], \quad \text{for } k = 1, \dots, N. \end{cases} \end{equation}\]
This defines a function that is constant on each interval \((t_{k-1}, t_k]\), and takes the value of the approximation at the right endpoint \(t_k\) of the interval.
See the Fuctionality
page for the implementation of the function
construct_piecewise_projection() that performs this
operation.
A function \(U_{\text{coarse}}(s)\) defined on a coarse mesh with \(N_{h}\) nodes can be projected onto a fine mesh with \(N_{h_{\text{ok}}}\) nodes by doing \(U_{\text{fine}}(s) = \boldsymbol{\Psi} U_{\text{coarse}}(s)\), where \(\boldsymbol{\Psi}\) is a matrix with entries \(\boldsymbol{\Psi}_{ij}=\psi^j_h(s_i)\) for \(j=1\dots, N_h\) and \(i = 1,\dots N_{h_{\text{ok}}}\).
# Parameters
T_final <- 2
kappa <- 15
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
AAA = 1
OMEGA = pi
# Overkill parameters
overkill_time_step <- 0.02
overkill_h <- 0.02
# Finest time and space mesh
overkill_time_seq <- seq(0, T_final, length.out = ((T_final - 0) / overkill_time_step + 1))
overkill_graph <- gets.graph.tadpole(h = overkill_h)
# Compute the weights on the finest mesh
overkill_graph$compute_fem() # This is needed to compute the weights
overkill_weights <- overkill_graph$mesh$weights
alpha <- 0.5 # from 0.5 to 2
beta <- alpha / 2
# Compute the eigenvalues and eigenfunctions on the finest mesh
overkill_eigen_params <- gets.eigen.params(N_finite = N_finite, kappa = kappa, alpha = alpha, graph = overkill_graph)
EIGENVAL_ALPHA <- overkill_eigen_params$EIGENVAL_ALPHA # Eigenvalues (they are independent of the meshes)
overkill_EIGENFUN <- overkill_eigen_params$EIGENFUN # Eigenfunctions on the finest mesh
# Compute the true solution on the finest mesh
overkill_U_true <- overkill_EIGENFUN %*%
outer(1:length(coeff_U_0),
1:length(overkill_time_seq),
function(i, j) (coeff_U_0[i] + coeff_FF[i] * G_sin(t = overkill_time_seq[j], A = AAA, lambda_j_alpha_half = EIGENVAL_ALPHA[i], omega = OMEGA)) * exp(-EIGENVAL_ALPHA[i] * overkill_time_seq[j]))
h <- 0.2
time_step <- 0.2
m <- 1
graph <- gets.graph.tadpole(h = h)
graph$compute_fem()
G <- graph$mesh$G
C <- graph$mesh$C
L <- kappa^2*C + G
eigen_params <- gets.eigen.params(N_finite = N_finite, kappa = kappa, alpha = alpha, graph = graph)
EIGENFUN <- eigen_params$EIGENFUN
U_0 <- EIGENFUN %*% coeff_U_0 # Compute U_0 on the current mesh
A <- graph$fem_basis(overkill_graph$get_mesh_locations())
time_seq <- seq(0, T_final, length.out = ((T_final - 0) / time_step + 1))
my_op_frac <- my.fractional.operators.frac(L, beta, C, scale.factor = kappa^2, m = m, time_step)
INT_BASIS_EIGEN <- t(overkill_EIGENFUN) %*% overkill_graph$mesh$C %*% A
# Compute matrix F with columns F^k
FF_approx <- t(INT_BASIS_EIGEN) %*%
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))
U_approx <- matrix(NA, nrow = nrow(C), ncol = length(time_seq))
U_approx[, 1] <- U_0
for (k in 1:(length(time_seq) - 1)) {
U_approx[, k + 1] <- as.matrix(my.solver.frac(my_op_frac, my_op_frac$C %*% U_approx[, k] + time_step * FF_approx[, k + 1]))
}
projected_U_approx <- A %*% U_approx
projected_U_piecewise <- construct_piecewise_projection(projected_U_approx, time_seq, overkill_time_seq)graph.plotter.3d(overkill_graph, overkill_time_seq, overkill_time_seq, list(overkill_U_true = overkill_U_true, projected_U_piecewise = projected_U_piecewise))
Figure 3: Illustration of the projection of the approximated solution
\(U_{\text{approx}}\) (called
projected_U_piecewise) onto a fine space-time mesh.
overkill_U_true is the true solution \(U_{\text{true}}\).
start_idx <- which.min(abs(overkill_time_seq - 0))
end_idx <- which.min(abs(overkill_time_seq - T_final))
idx <- seq(start_idx, end_idx, by = 10)
graph.plotter.3d.two.meshes.time(overkill_graph, graph, time_seq, time_seq, fs_finer = list(projected_U_approx = projected_U_approx, overkill_U_true = overkill_U_true[,idx]), fs_coarser = list(U_approx = U_approx))
Figure 4: U_approx is the coarse approximated solution on
the coarse mesh, projected_U_approx is the coarse solution
on the fine mesh, and overkill_U_true is the true solution
on the fine mesh. We use the coarse time sequence.
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).