[R] adaptive optimization of mesh size

From: baptiste Auguié <ba208_at_exeter.ac.uk>
Date: Sun, 04 May 2008 16:39:43 +0100


DeaR list,

I'm running an external program that computes some electromagnetic response of a scattering body. The numerical scheme is based on a discretization with a characteristic mesh size "y". The smaller y is, the better the result (but obviously the computation will take longer).

A convergence study showed the error between the computed values and the exact solution of the problem to be a quadratic in y, with standard error increasing as y^3. I wrote the interface to the program in R, as it is much more user friendly and allows for post- processing analysis. Currently, it only runs with user-defined discretization parameter. I would like to implement an adaptive scheme [1] and provide the following improvements,

  1. obtain an estimate of the error by fitting the result against a series of mesh sizes with the quadratic model, and extrapolate at y = 0. (quite straight forward)
  2. adapt "dynamically" the set of mesh sizes to fulfill a final accuracy condition, between a starting value (a rule-of thumb estimate is given by the problem values). The lower limit of y should also be constrained by the resources (again, an empirical rule dictates the computation time and memory usage).

I'm looking for advice on this second point (both on the technical aspect, and whether this is sound statistically):

Below are a few lines of code to illustrate the problem,

> program.result <- function(x, p){ # made up function that mimicks
> the results of the real program
>
> y <- p[3]*x^2 + p[2]*x + p[1]
> y * (1 + rnorm(1, mean=0, sd = 0.1 * y^3))
> }
>
>
> p0 <- c(0.1, 0.1, 2) # set of parameters
>
> ## user defined limits of the y parameter (log scale)
> limits <- c(0.1, 0.8)
> limits.log <- (10^limits)
> y.log <- seq(limits.log[1], limits.log[2], l=10)
>
> y <- log10(y.log)
>
> result <- sapply(y, function(x) program.result(x, p0)) # results of
> the program
>
> #### fitting and extrapolation procedure ####
> library(gplots) # plot with CI
> plotCI(y, result, y^3, xlim=c(0, 1), ylim=c(0, 2)) # the data
> with y^3 errors
>
> my.data <- data.frame(y = y, value = result)
>
> fm <- lm(value ~ poly(y, degree=2, raw=TRUE), data = my.data ,
> weights = 1/y^3)
> lines(y, predict(fm, data.frame(y=y)), col = 2)
>
> extrap <- summary(fm)$coefficients[1,] # intercept and error on it
> plotCI(0,extrap[1], 2 * extrap[2], col = 2, add=T)
>
>
> ### my naive take on adaptive runs... ##
>
> objective <- 1e-3 # stop when the standard error of the
> extrapolated value is smaller than this
>
> err <- extrap[2]
> my.color <- 3
> while (err > objective){
>
> new.value <- min(y)/2 # i don't know how to choose this optimally
> y <- c(new.value, y)
> new.result <- program.result(new.value, p0)
> result <- c(new.result, result)
> points(new.value, new.result, col= my.color)
> my.data <- data.frame(y = y, value = result)
> fm <- lm(value ~ poly(y, degree=2, raw=TRUE), data = my.data ,
> weights = 1/y^3)
> lines(y, predict(fm, data.frame(y=y)), col = my.color)
>
> extrap <- summary(fm)$coefficients[1,] # intercept and error on it
> err <- extrap[2]
> print(err)
> plotCI(0,extrap[1], 2 * err, col = 2, add=T)
> my.color <- my.color + 1
>
> }
> err

Many thanks in advance for your comments,

baptiste

[1]: Yurkin et al., Convergence of the discrete dipole approximation. II. An extrapolation technique to increase the accuracy. J. Opt. Soc. Am. A / Vol. 23, No. 10 / October 2006


Baptiste Auguié

Physics Department
University of Exeter
Stocker Road,
Exeter, Devon,
EX4 4QL, UK

Phone: +44 1392 264187

http://newton.ex.ac.uk/research/emag
http://projects.ex.ac.uk/atto



R-help_at_r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Sun 04 May 2008 - 16:38:59 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Tue 06 May 2008 - 10:30:33 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.

list of date sections of archive