Skip to content

Instantly share code, notes, and snippets.

@Mathias-Fuchs
Last active August 19, 2018 08:50
Show Gist options
  • Save Mathias-Fuchs/b127324f67135ee23b5628ea15924c0f to your computer and use it in GitHub Desktop.
Save Mathias-Fuchs/b127324f67135ee23b5628ea15924c0f to your computer and use it in GitHub Desktop.
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
*~
.Rhistory
Rplot*.*
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)))
}
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"))
}
@Mathias-Fuchs
Copy link
Author

screenshot_2018-08-15_17-28-21

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment