2 added 1110 characters in body

Let $c(r)$ be a nice, continuous function with compact support. For example, $c(r) = \tfrac 1 5 (1-r)^{11} \big( 5 + 55r + 239 r^2 + 429 r^3 \big)$ for $r \in [0,1]$, and $c(r) = 0$ otherwise.

On a rectangular domain $D$ in $\mathbb R^2$, we define $c(x,y) = c(|x-y|)$, which we interpret as the kernel of an integral operator. Consider the integral equation $$\iint_D c(x,y) f(y) \ \mathrm{d}y = \lambda f(x).$$ By Mercer's theorem, we can find a sequence of non-negative eigenvalues $\lambda_k$ and eigenfunctions $f_k(x)$ such that:

• The sequence $f_k$ is an orthonormal basis for $L^2(D)$,

• The eigenfunctions $f_k$ corresponding to non-negative eigenvalues are continuous on $D$, and

• We have the representation $c(x,y) = \sum \lambda_k f_k(x) f_k(y)$.

Question: How can we numerically approximate the first $K$ eigenfunctions $f_1, \cdots, f_K$ and keep the total error within a desired error? (Polynomial or Fourier approximations would suffice)

Edit:A number of commentators have described the natural technique: discretize space, turn the covariance into a matrix, easily solve for the eigenvectors, then reconstruct the functions $f_k(t)$ by interpolation.

Our motivation is to generate a Gaussian random field $\xi$ using the Karhunen-Loève transformation: let $Z_1, \cdots, Z_K$ be independent Gaussian RVs, and define the random field $$\xi(t) = \sum_{k=1}^K Z_k f_k(t).$$ Since the covariance is so smooth, the eigenfunctions will be smooth, and the Gaussian field will be smooth. I'm sure that a competent numerical analyst could use the discretization technique to reconstruct the eigenfunctions up to arbitrary derivatives. Unfortunately, I am an incompetent numerical analyst, and I would surely introduce systematic errors along the way. Yet the problem is so common that I am still hoping somebody could point us to a simple reference which does this well.

To make the question more mathematical precise: the numerical approximations $f_k$ should be within an arbitrarily small $C^3$-error of the actual eigenfunctions.

1

# Numerically finding a Mercer expansion for a given covariance kernel

Let $c(r)$ be a nice, continuous function with compact support. For example, $c(r) = \tfrac 1 5 (1-r)^{11} \big( 5 + 55r + 239 r^2 + 429 r^3 \big)$ for $r \in [0,1]$, and $c(r) = 0$ otherwise.

On a rectangular domain $D$ in $\mathbb R^2$, we define $c(x,y) = c(|x-y|)$, which we interpret as the kernel of an integral operator. Consider the integral equation $$\iint_D c(x,y) f(y) \ \mathrm{d}y = \lambda f(x).$$ By Mercer's theorem, we can find a sequence of non-negative eigenvalues $\lambda_k$ and eigenfunctions $f_k(x)$ such that:

• The sequence $f_k$ is an orthonormal basis for $L^2(D)$,

• The eigenfunctions $f_k$ corresponding to non-negative eigenvalues are continuous on $D$, and

• We have the representation $c(x,y) = \sum \lambda_k f_k(x) f_k(y)$.

Question: How can we numerically approximate the first $K$ eigenfunctions $f_1, \cdots, f_K$ and keep the total error within a desired error? (Polynomial or Fourier approximations would suffice)