I am running physics simulations. Sometimes a simulation has insufficient resolution and produces discontinuities in the output variables. These are very easy to spot by eye in a scatterplot. I would like to automate their detection, and rerun any simulation with a detected discontinuity using a finer resolution.
For the sake of simplicity, let's say there is one variable $x$ that varies with time $t$. The variable $x$ can go up and down, but it should vary continuously and smoothly with $t$. However, I only get to sample $x$ at irregular intervals $t_1, \ldots, t_n$. Here are examples of ok and not ok series:
t <- sort(runif(20, 0, pi))
x <- t**2 # ok
x <- c(sin(t[1:10]), cos(t[11:20])) # not ok
The first thing I tried doing is to examine the differences in $x$ and see if any of them have changes that are much bigger than the changes between other points. In R, this corresponds to:
length(boxplot.stats(diff(x), coef=10)$out) > 0
This works ok, but unfortunately, some time steps, especially at the beginning and the end of the simulation, can be large; meaning that even when it works, there are large jumps in $x$ that don't arise from discontinuities.
The next thing I tried to do was fit splines to $x$ and look at derivatives, for example:
new_ts <- seq(min(t), max(t), length=1000)
splinefun(t, x)(new_ts, 3) # obtains the jerk of x with respect to t
I then look at outliers in the same way as before. This works sometimes, but it seems very difficult to calibrate how far the jerk needs to be from the rest of the sample in order to say that there has been a jump.
The next thing I'm trying is to fit a function to $x$ then then see if there are outliers in the residuals. However, I don't know what form $x$ will take analytically -- that's why I'm running the simulation! I could use a black box like a neural network, but it might just learn the discontinuity instead of detecting it.
So, does anyone have any advice?
