Convergence in τ for the optimal control variable

Go back to the Contents page.

---

Press Show to reveal the code chunks.

---

# Create a clipboard button on the rendered HTML page
source(here::here("clipboard.R")); clipboard

# 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)
}

library(MetricGraph)
library(ggplot2)
library(reshape2)
library(plotly)
library(patchwork)
library(slackr)
source("keys.R")
slackr_setup(token = token) # token comes from keys.R

## [1] "Successfully connected to Slack"

capture.output(
  knitr::purl(here::here("functionality.Rmd"), output = here::here("functionality.R")),
  file = here::here("purl_log.txt")
)
source(here::here("functionality.R"))

# Parameters
T_final <- 2
kappa <- readRDS("old/kappa.RDS") #4
divider <- readRDS("old/divider.RDS") 
mu <- 0.1
a <- - 0.5
b <- 0.5
N_finite = 4 # choose even
adjusted_N_finite <- N_finite + N_finite/2 + 1
# Coefficients for f and g
coeff_elliptic_g <- 20*(1:adjusted_N_finite)^-1
coeff_elliptic_g[-2] <- 0
coeff_elliptic_f <- rep(0, adjusted_N_finite)
coeff_elliptic_f[2] <- 10
# Time step and mesh size
POWERS <- seq(6.125, 5.625, by = -0.1)
tau_vector <- 0.1 * 2^-POWERS

# Overkill parameters
overkill_time_step <- 0.1 * 2^-14
overkill_h <- (0.1 * 2^-14)^(1/2)

# 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

overkill_A <- matrix(a, nrow = length(overkill_weights), ncol = length(overkill_time_seq))
overkill_B <- matrix(b, nrow = length(overkill_weights), ncol = length(overkill_time_seq))

overkill_psi <- cos(overkill_time_seq)
overkill_psi_prime <- - sin(overkill_time_seq)
overkill_phi <- sin(T_final - overkill_time_seq)
overkill_phi_prime <- - cos(T_final - overkill_time_seq)
# psi_rate <- -2
# phi_rate <- -2
# phi_coefficient <- 1
# overkill_psi <- exp(psi_rate*overkill_time_seq)
# overkill_psi_prime <- psi_rate * exp(psi_rate*overkill_time_seq)
# overkill_phi <- phi_coefficient * exp(psi_rate*overkill_time_seq) - phi_coefficient*exp(psi_rate*T_final)
# overkill_phi_prime <- phi_coefficient * psi_rate * exp(psi_rate*overkill_time_seq)

m_values <- c()
alpha_vector <- seq(1, 1.8, by = 0.2)

