Another page [1] on this site describes how to design "Chebyshevoptimal" filters to calculate numerical estimates of the first derivative from samples of a function. This page extends that approach to estimate the second derivative. These filters are also optimal in the sense of minimizing the maximum deviation from the "ideal design" specification. Despite promising results, the designs are not the best choice for all applications.
Adapting the solution method
While the approach used for designing firstderivative approximators applies nicely, there are several details to take into account when applying the method to secondderivative approximators.

Filter response characteristic. Since the second derivative estimator is the equivalent of applying a firstderivative estimator twice, the frequencydomain response of the secondderivative filter will be the equivalent of two multiplications by
w
and two phase shifts byΠ/2
. Thus, in the critical lowfrequency region, the response curve is an upwardcurving quadratic curve. The combination of the two phase shifts is the equivalent of multiplying by a scalar factor1.0
. This is a somewhat harder problem numerically because the upward curvature must be countered by downward curvature to return the response curve somewhat close to zero at high frequencies, thus avoiding excessive noise response. The filters tend to need a few more terms than the firstderivative approximators to achieve similar accuracy. 
Basis functions. A second derivative filter will have even symmetry. Cosine waves have the desired properties. So, instead of using a sineseries approximation as in the firstderivative case, use a cosine series approximation.

No natural zero constraint. The second derivative of a constant sequence is zero, much like the first derivative. However, cosine "basis functions" do not have a natural property of passing through zero. Thus, a constraint that the combination of cosine functions must equal 0 at 0 frequency must be included, with rigorous tolerance on this requirement. Since the value of a cosine function is 1.0 when evaluated at 0, this constraint is equivalent to saying that the sum of the filter coefficients must be equal to 0.

Basis size. The number of terms in each of the basis vectors of matrix
B
isN/2  Nt
terms rather thanN/2  1  Nt
terms, to account for the added constraint at location 0.
After taking these adjustments into account, the solution method is fundamentally the same as for the firstderivative case. Full details of the method are provided on the page that covers the firstderivative case [1].
Here is the result of a typical design (see the 15term design below).
Design script
An mfile design script for Octave shows the numerical details of constructing the matrix terms, calculating the solution, and plotting the results for the secondderivative filter design. If you wish to experiment with this script, download a copy of the minmax2.m file. The Octave solution is based on the GLPK package, but you can use the equivalent solver from the Optimization Toolbox for Matlab [2].
Comparison to other designs
To see the practical benefits of the minmax designs, the 15term filter (as tabulated below) is compared to the 9term "Robust Low Noise" design from Pavel Holoborodko's site. While this comparison seems inherently unfair — the minmax design is allowed 60% more parameters to work with — this choice is pretty much forced.

The ninthorder "Robust Low Noise" designs is the highest order filter listed.

The "Robust Low Noise" designs deliver worse accuracy as the filter order increases, as you can clearly see from the frequency response plots for estimators of various orders — more parameters only makes the comparison more unfair.
First, here is the performance delivered by the RLN method in the presence of white random noise at 5% of the signal level, at a frequency of about 8% of the sampling rate (approximately 12 terms span a complete waveform cycle). A mathematically perfect result in green is plotted on top of the filter estimates in blue.
Next, here is the performance delivered by the minmax method, same level of random noise, same signal. (The apparent different number of terms is only an edge effect artifact, nothing of significance. When a longer filter is convolved with the same length data block, there are more invalid values at the ends where the filter "overlaps" beyond the ends of the data block.)
At first glance, it might seem like the minmax design is undisputed winner — but unfortunately, the choice is not that clear. While the minmax design tracks changing curvature and higher frequencies better, in situations where the curvature is small and changing slowly, a "Robust Low Noise" design produces superior results. The "RLN" approach concentrates maximum effort at matching low frequencies accurately. In contrast, the "minmax" design tolerates a small absolute error at every frequency, and this can result in a significant "relative error" when the size of the desired result is also small. For now, the minmax designs are not a complete and general solution, but very good a lot of the time.
Some tabulated results
The following table provides some design results that you can use directly without resorting to the design script. For each target number of filter terms, the design parameters, the resulting performance specs, and filter coefficients are listed. If you are not sure about where to start, try the 15term filter design with its 0.08·Fs
guaranteed accurate band. If you need more economy or better performance, switch to one of the other designs.
11term — 5% accurate — 8% useful 
Design: pass: 0.05 transition: 0.185 sensitivity: 100.0 
pass ripple: 0.001 stopband edge: 0.24 stopband ripple: 0.10 
filter: 0.024705 0.064369 0.070200 0.018820 0.070707 0.115954 0.070707 0.018820 0.070200 0.064369 0.024705

13term — 6% accurate — 9% useful 
Design: pass: 0.08 transition: 0.135 sensitivity: 100.0 
pass ripple: 0.001 stopband edge: 0.21 stop ripple: 0.104 
filter: 0.0303998 0.0062139 0.0870660 0.0549440 0.0219066 0.0862821 0.1068973 0.0862821 0.0219066 0.0549440 0.0870660 0.0062139 0.0303998

15term — 8% accurate — 11% useful 
Design: pass: 0.08 transition: 0.14 sensitivity: 150 
pass ripple: 0.00075 stopband edge: 0.22 stopband ripple: 0.1 
filter: 0.010053 0.051074 0.018345 0.070910 0.081264 0.017178 0.081808 0.129738 0.081808 0.017178 0.081264 0.070910 0.018345 0.051074 0.010053

17term — 10% accurate — 13% useful 
Design: pass: 0.10 transition: 0.16 sensitivity: 420 
pass ripple: 0.0002 stopband edge: 0.26 stopband ripple: 0.082 
filter: 0.0020745 0.0274629 0.0434860 0.0404178 0.0523018 0.1458368 0.0842230 0.1122293 0.2232337 0.1122293 0.0842230 0.1458368 0.0523018 0.0404178 0.0434860 0.0274629 0.0020745

Footnotes:
[1] The page on this site about firstderivative filter approximation has the details about to set up and solve the approximation problem.
[2] The "robust low noise estimators" are described at Pavel Holoborodko's web site, http://www.holoborodko.com/pavel/numericalmethods/numericalderivative/smoothlownoisedifferentiators/.
[3] You can find more information about the Optimization Toolkit at
http://www.mathworks.com/products/optimization/index.html .