Example on the sphere

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"))

library(fmesher)
library(Matrix)

We consider the following equation on the sphere

\[\begin{equation} \label{eq:maineq} \tag{1} \left\{ \begin{aligned} \partial_t u(s,t) + (\kappa^2 - \Delta_{\mathbb{S}^2})^{\alpha/2} u(s,t) &= f(s,t), && \quad (s,t) \in \mathbb{S}^2 \times (0, T), \\ u(s,0) &= u_0(s), && \quad s \in \mathbb{S}^2, \end{aligned} \right. \end{equation}\]

where \(\Delta_{\mathbb{S}^2}\) is the Laplace-Beltrami operator on the sphere \(\mathbb{S}^2\).

The eigenvalues and eigenfunctions of the Laplace-Beltrami operator on the sphere are known explicitly. The eigenvalues are given by \(\lambda_\ell = \ell(\ell+1)\), with multiplicity \(2\ell + 1\), and the eigenfunctions are the spherical harmonics \(Y_{\ell,m}\), for \(\ell = 0, 1, 2, \ldots\) and \(m = -\ell, \ldots, \ell\), which satisfy the eigenvalue equation

\[\begin{equation} -\Delta_{\mathbb{S}^2} Y_{\ell,m} = \lambda_\ell Y_{\ell,m}. \end{equation}\]

Below we use the eigenpairs to construct an exact solution to equation \(\eqref{eq:maineq}\) and compare it with the numerical approximation. We take \(u_0(s) = x_3 Y_{3,0}(s)\) with \(x_3 = 1\) and \(f(s,t) = y_2 \sin(\pi t) Y_{2,0}(s)\) with \(y_2 = 10\). Then the exact solution is given by \[\begin{equation} u(s,t) = x_3\text{e}^{-\lambda^{\alpha/2}_3t}Y_{3,0}(s)+y_2 G_2(t)\text{e}^{-\lambda^{\alpha/2}_2t}Y_{2,0}(s), \end{equation}\] where \(G_2(t)= \displaystyle\int_0^t \text{e}^{\lambda^{\alpha/2}_2r}\sin(\pi r)dr\). Below the exact solution is computed as U_true and compared with the numerical approximation U_approx.

We take parameters \(T = 2\), \(\kappa = 4\), \(\alpha = 1.8\), \(m=4\), \(h = 0.05\), and time step size \(\tau = 0.05\).

T_final <- 2
kappa <- 4


L_max <- 10
rot.inv <- TRUE


coeff_U_0 <- rep(0, L_max+1)
coeff_U_0[3] <- 1
coeff_FF <- rep(0, L_max+1)
coeff_FF[2] <- 10


AAA = 1
OMEGA = pi

alpha <- 1.8
m = 4
beta <- alpha/2


time_step <- 0.05
time_seq <- seq(0, T_final, length.out = ((T_final - 0) / time_step + 1))
h <- 0.05
globe_coarse <- globe_from_h(h)
globe_fine <- 2*globe_coarse

mesh_coarse <- fm_rcdt_2d(globe = globe_coarse)
FEM_coarse <- fm_fem(mesh_coarse, order = 1)
C <- FEM_coarse$c1
G <- FEM_coarse$g1
L <- kappa^2 * C + G

eig <- compute_true_eigen_sphere(
  mesh = mesh_coarse, 
  kappa = kappa, 
  alpha = alpha,
  L_max = L_max, 
  rot.inv = rot.inv)
EIGENVAL_ALPHA <- eig$EIGENVAL_ALPHA
EIGENFUN <- eig$EIGENFUN
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]))

mesh_fine <- fm_rcdt_2d(globe = globe_fine)
FEM_fine <- fm_fem(mesh_fine, order = 1)
Psi <- fm_basis(mesh_coarse, mesh_fine$loc)

eig_fine <- compute_true_eigen_sphere(
  mesh = mesh_fine, 
  kappa = kappa, 
  alpha = alpha, 
  L_max = L_max, 
  rot.inv = rot.inv)
EIGENFUN_fine <- eig_fine$EIGENFUN

FF_approx <- t(Psi) %*% 
  FEM_fine$c1 %*%
  EIGENFUN_fine %*%
      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))

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_approx)


diff <- U_true - U_approx
error <- sqrt(as.double(time_step * sum(diff * (C %*% diff))))
error

## [1] 0.1540148

Below we visualize the first few eigenfunctions on the sphere.

plot_3d_slider_sphere(mesh_coarse, eig$EIGENVAL, EIGENFUN, fixed_colorscale = TRUE)

Figure 1: Eigenfunctions on the sphere.

Below we show the absolute value of the difference \(|u(s, 0.5) - U^\tau_{h,m}(s,0.5)|\) and the function value \(U^\tau_{h,m}(s, 0.5)\) at \(t=0.5\).

idx <- which.min(abs(time_seq - 0.5))
error_sphere <- plot_3d_sphere_surface_onecol(mesh_coarse,  abs(diff[,idx]))
save(error_sphere, file = here::here("data_files/exp_1_error_sphere.RData"))
error_sphere

Figure 2: Absolute difference between the exact and approximate solution at \(t_0=0.5\).

f_sphere <- plot_3d_sphere_surface_onecol(mesh_coarse,  U_approx[,idx])
f_sphere

Figure 3: Approximate solution at \(t_0=0.5\).

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