# Create a matrices to store the errors
errors_u_bar <- matrix(NA, nrow = length(tau_vector), ncol = length(alpha_vector))
errors_p_bar <- matrix(NA, nrow = length(tau_vector), ncol = length(alpha_vector))
errors_z_bar <- matrix(NA, nrow = length(tau_vector), ncol = length(alpha_vector))
for (j in 1:length(alpha_vector)) {
  alpha <- alpha_vector[j] 
  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_MINUS_ALPHA <- overkill_eigen_params$EIGENVAL_MINUS_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_elliptic_f <- as.vector(overkill_EIGENFUN %*% coeff_elliptic_f)
  overkill_elliptic_g <- as.vector(overkill_EIGENFUN %*% coeff_elliptic_g)
  # Construct the corresponding elliptic solution u and v on the integration mesh
  overkill_elliptic_u <- as.vector(overkill_EIGENFUN %*% (coeff_elliptic_f * EIGENVAL_MINUS_ALPHA))
  overkill_elliptic_v <- as.vector(overkill_EIGENFUN %*% (coeff_elliptic_g * EIGENVAL_MINUS_ALPHA))
  overkill_u_bar <- outer(overkill_elliptic_u, overkill_psi)
  overkill_p_bar <- - mu * outer(overkill_elliptic_v, overkill_phi)
  overkill_z_bar <- pmax(overkill_A, pmin(overkill_B, - overkill_p_bar / mu))
  for (i in 1:length(tau_vector)) {
    time_step <- tau_vector[i]
    h <- (time_step^(1/alpha))/divider
    m <- min(4, ceiling((log(h))^2 / (pi^2 * (1 - alpha / 2)))) # min(4, ceiling(25 * (log(h))^2 / (16 * pi^2 * (1 - alpha / 2))))
    m_values <- c(m_values, m)
    time_seq <- seq(0, T_final, length.out = ((T_final - 0) / time_step + 1))
    graph <- gets.graph.tadpole(h = h)
    graph$compute_fem()
    G <- graph$mesh$G
    C <- graph$mesh$C
    L <- kappa^2*C + G
    # Construct the fractional operator, which is shared for the forward and adjoint problems
    my_op_frac <- my.fractional.operators.frac(L, 
                                               beta, 
                                               C, 
                                               scale.factor = kappa^2, 
                                               m = m, 
                                               time_step)
    eigen_params <- gets.eigen.params(N_finite = N_finite, 
                                      kappa = kappa, 
                                      alpha = alpha, 
                                      graph = graph)
    EIGENFUN <- eigen_params$EIGENFUN
    # Construct the right hand side functions f and g for the elliptic problem
    elliptic_f <- as.vector(EIGENFUN %*% coeff_elliptic_f)
    elliptic_g <- as.vector(EIGENFUN %*% coeff_elliptic_g)
    # Construct the corresponding elliptic solution u and v
    elliptic_u <- as.vector(EIGENFUN %*% (coeff_elliptic_f * EIGENVAL_MINUS_ALPHA))
    elliptic_v <- as.vector(EIGENFUN %*% (coeff_elliptic_g * EIGENVAL_MINUS_ALPHA))
    
    psi <- cos(time_seq)
    psi_prime <- - sin(time_seq)
    phi <- sin(T_final - time_seq)
    phi_prime <- - cos(T_final - time_seq)
    # psi <- exp(psi_rate*time_seq)
    # psi_prime <- psi_rate * exp(psi_rate*time_seq)
    # phi <- phi_coefficient * exp(phi_rate*time_seq) - phi_coefficient * exp(phi_rate*T_final)
    # phi_prime <- phi_coefficient * phi_rate * exp(phi_rate*time_seq)
    # Construct the projection matrix
    Psi <- graph$fem_basis(overkill_graph$get_mesh_locations())
    R <- t(Psi) %*% overkill_graph$mesh$C
    
    f_plus_z_bar <- outer(elliptic_u, psi_prime) + outer(elliptic_f, psi)
    f_plus_z_bar_innerproduct <- R %*% Psi %*% f_plus_z_bar
    u_0 <- elliptic_u
    # Solve the forward problem
    u_bar <- solve_fractional_evolution(my_op_frac,
                                        time_step,
                                        time_seq,
                                        val_at_0 = u_0,
                                        RHST = f_plus_z_bar_innerproduct)
    # u_bar <- outer(elliptic_u, psi)
    
    u_d <- outer(elliptic_u, psi) -
      mu * outer(elliptic_v, phi_prime) +
      mu * outer(elliptic_g, phi)
    v_d <- reversecolumns(R %*% Psi %*% u_d)
    v_bar <- reversecolumns(R %*% Psi %*% u_bar)
    # Solve the adjoint problem
    q_bar <- solve_fractional_evolution(my_op_frac, 
                                        time_step, 
                                        time_seq, 
                                        val_at_0 = u_0*0, 
                                        RHST = v_bar - v_d)
    p_bar <- reversecolumns(q_bar)
    # Compute the control variable
    A <- matrix(a, nrow = nrow(C), ncol = length(time_seq))
    B <- matrix(b, nrow = nrow(C), ncol = length(time_seq))
    z_bar <- pmax(A, pmin(B, - p_bar / mu))
    
    projected_u_bar_piecewise <- construct_piecewise_projection(Psi %*% u_bar, time_seq, overkill_time_seq)
    projected_p_bar_piecewise <- construct_piecewise_projection(Psi %*% p_bar, time_seq, overkill_time_seq)
    projected_z_bar_piecewise <- construct_piecewise_projection(Psi %*% z_bar, time_seq, overkill_time_seq)
    
    errors_u_bar[i,j] <- sqrt(as.double(t(overkill_weights) %*% (overkill_u_bar - projected_u_bar_piecewise)^2 %*% rep(overkill_time_step, length(overkill_time_seq))))
    errors_p_bar[i,j] <- sqrt(as.double(t(overkill_weights) %*% (overkill_p_bar - projected_p_bar_piecewise)^2 %*% rep(overkill_time_step, length(overkill_time_seq))))
    errors_z_bar[i,j] <- sqrt(as.double(t(overkill_weights) %*% (overkill_z_bar - projected_z_bar_piecewise)^2 %*% rep(overkill_time_step, length(overkill_time_seq))))
    print(paste0("m =", m, ", alpha =", alpha, ", h =", h, ", time_step =", time_step))
    slackr_msg(text = paste0("m =", m, ", alpha =", alpha, ", h =", h, ", time_step =", time_step), channel = "#research")
  }
}

