6
$\begingroup$

Using the following data, which I have used before in other posts (the y-axis is heat flow (mW) and x-axis is Temperature in °C (not time)):

Import["https://pastebin.com/raw/SMKZUtbQ", "Package"]

which plotted using ListLinePlot[data, PlotRange -> {{50, 100}, {-0.1, 1}}] gives:

image

Question: How can I find the value (x and y coordinate) of the inflection(s) point to the left of each peak as shown in the figure (done with power point)?.

For this I have been using the code provided by MarcoB here: Find onset and peak temperatures

and additionally using the following approach:

start = 55;
end = 95;
region = Select[data, start <= #[[1]] <= end &];
fint = Interpolation[region];

which after using something like: infp = {x, fint[x]} /. FindRoot[fint''[x] == 0, {x, 59.4, 61}] finds the value {60.2085, 0.766843} which may be a possible inflection point for the first peak (left side) but it does not find all the inflection points of that peak. So, I am looking for a way to find the inflection(s) point(s) and a way to evaluate if they are indeed inflection(s) point(s).

Here's the point and the plot together using:

Show[Plot[fint[x], {x, start + 0.1, end - 0.1}, 
  PlotRange -> {{start, end}, {-0.5, 2.5}}, PlotStyle -> {Blue}, 
  AspectRatio -> aspect, Frame -> True, FrameStyle -> 14, 
  Axes -> False, GridLines -> Automatic, 
  GridLinesStyle -> Lighter[Gray, .8], 
  FrameTicks -> {Automatic, Automatic}, 
  LabelStyle -> {Black, Bold, 10}], 
 ListPlot[{infp}, PlotStyle -> Red]]

enter image description here

I appreciate your input

EDIT:

When I use infp = {x, fint[x]} /.FindRoot[fint''[x] == 0, {x, 59.8,61}] I find the following values {59.8211, 0.589037} which I check if this is a inflection point using (a very bad code):

belowinfp = fint''[infp[[1]] - 0.001]

aboveinfp = fint''[infp[[1]] + 0.001]

so that if belowinfp and aboveinfp have different signs, then I conclude that it is indeed an inflection point (in this case it is). If I used for instances infp = {x, fint[x]} /.FindRoot[fint''[x] == 0, {x, 59.4,61}] I find the values {60.2085, 0.766843} which using the same test I conclude that this is also an inflection point.

So, I guess the problem reduces to doing this automatically so that it will give me all the values of the inflections points using this test or any other.

$\endgroup$
  • $\begingroup$ Your arrows seem to point to the maxima of each peak, rather than the inflection points. As a start, look for zeroes in the second derivative? I think this goes back to the same problems outline in my answer here, doesn't it? $\endgroup$ – MarcoB Jun 18 at 18:27
  • $\begingroup$ MarcoB! haha yes, I am sorry for the bad drawing but you understand!. It does go back to the answer you provided but I am looking here just for the number which I guess is even much more simpler than your answer in that post. From that answer for instance, how can I extract just that one value (x and y coordinate) of the inflection point? $\endgroup$ – John Jun 18 at 18:33
  • $\begingroup$ I think that using MarcoB’s answer to that other question, you should be able to ask a more specific question than this. For example, he shows how to smoothen the second derivative so that the zeros can be found. So you know how to do this. We don’t need to repeat that here. So what you should ask is, given the second derivative (provide the curve here somehow and link to the answer), how can I find the zeros? That is where you are in your problem solving now. No need to start over. $\endgroup$ – C. E. Jun 18 at 20:41
  • $\begingroup$ If you reframe the question thus, also be prepared for the question ”what have you tried?” The problem: given a curve, find the zeros. Can you find any functions or an approach that might be able to do that? $\endgroup$ – C. E. Jun 18 at 20:49
  • $\begingroup$ @C.E thanks for your suggestion. I edited the question to make it better. At first I was simplying trying something simple like: FindRoot[int''[x] == 0, {x, 56, 60}] (for the first peak) and the same for the other peaks. But this approach does not seem to work very well. $\endgroup$ – John Jun 18 at 22:25
6
$\begingroup$

Differentiating experimental data twice will blow up the noise, so you will probably need to smooth the data to get something usable. @halirutan's answer here applies a GaussianFilter to smooth the data.

To detect the zero crossings, we can use @Daniel Lichtblau's answer here.

The following workflow shows one possible approach that may point you in the right direction.

Import["https://pastebin.com/raw/SMKZUtbQ", "Package"]
start = 55;
end = 95;
region = Select[data, start <= #[[1]] <= end &];
fint = Interpolation[region];
(* Use halirutan's GaussianFilter answer to smooth data *)
ApplyGaussianFilter[data_, r_] := 
  Transpose[{#1, GaussianFilter[#2, r]}] & @@ Transpose[data];
data = ApplyGaussianFilter[data, 2];
(* Use BSplineFunction to Smooth and Resample Data on uniform x scale \
*)
bsf = BSplineFunction[data];
resampleddata = bsf[#] & /@ Subdivide[0, 1, 1000];
(* Create interpolation function *)
ifun = Interpolation[resampleddata, Method -> "Hermite"];
(* Use Daniel Lichtblau's Answer to Find Zeros using NDSolve *)
zeros = Reap[
    NDSolve[{y'[x] == D[ifun''[x], x], 
      WhenEvent[y[x] == 0, Sow[{x, y[x]}]], 
      y[start + 0.1] == ifun''[start + 0.1]}, {}, {x, start + 0.1, 
      end - 0.1}]][[-1, 1]];
pointsOnCurve = {#, ifun[#]} & /@ zeros[[All, 1]];
Plot[{fint[x], ifun[x]}, {x, start + 0.1, end - 0.1}, 
 Epilog -> {PointSize[Medium], Red, 
   Point[pointsOnCurve[[1 ;; -1 ;; 2]]], Green, 
   Point[pointsOnCurve[[2 ;; -1 ;; 2]]]}, PlotRange -> {-0.5, 1}]
Plot[{ifun[x], ifun''[x]}, {x, start + 0.1, end - 0.1}, 
 Epilog -> {PointSize[Medium], Red, Point[zeros[[1 ;; -1 ;; 2]]], 
   Green, Point[zeros[[2 ;; -1 ;; 2]]]}, PlotRange -> {-1, 1}, 
 PlotLabel -> "Smoothed"]
Plot[{fint[x], fint'''[x]}, {x, start + 0.1, end - 0.1}, 
 Epilog -> {PointSize[Medium], Red, Point[zeros[[1 ;; -1 ;; 2]]], 
   Green, Point[zeros[[2 ;; -1 ;; 2]]]}, PlotRange -> {-1, 1}, 
 PlotLabel -> "Unsmoothed"]

Inflection Point Detection

It did a pretty good job at detecting inflection points. Without smoothing, you get lots of false detections.

| improve this answer | |
$\endgroup$
  • $\begingroup$ Tim Laska thank you! This is a very good answer and very adequate to what I needed !. I appreciate very much your help ! $\endgroup$ – John Jun 19 at 18:18
  • $\begingroup$ @John You are welcome. I left out the definition of ApplyGaussianFilter, but I edited the post to include it. $\endgroup$ – Tim Laska Jun 19 at 18:51
  • $\begingroup$ Tim Laska! Thank you! I accepted your answer!. One quick question: Why is it necessary to resample the data and use a second interpolation function with BSplineFunction? $\endgroup$ – John Jun 19 at 20:46
  • 1
    $\begingroup$ @John You are welcome. It may not be necessary. Initially, I was using ListConvolve to smooth the data and I think that works better on equally spaced points. BSplineFunction will provide some additional smoothing of the data and it does require the curve to pass through every point like Interpolation does. $\endgroup$ – Tim Laska Jun 19 at 21:41
5
$\begingroup$

This is a testbed case for the function RegularisedInterpolation !

Import["https://pastebin.com/raw/SMKZUtbQ", "Package"]
fit = RegularisedInterpolation[data,
         FitRegularization->{"Curvature", 0.1}]

picture of output RegularisedInterpolatingFunction

Show[
  ListPlot[data, PlotRange -> {{50, 100}, Automatic}],
  Plot[fit[x], {x, 50, 100},PlotStyle-> Directive[Red, Dashed]]
]

plot of data and interpolation

Thanks to the regularisation it can be differentiated twice.

d2fit[x_] = D[fit[x], x, x];
Plot[d2fit[x], {x, 60, 70}]

plot of the second derivative

Then you can bracket the zeros:

FindRoot[d2fit[x] == 0, {x, 62, 64}]
FindRoot[d2fit[x] == 0, {x, 64, 66}]

(* 
 {x->62.3478}
 {x->64.4095} 
*) 

or use Daniel Lichtblau's zero crossings.

Validation

We can check that the result is fairly robust to the strength of smoothing

Table[
  fit = RegularisedInterpolation[data, 
    FitRegularization -> {"Curvature", 10^i}];
  d2fit[x_] = D[fit[x], x, x];
  x /. {FindRoot[d2fit[x] == 0, {x, 62, 64}],
  FindRoot[d2fit[x] == 0, {x, 64, 66}]},
  {i, -3, 1}]

(* {
{62.227,  64.4562},
{62.289,  64.4582},
{62.3478, 64.4095},
{62.3464, 64.413},
{62.2796, 64.4675}
} *)
| improve this answer | |
$\endgroup$
  • $\begingroup$ Very nice. I was not aware of this function! (+1) $\endgroup$ – MarcoB Jun 19 at 16:08
  • $\begingroup$ chris thank you very much! This is an amazing answer!. May I ask that for some reason in my Mathematica (version 12.1) the function RegularisedInterpolationdoes not work and I cannot get information about it in the regular Mathematica documentation. Is there anything else I need to do to use that function?. Thanks again $\endgroup$ – John Jun 19 at 16:52
  • 2
    $\begingroup$ @John its because I wrote it. Follow the link? $\endgroup$ – chris Jun 19 at 16:53
  • $\begingroup$ chris thank you I was not aware of that. Now I copied and pasted the link and it worked! and I was able to have the same first plot you showed. However, the line Plot[d2fit[x_] = D[fit[x], x, x] // Evaluate, {x, 60, 70}] does not run for me. I think I am simply still too naive on how things here work or if I still need to add something extra to run that line as well. $\endgroup$ – John Jun 19 at 17:15
  • $\begingroup$ @John sorry my bad. It should work with d2fit defined before the plot. $\endgroup$ – chris Jun 20 at 6:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.