`first_order` only differs from `second_order` at the endpoints

[pavelkomarov]

To calculate the first location, first_order does a diff diff, duplicates the first and last locations with an hstack, and then takes the mean of [:-1] and [1:] indexings of this array.

That gives the first location of the result $= \frac{\frac{x_1 - x_0}{dt} + \frac{x_1 - x_0}{dt}}{2} = \frac{x_1 - x_0}{dt}$.

But then the not-endpoints are doing second-order center-differencing, because

$$\frac{\frac{x_n - x_{n-1}}{dt} + \frac{x_{n+1} - x_n}{dt}}{2} = \frac{x_{n+1} - x_{n-1}}{2dt}$$.

This is the same as what second_order is doing for its central points. So first_order isn't truly first order, except at the endpoints.

[pavelkomarov]

Yeah, run this script:

from pynumdiff.finite_difference import first_order, second_order
import numpy as np

t = np.linspace(0, 3, 31)
x = t**2 - t + 1

print(first_order(x, 0.1)[1])
print(second_order(x, 0.1)[1])

You get back

[-0.9 -0.8 -0.6 -0.4 -0.2  0.   0.2  0.4  0.6  0.8  1.   1.2  1.4  1.6
  1.8  2.   2.2  2.4  2.6  2.8  3.   3.2  3.4  3.6  3.8  4.   4.2  4.4
  4.6  4.8  4.9]
[-1.  -0.8 -0.6 -0.4 -0.2  0.   0.2  0.4  0.6  0.8  1.   1.2  1.4  1.6
  1.8  2.   2.2  2.4  2.6  2.8  3.   3.2  3.4  3.6  3.8  4.   4.2  4.4
  4.6  4.8  5. ]

[pavelkomarov]

I've replaced the slightly-strangely-structured code in the finite difference module with my own. See #87. But now I'm seeing something funny in the tutorial notebook when I iterate.

Here are answers without noise, using standard 2nd order FD equations, first not iterated:

then iterated:

And here is the same for fourth order, not iterated:

and iterated:

I have confidence my FD implementations are correct. They match this tool; RMS of the above are super low; and the unit test plots look right.

Now let's add back in the noise. Here's second order:

iterated second order:

fourth order:

iterated fourth order:

Gross! Look at those edge effects. What's going on?

Well, the old first_order code that was actually second order, which was the thing getting iterated, contains a clue: The second-order endpoints are replaced with first-order ones. If I replace the endpoint formulas in my FD code with first-order diffs, here are the plots for second order:

iterated second order:

fourth order:

iterated fourth order:

Voila! But why!? I hear you cry. I believe the answer is that the endpoint formulas involve large constants, for example dxdt_hat[-1] = (25*x_hat[-1] - 48*x_hat[-2] + 36*x_hat[-3] - 16*x_hat[-4] + 3*x_hat[-5])/(12*dt) in 4th order. 48 > 12, so we're amplifying noise from the penultimate position! The center-difference formulas don't suffer this weakness, because they only contain constants lower than that in the denominator!

What do? I think I'll keep the fancy formulas around but only use them on the last iteration. This is an important caveat for these methods, that they're not truly $\nu^\text{th}$ order at the edges, although due to our trapezoidal integration, I'm not sure we can really say they're more than 1st order anyway.

[pavelkomarov]

Funny enough, iterated fourth-order is among the best methods in the repo for the cruise control data, competitive with TVR jerk and Kalman.

[pavelkomarov]

It's worth mentioning that I tried integrating with cumulative_simpson instead of cumulative_trapezoid, because second order is implicitly finding a quadratic interpolation, just as Simpson's rule does, and it made some sense to try to match the assumptions. But the results seem to be worse, even for second_order for most choices of num_iterations, and especially for fourth_order. I'm not entirely sure why.

Here are the iterated FD results from notebook 1 (num_iterations=50) with trapezoidal integration for 2nd order:

and 4th order:

And here they are with Simpson's rule for 2nd order: (pretty good, actually)

and 4th order:

Here are best results from notebook 2a using trapezoidal rule for 2nd order:
Optimal parameters: {'num_iterations': 32}
RMS error in velocity: 0.9787870657825298

and 4th order:
Optimal parameters: {'num_iterations': 72}
RMS error in velocity: 0.6732352858435942

And here they are with Simpson's rule for 2nd order:
Optimal parameters: {'num_iterations': 47}
RMS error in velocity: 0.9659739288632198

and 4th order:
Optimal parameters: {'num_iterations': 1000}
RMS error in velocity: 0.773696957008726

It's curious that 4th order can still do well, but it takes hitting the upper rail of num_iterations.

In the end, I believe what's going on is: "When the data is noisy, Simpson’s rule may overreact to curvature that isn’t real." This would explain why it takes more iterations to smooth out the answer with Simpson's rule.

[pavelkomarov]

Ahha. I'm reinvestigating this for the paper a bit, and I've found that essentially the FD->quadrature loop is implementing a low-pass IIR filter. The coefficients of this filter depend on the order of the schemes used. Simpson's rule is $\frac{1}{3}\Delta t (y_{0}+4y_{1}+y_{2})$, which has a 4/3 in there, which is >1 and risks amplifying/concentrating noise. The trapezoidal rule has only coefficients 1/2 and always smooths noise. Really it's the interplay of whichever orders of FD and quadrature you choose that make the filter input kernel, but maybe better to stay away from things that risk any amplification.