## [1] "m =4, alpha =1, h =0.00014328188175073, time_step =0.0014328188175073"
## [1] "m =4, alpha =1, h =0.000153565718522695, time_step =0.00153565718522695"
## [1] "m =4, alpha =1, h =0.000164587661867943, time_step =0.00164587661867943"
## [1] "m =4, alpha =1, h =0.000176400688250958, time_step =0.00176400688250958"
## [1] "m =4, alpha =1, h =0.000189061576440515, time_step =0.00189061576440515"
## [1] "m =4, alpha =1, h =0.00020263118041422, time_step =0.0020263118041422"
## [1] "m =4, alpha =1.2, h =0.000426736113080764, time_step =0.0014328188175073"
## [1] "m =4, alpha =1.2, h =0.000452111162839404, time_step =0.00153565718522695"
## [1] "m =4, alpha =1.2, h =0.000478995091576214, time_step =0.00164587661867943"
## [1] "m =4, alpha =1.2, h =0.00050747762190425, time_step =0.00176400688250958"
## [1] "m =4, alpha =1.2, h =0.000537653811620773, time_step =0.00189061576440515"
## [1] "m =4, alpha =1.2, h =0.000569624370953814, time_step =0.0020263118041422"
## [1] "m =4, alpha =1.4, h =0.00093048216814821, time_step =0.0014328188175073"
## [1] "m =4, alpha =1.4, h =0.000977710315330173, time_step =0.00153565718522695"
## [1] "m =4, alpha =1.4, h =0.00102733560451291, time_step =0.00164587661867943"
## [1] "m =4, alpha =1.4, h =0.00107947970656676, time_step =0.00176400688250958"
## [1] "m =4, alpha =1.4, h =0.00113427046796645, time_step =0.00189061576440515"
## [1] "m =4, alpha =1.4, h =0.00119184222424404, time_step =0.0020263118041422"
## [1] "m =4, alpha =1.6, h =0.00166962734076687, time_step =0.0014328188175073"
## [1] "m =4, alpha =1.6, h =0.00174354805838684, time_step =0.00153565718522695"
## [1] "m =4, alpha =1.6, h =0.0018207415257756, time_step =0.00164587661867943"
## [1] "m =4, alpha =1.6, h =0.00190135263994435, time_step =0.00176400688250958"
## [1] "m =4, alpha =1.6, h =0.00198553271304304, time_step =0.00189061576440515"
## [1] "m =4, alpha =1.6, h =0.00207343975638281, time_step =0.0020263118041422"
## [1] "m =4, alpha =1.8, h =0.00263091165418448, time_step =0.0014328188175073"
## [1] "m =4, alpha =1.8, h =0.00273419920948592, time_step =0.00153565718522695"
## [1] "m =4, alpha =1.8, h =0.00284154175426721, time_step =0.00164587661867943"
## [1] "m =4, alpha =1.8, h =0.00295309848427691, time_step =0.00176400688250958"
## [1] "m =4, alpha =1.8, h =0.00306903484516543, time_step =0.00189061576440515"
## [1] "m =4, alpha =1.8, h =0.00318952277785138, time_step =0.0020263118041422"

#save(errors_u_bar, errors_p_bar, errors_z_bar, file = here::here("data_files/control_error_tau.RData"))

1 Convergence results

# Load the errors data
#load(here::here("data_files/control_error_tau.RData"))

observed_rates_u_bar <- numeric(length(alpha_vector))
observed_rates_p_bar <- numeric(length(alpha_vector))
observed_rates_z_bar <- numeric(length(alpha_vector))

for (i in 1:length(alpha_vector)) {observed_rates_u_bar[i] <- coef(lm(log10(errors_u_bar[, i]) ~ log10(tau_vector)))[2]}
for (i in 1:length(alpha_vector)) {observed_rates_p_bar[i] <- coef(lm(log10(errors_p_bar[, i]) ~ log10(tau_vector)))[2]}
for (i in 1:length(alpha_vector)) {observed_rates_z_bar[i] <- coef(lm(log10(errors_z_bar[, i]) ~ log10(tau_vector)))[2]}
theoretical_rates <- rep(1, length(alpha_vector))

p_u_bar <- error.convergence.plotter(x_axis_vector = tau_vector, 
                               alpha_vector, 
                               errors_u_bar, 
                               theoretical_rates, 
                               observed_rates_u_bar,
                               line_equation_fun = loglog_line_equation,
                               fig_title = expression(italic(bar(u))),
                               x_axis_label = expression(italic(tau)))

