# Derivative of contour in ContourPlot

Sorry to bring this question up again, since there are many similar questions on the site. I use to think this is a easy job to do, because for the worst case I can follow the answers on this site. But after several tests, I found it is not as easy as I expected.

First is the contour plot:

Eq = (Sinh[
b] (80 - 80 Cos[b] + 8 b^4 Cos[b] - 8 b^4 Cos[b (1 - 2 xm)] +
10 b Sin[b] + 32 b^3 Sin[b] - 10 b Cosh[b] Sin[b] +
20 b Cosh[b/2]^2 (-1 + 2 Cosh[2 b xm]) Sin[b] -
20 b Sin[b (1 - 2 xm)] - 20 b Sin[2 b xm] -
20 b Cos[b] Sinh[b] + 20 b Cos[b (1 - 2 xm)] Sinh[b] -
80 Sin[b] Sinh[b] -
160 b Cosh[b/2] Sin[b/2] Sin[b xm] Sinh[b xm] +
20 b Sinh[2 b xm] - 20 b Cos[b] Sinh[2 b xm] +
8 b^4 Sin[b] Sinh[2 b xm] -
20 b Sin[b] Sinh[b] Sinh[2 b xm]) +
Cosh[b] (20 b - 50 b Cos[b] + 30 b Cos[b (1 - 2 xm)] +
10 b Cos[b] Cosh[b] - 10 b Cos[b (1 - 2 xm)] Cosh[b] -
120 Sin[b] + 8 b^4 Sin[b] + 40 Cosh[b] Sin[b] -
2 b Cosh[b xm]^2 (10 - 10 Cos[b] + 4 b^3 Sin[b]) +
20 b Sin[b] Sinh[b] - 20 b Cosh[2 b xm] Sin[b] Sinh[b] +
160 b Sin[b/2] Sin[b xm] Sinh[b/2] Sinh[b xm] -
20 b Sinh[b xm]^2 + 20 b Cos[b] Sinh[b xm]^2 -
8 b^4 Sin[b] Sinh[b xm]^2 - 30 b Sin[b] Sinh[2 b xm] +
10 b Cosh[b] Sin[b] Sinh[2 b xm] +
20 Cosh[b/2]^2 (b Cos[b] - b Cos[b (1 - 2 xm)] + 4 Sin[b] +
b Sin[b] Sinh[2 b xm])));

p = ContourPlot[Eq == 0, {xm, 0, 1/2}, {b, 2, 4}, PlotPoints -> 100]


Following suggestion from Mr.Wizard♦ in https://mathematica.stackexchange.com/a/31169/11867, I set PlotPoints to 100.

Then, extract the points from the plot:

