a snippet of R code showing the convergence of a Fourier series against a sum of Dirac measures. The spatstat package is only used for visualisation

.gitignore

*~
.Rhistory
Rplot*.*

Goodness2.R

require(spatstat)
window <- owin(c(0, 1), c(0, 1))
maxM <- 30
eig <- function(k, m, x, y) sin(k * pi * x) * sin(m * pi * y)
MakePredictedFunction <- function(chi) {
    f <- function(x, y) {
        a <- 0
        for (k in 1:maxM) {
            for (m in 1:maxM) {
                a <- a + chi[k, m] * eig(x = x, y = y, k = k, m = m)
            }
        }
        a
    }
    f
}

circleParamX <- function(t) cos(2 * pi * t)/3 + 1/2
circleParamY <- function(t) sin(2 * pi * t)/3 + 1/2
smallCircleLeftParamX <- function(t) cos(2 * pi * t)/6 + 1/4
smallCircleLeftParamY <- function(t) sin(2 * pi * t)/6 + 1/2
smallCircleRightParamX <- function(t) cos(2 * pi * t)/6 + 3/4
smallCircleRightParamY <- function(t) sin(2 * pi * t)/6 + 1/2
twoCirclesParamX <- function(t) {
    if (length(t) == 1) {
        if (t<1/2) {return(smallCircleLeftParamX(2*t))} else {return(smallCircleRightParamX(2*t-1))} }
    r <- rep(NA, length(t))
    ww <- t < 1/2
    r[ww] <- smallCircleLeftParamX(2*t[ww])
    r[!ww] <- smallCircleRightParamX(2*t[!ww])
    r
}
twoCirclesParamY <- function(t) {
    if (length(t) == 1) {
        if (t<1/2) {return(smallCircleLeftParamY(2*t))} else {return(smallCircleRightParamY(2*t-1))} }
    r <- rep(NA, length(t))
    ww <- t < 1/2
    r[ww] <- smallCircleLeftParamY(2*t[ww])
    r[!ww] <- smallCircleRightParamY(2*t[!ww])
    r
}

                                        # start of homotopy
integrand <- function(t, k, m) eig(k=k, m=m, x=twoCirclesParamX(t), twoCirclesParamY(t))
chi0 <- matrix(0.0, ncol = maxM, nrow = maxM)
for (k in 1:maxM) {
    for (m in 1:maxM) {
        chi0[k, m] <- integrate(integrand, 0, 1, k=k, m=m)$value
    }
}
                                        # try out this one: plot(as.im(chi0))
predictedFunction <- MakePredictedFunction(chi0)
pf <- as.im(predictedFunction, W = window)
pp <- ppp(x=twoCirclesParamX(seq(0,1,length.out=1000)), y=twoCirclesParamY(seq(0,1,length.out=1000)), W = window)
plottees <- list(pf,pp)
plot(as.layered(as.solist(plottees)), main = paste("partial sum of Fourier series comprising", maxM^2, "terms"))



                                        # end of homotopy
integrand <- function(t, k, m) eig(k=k, m=m, x=circleParamX(t), circleParamY(t))
chi1 <- matrix(0.0, ncol = maxM, nrow = maxM)
for (k in 1:maxM) {
    for (m in 1:maxM) {
        chi1[k, m] <- integrate(integrand, 0, 1, k=k, m=m)$value
    }
}
                                        # try out this one: plot(as.im(chi1))
predictedFunction <- MakePredictedFunction(chi1)
pf <- as.im(predictedFunction, W = window)
pp <- ppp(x=circleParamX(seq(0,1,length.out=1000)), y=circleParamY(seq(0,1,length.out=1000)), W = window)
plottees <- list(pf,pp)
plot(as.layered(as.solist(plottees)), main = paste("partial sum of Fourier series comprising", maxM^2, "terms"))




for (j in seq(0,1,length.out=40)) {
    chi <- (1-j) * chi0 + j*chi1
    predictedFunction <- MakePredictedFunction(chi)
    pf <- as.im(predictedFunction, W = window)
    plottees <- list(pf)
    plot(as.layered(as.solist(plottees)))
}

GoodnessExperiment.R

require(spatstat)
eig <- function(k, m, x, y) sin(k * pi * x) * sin(m * pi * y)
xs <- runif(14)
ys <- runif(14)
D <- data.frame(x = xs, y = ys)
p <- ppp(xs, ys)

par(mfrow = c(3, 3))

for (i in 2:10) {
    maxM <- i
    chi <- matrix(0.0, ncol = maxM, nrow = maxM)
    for (k in 1:maxM) {
        for (m in 1:maxM) {
            chi[k, m] <- mean(apply(D, 1, function(z) eig(x = z[1], y = z[2], k = k, m = m)))
        }
    }
    predictedFunction <- function(x, y) {
        a <- 0
        for (k in 1:maxM) {
            for (m in 1:maxM) {
                a <- a + chi[k, m] * eig(x = x, y = y, k = k, m = m)
            }
        }
        a
    }
    pf <- as.im(predictedFunction, W = Window(p))
    plottees <- list(pf, p)
    plot(as.layered(as.solist(plottees)), main = paste("partial sum of Fourier series comprising", maxM^2, "terms"))
}