p_p_bar <- error.convergence.plotter(x_axis_vector = tau_vector, 
                               alpha_vector, 
                               errors_p_bar, 
                               theoretical_rates, 
                               observed_rates_p_bar,
                               line_equation_fun = loglog_line_equation,
                               fig_title = expression(italic(bar(p))),
                               x_axis_label = expression(italic(tau)))

p_z_bar <- error.convergence.plotter(x_axis_vector = tau_vector, 
                               alpha_vector, 
                               errors_z_bar, 
                               theoretical_rates, 
                               observed_rates_z_bar,
                               line_equation_fun = loglog_line_equation,
                               fig_title = expression(italic(bar(z))),
                               x_axis_label = expression(italic(tau)))

p_all_tau <- (p_u_bar | p_p_bar | p_z_bar) + 
  plot_annotation(
    title = expression("         Convergence in " * italic(tau)),
    theme = theme(plot.title = element_text(size = 18, face = "bold", hjust = 0.5))
  )
p_all_tau

Figure 1: Comparison of theoretical and observed convergence behavior for the \(L_2((0,T);L_2(\Gamma))\)-error with respect to \(\tau\) on a \(\text{log}_{10}\)–\(\text{log}_{10}\) scale. Dashed lines indicate the theoretical rates, and solid lines represent the observed error curves. The legend below each plot shows the value of \(\alpha\) along with the corresponding theoretical (‘theo’), and observed (‘obs’) rates for each case.

# ggsave(here::here("data_files/control_conv_rates_tau_u_bar.png"), width = 4, height = 5, plot = p_u_bar, dpi = 300)
# ggsave(here::here("data_files/control_conv_rates_tau_p_bar.png"), width = 4, height = 5, plot = p_p_bar, dpi = 300)
# ggsave(here::here("data_files/control_conv_rates_tau_z_bar.png"), width = 4, height = 5, plot = p_z_bar, dpi = 300)
ggsave(here::here("data_files/control_conv_rates_tau_all.png"), width = 12, height = 6, plot = p_all_tau, dpi = 300)

2 References

cite_packages(output = "paragraph", out.dir = ".")

We used R version 4.5.0 (R Core Team 2025) and the following R packages: 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), knitr v. 1.50 (Xie 2014, 2015, 2025), 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), patchwork v. 1.3.1 (Pedersen 2025), plotly v. 4.10.4 (Sievert 2020), RColorBrewer v. 1.1.3 (Neuwirth 2022), renv v. 1.0.7 (Ushey and Wickham 2024), reshape2 v. 1.4.4 (Wickham 2007), rmarkdown v. 2.29 (Xie, Allaire, and Grolemund 2018; Xie, Dervieux, and Riederer 2020; Allaire et al. 2024), rSPDE v. 2.5.1.9000 (Bolin and Kirchner 2020; Bolin and Simas 2023; Bolin, Simas, and Xiong 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), xaringanExtra v. 0.8.0 (Aden-Buie and Warkentin 2024).

Aden-Buie, Garrick, and Matthew T. Warkentin. 2024. xaringanExtra: Extras and Extensions for “xaringan” Slides. https://doi.org/10.32614/CRAN.package.xaringanExtra.

Allaire, JJ, Yihui Xie, Christophe Dervieux, Jonathan McPherson, Javier Luraschi, Kevin Ushey, Aron Atkins, et al. 2024. rmarkdown: Dynamic Documents for r. https://github.com/rstudio/rmarkdown.

Bates, Douglas, Martin Maechler, and Mikael Jagan. 2025. Matrix: Sparse and Dense Matrix Classes and Methods. https://doi.org/10.32614/CRAN.package.Matrix.

Bolin, David, and Kristin Kirchner. 2020. “The Rational SPDE Approach for Gaussian Random Fields with General Smoothness.” Journal of Computational and Graphical Statistics 29 (2): 274–85. https://doi.org/10.1080/10618600.2019.1665537.

Bolin, David, Mihály Kovács, Vivek Kumar, and Alexandre B. Simas. 2024. “Regularity and Numerical Approximation of Fractional Elliptic Differential Equations on Compact Metric Graphs.” Mathematics of Computation 93 (349): 2439–72. https://doi.org/10.1090/mcom/3929.

