Smoothing splines with a boundary constraint Is there a way to set a boundary constraint when one uses smoothing splines? 
For example, I usually fit the data with the following code using the R package mgcv :
library(mgcv)
x <- runif(100, 1, 9)
y <- sin((x^2)/10)+rnorm(100)
plot(gam(y~s(x)))

Suppose I would like to set a boundary constraint at the two points of x at 0 and 10 so that y = 0 when x = 0, 10. Notice that the two end points are outside the range [1, 9] of x. Is this achievable with mgcv?
 A: As pointed in the comments, the pc argument of the s() function included in the mgcv package does not allow for multiple constraint points. This is unfortunate but I think should not be too complicated to achieve the objective outside the realm of the specific package.
Intro
I think we can obtain the result we wish using two strategy:


*

*approximate the conditions (we will not get exactly $y = 0$ at $x = [0, 10]$)

*exact constraints


The first strategy has the advantage to allow for easy inference and can also be easily translate in Bayesian settings if one wish so (and might also be possible to achieve within mgcv but I am not a super expert of the package). However I will not go so much into the details but I will  point to some reference.
I will discuss both solutions using P-splines smoothing as introduced by Eilers and Marx, 1991 (option bs = ps in s() function). P-splines combine B-spline bases and finite difference penalties (you can read more about this here and here...please give a look at the extrapolation properties of the P-splines because it is relevant in your case).
In what follows I will indicate with $B$ a matrix of B-spline bases, with $P$ a finite difference penalty matrix and with $\lambda$ the smoothing parameter (I will keep it fixed for convenience in the codes). 
Strategy 1 - extra penalty
This 'trick' consists in adding an extra penalty term to the penalized problem. The penalized least squares problem becomes then
$$
\min_{c} S_{p} = \|y - B c\|^{2} + \lambda c^{\top} P c + \kappa (\Gamma c - v(x_{0}))^{\top}  (\Gamma c - v(x_{0}))
$$
where $\Gamma$ is the B-spline functions evaluated at the boundary points, $\kappa$ is a large constant (say $10^8$) and $v(x_{0})$ are the boundary abscissa (n your case a vector of zero of dim 2).
Strategy 2 - Lagrange multipliers
The previous strategy gives only a sort of 'soft' approximations. We can obtain an exact match using Lagrange multipliers (a reference in this context is here). In this case the penalized least squares problem is slightly different:
$$
\min_{c} S_{l} = \|y - B c\|^{2} + \lambda c^{\top} P c + \gamma^{\top} (\Gamma c - v(x_{0}))
$$ 
where $\gamma$ is a vector of Lagrange multipliers to be estimated.
A small R code
I will now use both strategies to smooth your data. I hope the code is clear enough (anyway I left some comments to guide you). The code supposes that you have a function to compute the B-splines $B$ (see for example Eilers and Marx, 2010). 
  rm(list =ls()); graphics.off()

  # Simulate some data
  set.seed(2020)
  xmin = 1
  xmax = 9
  m    = 200
  x    = seq(xmin, xmax, length = m)
  ys   = sin((x^2)/10)
  y    = ys+rnorm(m) * 0.2

  # Boundary conditions
  bx = c(0, 10)
  by = c(0, 0)

  # Compute bases for function, first and second derivative 
  bdeg = 3
  nseg = 50
  B0 = bbase(x, bx[1], bx[2], nseg, bdeg)
  nb = ncol(B0)
  Gi = bbase(bx, bx[1], bx[2], nseg, bdeg)


  # Set syste penalty and  with extra penalty
  D   = diff(diag(nb), diff = 2)
  P   = t(D) %*% D
  Bb  = t(B0) %*% B0
  Ci  = t(Gi) %*% Gi
  lam = 1e1
  kap = 1e8

  # Solve system strategy 1
  cof_p = solve(Bb + lam * P + kap * Ci) %*% (t(B0) %*% y + kap * t(Gi) %*% by)

  # Solve system strategy 2
  LS = rbind((Bb + lam * P), Gi)
  RS = rbind(t(Gi), 0 * diag(0, nrow(Gi)))
  cof_l = solve(cbind(LS, RS)) %*% c(t(B0) %*% y, by)

  # Plot results
  plot(x, y, xlim = bx, pch = 16)
  lines(x, ys, col = 8, lwd = 2)
  points(bx, by, pch = 15)

  # Strategy 1
  lines(x, B0 %*% cof_p, lwd = 2, col = 2)
  points(bx[1], (Gi %*% cof_p)[1], col = 2, pch = 16)
  points(bx[2], (Gi %*% cof_p)[2], col = 2, pch = 16)

  # Strategy 2
  lines(x, B0 %*% cof_l[1:nb], lwd = 2, col = 3, lty = 2)
  points(bx[1], (Gi %*% cof_l[1:nb])[1], col = 3, pch = 16, cex = 0.75)
  points(bx[2], (Gi %*% cof_l[1:nb])[2], col = 3, pch = 16, cex = 0.75)

  legend('bottomleft', c('data', 'signal', 'strategy1', 'strategy2'), col = c(1, 8, 2, 3), pch = 16)

The final results should look like this:

I hope this helps somehow.