lines = Cases[p, _Line, Infinity];
points = p[[1, 1]];
l1 = Map[Part[points, #] &, lines[[1, 1]]];


I think the curve is pretty smooth, so I use Interpolation to interpolate the points:

Plot[Interpolation[l1]'[x], {x, 0, 0.5}]


It is clear jitters in the data cause the derivative not smooth. Therefore, following stevenvh's advice in (http : // mathematica.stackexchange.com/a/10987/11867), I reduce the InterpolationOrder to 2 and 1.

Plot[Interpolation[l1, InterpolationOrder -> 2]'[x], {x, 0, 0.5}]
Plot[Interpolation[l1, InterpolationOrder -> 1]'[x], {x, 0, 0.5}]


It is better obviously, but not ideal. I think answer from chris make sense, BSpline should be a better solution, but not for derivative.

bs = BSplineFunction[l1];
ParametricPlot[bs[t], {t, 0, 1}]


ParametricPlot[bs'[t], {t, 0, 1}]


After that, I recall s.s.o mentioned Wavelet can be used to smooth the data in https://mathematica.stackexchange.com/a/65590/11867, but I suddenly realize that the data points are not evenly sampled. Therefore, I try to use other smooth method, like MovingAverage.

Interpolation[MovingAverage[l1, 51], InterpolationOrder -> 2]
Plot[%'[x], {x, 0, 1/2}]


This is pretty close to smooth curve, but not ideal. Then, I think MovingMedian may remove the jitter noise:

Interpolation[MovingMedian[l1, 15], InterpolationOrder -> 2]
Plot[%'[x], {x, 0, 1/2}]


I would like to narrow this problem to processing the data points in ContourPlot, because for this problem the derivative actually can be obtained through the analytical equation. I am curious about how to solve it by processing the data points in the figure.

## Summary

I would not like to call this a conclusion, because I think there may be a better solution out there. Untill now, we have mainly two possible solutions for this problem.

### Low pass filter based methods

@Jens and @Michael E2 both provide a solution based on low pass filter. It should be noted LowpassFilter and DifferentiatorFilter support only even sampled data for Mathematica lower than v10.1:

LowpassFilter currently only supports regularly sampled time series \ inputs. Use TimeSeriesResample or TemporalRegularity to make the \ input regularly sampled. >>

You can use code like following:

lp = LowpassFilter[TimeSeries@TimeSeriesResample[l1, 0.001], 0.5];


Note: the cut-off frequency in both filters will reduce the absolute value of derivative due to filter process. You should try several cut-off frequencies and make trade-off between smoothness and the reduction.

Moreover, I find MovingAverage may also provide smooth result if you increase PlotPoints to 300. The following command provides a decent result:

Interpolation[MovingAverage[l1, 30], InterpolationOrder -> 2];
Plot[%'[x], {x, 0, 1/2}]


Note: MovingAverage is also one kind of low-pass filter, and thus has the same problem as the previous two methods. However, it is fast and easy to understand.

### Spline based methods

There are two different implementations in the answer of @Alexey Popkov.

The first one is not strictly a spline based method, the idea is interpolate the data in sub-segements, and then fit the data with NonlinearModelFit.

The second one is based on QuantileRegression, which can be downloaded from here. The performance of the package looks quite decent.

BTW, the correct way to use BSplineFunction is shown in the comment by @Michael E2 as

dbs[t_?NumericQ] := First@Ratios[bs'[t]];
xbs[t_?NumericQ] :=
First@bs[t]; ParametricPlot[{xbs[t], dbs[t]}, {t, 0, 1},
AspectRatio -> 1/GoldenRatio]


The above codes produce following figure for PlotPoints->300

### Other thoughts

Piecewise polynomial should be a candidate for this problem, but as I commented in answer of @Alexey Popkov, implementation from Piecewise Polynomial Interpolation by @Michael E2 does not provide an ideal result for this problem. loess in R handles problems similar to this one pretty well. I find one implementation here, but I cannot make it work yet.

Currently, the best choices for this problem are QuantileRegression, LowpassFilter and MovingAverage.

• You'll need to smooth further. You may try ND[ ] with higher Terms option or some filter like GaussianFilter[ ] – Dr. belisarius Jun 5 '15 at 4:38
• I agree that this is probably a smoothing problem, but I also wonder: your l1 variable has exactly the same value as your points variable. Also, have you looked at the contents of lines? It is quite uninteresting as well. Are those expressions really doing what you want them to do? – MarcoB Jun 5 '15 at 4:46
• Related: "Implementation of smoothing splines function." – Alexey Popkov Jun 6 '15 at 9:50
• The way to plot the derivative $dy/dx$ of bs1 is something like dbs[t_?NumericQ] := First@Ratios[bs'[t]]; xbs[t_?NumericQ] := First@bs[t]; ParametricPlot[{xbs[t], dbs[t]}, {t, 0, 1}, AspectRatio -> 1/GoldenRatio], but it has noise similar to the Interpolation method. – Michael E2 Jun 6 '15 at 15:57
• The output of bs[t] is an xy pair if t is numeric. The output of bs'[t] is the velocity vector, again if t is numeric, whose components are $dx/dt$ and $dy/dt$. In other words, bs[t] is a parametric curve, but you can get at the components only numerically and not symbolically. Another way to look at bs is that it is a data structure that Mathematica uses to compute a value when the input is numeric. – Michael E2 Jun 6 '15 at 16:56

This seems to be a perfect candidate for DifferentiatorFilter:

t = TimeSeries[l1[[All, 2]], {l1[[All, 1]]}]

d = DifferentiatorFilter[t, .01];

ListLinePlot[d]


I just had to adjust the cutoff frequency to weed out the fast oscillations. However, this filter scales the output differently for different cutoff frequencies.

Unfortunately, I believe the scaling of the derivative is a bug. If you go to the Help documentation, under "Scope," and evaluate the example for filtering a time series, the plot that Mathematica produces is not what the documentation shows as expected output before the evaluation. The expected output is the derivative, and it is scaled correctly. But after evaluation, you get a derivative that looks like a flat line in comparison to the original function. This is because the scaling is completely off. You have to fix it by hand.

Anyway, if it weren't for that bug, this would be the ideal solution, I think. Update: It has been confirmed by Wolfram Technical Support as a known bug. So if you're reading this in the future, you'll be able to use this answer...

• A pretty good solution, but the y-axis seems not that right. – Kattern Jun 5 '15 at 5:12
• Yes, I think that's actually a bug. I see no reason whatsoever why the y-axis is scaled the way it is. – Jens Jun 5 '15 at 5:14
• Thanks. Do you apply resampling on the list, I get an error like DifferentiatorFilter currently only supports regularly sampled time \ series inputs. Use TimeSeriesResample or TemporalRegularity to make \ the input regularly sampled. – Kattern Jun 5 '15 at 5:29
• I have filed a bug report on this [CASE:3355498]. Let's see if I hear back from them. – Jens Jun 5 '15 at 5:56
• @AlexeyPopkov I got confirmation from Wolfram that this is a known bug. – Jens Jun 9 '15 at 15:48

You can use polynomial regression approach as suggested in this answer:

lm = LinearModelFit[l1, x^Range[0, 6], x];
Plot[Evaluate@lm["Function"][x], {x, Min[l1[[;; , 1]]],
Max[l1[[;; , 1]]]}, PlotStyle -> Red, Prolog -> {Gray, Point[l1]},
PlotLabel -> "Comparison of the model with original data"]
ListPlot[lm["FitResiduals"], PlotRange -> All,
PlotLabel -> "Differences between the original data and the model"]
Plot[Evaluate[lm["Function"]'[x]], {x, Min[l1[[;; , 1]]],
Max[l1[[;; , 1]]]}, PlotLabel -> "First derivative of the model"]


And here is a spline fitting approach based on Hug's idea:

Clear[model, y, x]
nOfControlPoints = 9;
controlPoints =
Subdivide[#1, #2, nOfControlPoints - 1] & @@ MinMax[data[[;; , 1]]];
model[y : {__Real}] :=
Interpolation[Transpose[{controlPoints, y}], Method -> "Spline"]
nlm = NonlinearModelFit[data, model[Array[y, nOfControlPoints]][x],
Array[y, nOfControlPoints], x];
Plot[Evaluate@nlm["Function"][x], {x, Min[data[[;; , 1]]],
Max[data[[;; , 1]]]}, PlotStyle -> Red, Prolog -> {Point[data]},
PlotLabel -> "Comparison of the model with original data"]
ListPlot[nlm["FitResiduals"], PlotRange -> All,
PlotLabel -> "Differences between the original data and the model"]
Plot[Evaluate@nlm["Function"]'[x], {x, Min[data[[;; , 1]]],
Max[data[[;; , 1]]]}, PlotLabel -> "First derivative of the model"]


Another option is to use the Quantile regression with B-splines package by Anton Antonov as described here (or here):

Needs["QuantileRegression"]

qfunc = QuantileRegression[data, 30, {0.5},
InterpolationOrder -> 2][[1]]; // Quiet

Plot[qfunc[x], {x, Min[data[[All, 1]]], Max[data[[All, 1]]]},
Frame -> True, PlotStyle -> Red, Prolog -> Point[data],
PlotLabel -> "Comparison of the model with original data"]
ListPlot[data[[;; , 2]] - qfunc /@ data[[;; , 1]], PlotRange -> All,
Filling -> Axis,
PlotLabel -> "Differences between the original data and the model"]
Plot[qfunc'[x], {x, Min[data[[All, 1]]], Max[data[[All, 1]]]},
Exclusions -> None, PlotLabel -> "First derivative of the model"]


Playing with the number of knots and interpolation order allows to get even better result. For example, with 20 equally-spaced knots and interpolation order 3:

qfunc = QuantileRegression[data, 20, {0.5},
InterpolationOrder -> 3][[1]]; // Quiet
Plot[qfunc'[x], {x, Min[data[[All, 1]]], Max[data[[All, 1]]]},
Exclusions -> None, PlotLabel -> "First derivative of the model"]


• Thanks, but I do not think polynomial regression is a good choice for this job. They are not easy to control, and can be easily overfit. The devertaive calculated has obvious difference with the results posted above. – Kattern Jun 6 '15 at 6:39
• @kattern I updated the answer with another approach. – Alexey Popkov Jun 6 '15 at 7:08
• Thanks a lot! That latter one is definitely a good solution, I am disappoint by the result from BSplineFunction in Mathematica. This spline solution is absolutely the correct way to use spline in this question. – Kattern Jun 6 '15 at 7:41
• @kattern Added another option. – Alexey Popkov Jun 6 '15 at 10:36
• Thanks, I am trying PiecewisePolynomialInterpolation from here, surprisingly, not very hopeful. – Kattern Jun 6 '15 at 10:39

Signal-processing is not one of my strengths, but I think we used to use low-pass filters to remove such noise. It seems to work here:

l2 = LowpassFilter[TimeSeries[l1], .05];
ifn = Interpolation[l2];

plot1 = Plot[Evaluate@Interpolation[l1]'[x], {x, 0, 0.5}, PlotStyle -> {Gray}];
plot2 = Plot[Evaluate@ifn'[x], {x, 0, 0.5}, PlotStyle -> {Red}];

GraphicsRow[{plot2, Show[plot1, plot2]}]


• Yes, MovingAverage is also another kind of low-pass filter. MovingAverage has phase delay in the result, I do not know whether low-pass filter in Mathematica is free of phase delay. – Kattern Jun 5 '15 at 13:30
• Your answer does what I would have expected my answer to do - probably DifferentiatorFilter in my answer invokes LowPassFilter` (+1). – Jens Jun 5 '15 at 18:40
• @Jens I think so, looks like they have pretty close parameters. – Kattern Jun 6 '15 at 8:15