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.
# 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"))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.
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} \boldsymbol{U}^{k+1} = \boldsymbol{C}^{-1}\left(\sum_{k=1}^{m+1} a_k\left( \dfrac{\boldsymbol{L}\boldsymbol{C}^{-1}}{\kappa^2}-p_k\boldsymbol{I}\right)^{-1} + r\boldsymbol{I}\right) (\boldsymbol{C}\boldsymbol{U}^k+\tau \boldsymbol{F}^{k+1}) \end{equation}\] Let \(\boldsymbol{F}\) be the matrix with columns \(\boldsymbol{F}^k\). This matrix has entries \(\boldsymbol{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\boldsymbol{C}^{\text{ok}}\boldsymbol{e}_j\), where \(\boldsymbol{e}_j=[e_j(s_1)\dots e_j(s_{N_{h_{\text{ok}}}})]^\top\) and \(\boldsymbol{\psi}_i=[\psi^i_h(s_1)\dots \psi^i_h(s_{N_{h_{\text{ok}}}})]^\top\) are vectors of function evaluations on a fine spatial mesh with mesh size \(h_{\text{ok}}\) and nodes \(\{s_\ell\}_{\ell=1}^{N_{h_{\text{ok}}}}\). Matrix \(\boldsymbol{C}^{\text{ok}}\) contains the corresponding quadrature weights, with entries \(\boldsymbol{C}^{\text{ok}}_{nm} = (\psi_{h_\text{ok}}^n,\psi_{h_\text{ok}}^m)_{L_2(\Gamma)}\) for \(n,m = 1,\dots, N_{h_\text{ok}}\). 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_\text{ok}}\}_{\ell=1}^{N_{h_\text{ok}}}\) and mesh size \(h_\text{ok}\). If \(\boldsymbol{N}\) denotes the matrix with entries \(\boldsymbol{N}_{jk} = y_j g(t_k)\), then \(\boldsymbol{F}\) can be approximated as \([\boldsymbol{\psi}_1\dots \boldsymbol{\psi}_{N_h}]^\top \boldsymbol{C}^{\text{ok}} [\boldsymbol{e}_1\dots \boldsymbol{e}_{N_f}] \boldsymbol{N}\). In \(\texttt{R}\), \(\boldsymbol{F}\) and \(\boldsymbol{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()))}\\ \boldsymbol{C}^{\text{ok}} \iff&\quad \texttt{overkill_graph\$mesh\$C}\\ [\boldsymbol{e}_1\dots \boldsymbol{e}_{N_f}] \iff&\quad \texttt{overkill_EIGENFUN}\\ \boldsymbol{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))}\\ \boldsymbol{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))}\\ \boldsymbol{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.001
h <- 0.001
kappa <- 15
alpha <- 0.5
m = 1
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)
error <- sqrt(as.double(t(graph$mesh$weights) %*% (U_true - U_approx)^2 %*% rep(time_step, ncol(U_true))))The following is the error reported in the paper for the considered experiment.
## [1] "Total error = 0.004003178"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_approxstart_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)
graph.plotter.3d.single(aux_graph, FF_true[, idx], time_seq[idx])Figure 1: Time evolution of the right-hand side function \(f\).
abs_error <- abs(U_true[, idx] - U_approx[, idx])
graph.plotter.3d.single(aux_graph, abs_error, time_seq[idx])