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 now consider the following equation on the unit square \(D=[0,1]^2\), that is, \[\begin{equation} \label{eq:maineq_on_square_neumann} \tag{1} \left\{ \begin{aligned} \partial_t u(x,y,t) + (\kappa^2 - \Delta)^{\alpha/2} u(x,y,t) &= f(x,y,t), && \quad (x,y,t) \in D \times (0, T), \\ u(x,y,0) &= u_0(x,y), && \quad (x,y) \in D, \\ \partial_n u(x,y,t) &= 0, && \quad (x,y,t) \in \partial D \times (0,T), \end{aligned} \right. \end{equation}\] where \(\Delta\) denotes the Laplacian on \(D\) with homogeneous Neumann boundary conditions, and \(\partial_n\) denotes the outward normal derivative on \(\partial D\). The eigenvalues of \(-\Delta\) are known and given by \(\lambda_{m,n} = \pi^2(m^2+n^2)\) for \(m,n=0,1,2,\dots\), and the eigenfunctions are \(\phi_{m,n}(x,y)=\cos(m \pi x)\cos(n\pi y)\).

Below we use the eigenpairs to construct an exact solution to equation \(\eqref{eq:maineq_on_square_neumann}\) and compare it with the numerical approximation. We take \(u_0(x,y) = x_3 \phi_{1,0}(x,y)\) with \(x_3 = 1\) and \(f(x,y,t) = y_2 \sin(\pi t) \phi_{0,1}(x,y)\) with \(y_2 = 10\). Then the exact solution is given by \[\begin{equation} u(x,y,t) = x_3\text{e}^{-\lambda^{\alpha/2}_{1,0}t}\phi_{1,0}(x,y)+y_2 G_2(t)\text{e}^{-\lambda^{\alpha/2}_{0,1}t}\phi_{0,1}(x,y), \end{equation}\] where \(G_2(t)= \displaystyle\int_0^t \text{e}^{\lambda^{\alpha/2}_{0,1}r}\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

a <- b <- 1

m_max <- 1
n_max <- 1
N_eig <- (m_max + 1) * (n_max + 1)
coeff_U_0 <- rep(0, N_eig)
coeff_U_0[3] <- 1
coeff_FF <- rep(0, N_eig)
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
n_coarse <- ceiling(-1+sqrt(a*b)/h)
n_fine <- 2*n_coarse

mesh_and_FEM_coarse <- gets.mesh.and.FEM.on.rectangle(a, b, n_coarse)
mesh_coarse <- mesh_and_FEM_coarse$mesh
G <- mesh_and_FEM_coarse$G
C <- mesh_and_FEM_coarse$C
L <- kappa^2*C + G

eig <- compute_true_eigen_rectangle(
    a = a, 
    b = b, 
    loc = mesh_coarse$loc,
    kappa = kappa,
    alpha = alpha,
    m_max = m_max, 
    n_max = n_max)

m_vector = eig$m_vector
n_vector = eig$n_vector

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_and_FEM_fine <- gets.mesh.and.FEM.on.rectangle(a, b, n_fine)
mesh_fine <- mesh_and_FEM_fine$mesh
Psi <- fm_basis(mesh_coarse, mesh_fine$loc)

eig_fine <- compute_true_eigen_rectangle(
    a = a, 
    b = b, 
    loc = mesh_fine$loc,
    kappa = kappa,
    alpha = alpha,
    m_max = m_max, 
    n_max = n_max)
EIGENFUN_fine <- eig_fine$EIGENFUN

    # Compute matrix F with columns F^k
FF_approx <- t(Psi) %*% 
  mesh_and_FEM_fine$C %*%
  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.03092878

# Plot results
plot_3d_rectangle_scatter(mesh_coarse$loc, time_seq, abs(diff), fixed_colorscale = FALSE)

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

# Plot results
idx <- which.min(abs(time_seq - 0.5))
error_rectangle <- plot_3d_rectangle_onecol(mesh_coarse, as.vector(abs(diff[,idx])))
save(error_rectangle, file = here::here("data_files/exp_1_error_rectangle.RData"))
error_rectangle

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

f_rect <- plot_3d_rectangle_onecol(mesh_coarse, as.vector(U_approx[,idx]))
f_rect

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

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