Bolin, David, and Alexandre B. Simas. 2023. rSPDE: Rational Approximations of Fractional Stochastic Partial Differential Equations. https://CRAN.R-project.org/package=rSPDE.

Bolin, David, Alexandre B. Simas, and Jonas Wallin. 2023a. MetricGraph: Random Fields on Metric Graphs. https://CRAN.R-project.org/package=MetricGraph.

———. 2023b. “Statistical Inference for Gaussian Whittle-Matérn Fields on Metric Graphs.” arXiv Preprint arXiv:2304.10372. https://doi.org/10.48550/arXiv.2304.10372.

———. 2024. “Gaussian Whittle-Matérn Fields on Metric Graphs.” Bernoulli 30 (2): 1611–39. https://doi.org/10.3150/23-BEJ1647.

———. 2025. “Markov Properties of Gaussian Random Fields on Compact Metric Graphs.” Bernoulli. https://doi.org/10.48550/arXiv.2304.03190.

Bolin, David, Alexandre B. Simas, and Zhen Xiong. 2024. “Covariance-Based Rational Approximations of Fractional SPDEs for Computationally Efficient Bayesian Inference.” Journal of Computational and Graphical Statistics 33 (1): 64–74. https://doi.org/10.1080/10618600.2023.2231051.

Cheng, Joe, Carson Sievert, Barret Schloerke, Winston Chang, Yihui Xie, and Jeff Allen. 2024. htmltools: Tools for HTML. https://doi.org/10.32614/CRAN.package.htmltools.

Garnier, Simon, Ross, Noam, Rudis, Robert, Camargo, et al. 2023. viridis(Lite) - Colorblind-Friendly Color Maps for r. https://doi.org/10.5281/zenodo.4678327.

Kaye, Matt, Bob Rudis, Andrie de Vries, and Jonathan Sidi. 2025. slackr: Send Messages, Images, r Objects and Files to “Slack” Channels/Users. https://github.com/mrkaye97/slackr.

Müller, Kirill. 2020. here: A Simpler Way to Find Your Files. https://doi.org/10.32614/CRAN.package.here.

Neuwirth, Erich. 2022. RColorBrewer: ColorBrewer Palettes. https://doi.org/10.32614/CRAN.package.RColorBrewer.

Pedersen, Thomas Lin. 2025. patchwork: The Composer of Plots. https://doi.org/10.32614/CRAN.package.patchwork.

R Core Team. 2025. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. https://www.R-project.org/.

Sievert, Carson. 2020. Interactive Web-Based Data Visualization with r, Plotly, and Shiny. Chapman; Hall/CRC. https://plotly-r.com.

Ushey, Kevin, and Hadley Wickham. 2024. renv: Project Environments. https://doi.org/10.32614/CRAN.package.renv.

Van Boxtel, G.J.M., et al. 2021. gsignal: Signal Processing. https://github.com/gjmvanboxtel/gsignal.

Wickham, Hadley. 2007. “Reshaping Data with the reshape Package.” Journal of Statistical Software 21 (12): 1–20. http://www.jstatsoft.org/v21/i12/.

Wickham, Hadley, Mara Averick, Jennifer Bryan, Winston Chang, Lucy D’Agostino McGowan, Romain François, Garrett Grolemund, et al. 2019. “Welcome to the tidyverse.” Journal of Open Source Software 4 (43): 1686. https://doi.org/10.21105/joss.01686.

Wickham, Hadley, Thomas Lin Pedersen, and Dana Seidel. 2025. scales: Scale Functions for Visualization. https://doi.org/10.32614/CRAN.package.scales.

Xie, Yihui. 2014. “knitr: A Comprehensive Tool for Reproducible Research in R.” In Implementing Reproducible Computational Research, edited by Victoria Stodden, Friedrich Leisch, and Roger D. Peng. Chapman; Hall/CRC.

———. 2015. Dynamic Documents with R and Knitr. 2nd ed. Boca Raton, Florida: Chapman; Hall/CRC. https://yihui.org/knitr/.

———. 2025. knitr: A General-Purpose Package for Dynamic Report Generation in R. https://yihui.org/knitr/.

Xie, Yihui, J. J. Allaire, and Garrett Grolemund. 2018. R Markdown: The Definitive Guide. Boca Raton, Florida: Chapman; Hall/CRC. https://bookdown.org/yihui/rmarkdown.

Xie, Yihui, Christophe Dervieux, and Emily Riederer. 2020. R Markdown Cookbook. Boca Raton, Florida: Chapman; Hall/CRC. https://bookdown.org/yihui/rmarkdown-cookbook.