How can we fit a closed curve?

Question

Data={{15., -9.8861}, {14., -9.80713}, {13., -9.60669}, {12., -9.37226}, \
{11., -9.56058}, {10., -9.52826}, {9., -9.3295}, {8., -9.16436}, {7., \
-9.04076}, {6., -8.89642}, {5., -9.46254}, {4., -9.5109}, {3., \
-8.90791}, {2., -8.09209}, {1., -8.54887}, {0., -7.93869}, {0., \
-7.93869}, {-1., -6.92529}, {-2., -6.48322}, {-3., -5.36115}, {-4., \
-2.29101}, {-5., 1.52439}, {-6., 4.79197}, {-7., 7.2295}, {-8., 
  8.29515}, {-9., 9.19184}, {-10., 9.39874}, {-11., 9.56396}, {-12., 
  9.88203}, {-13., 9.62128}, {-14., 9.84259}, {-15., 9.89896}, {-15., 
  9.89896}, {-14., 9.81998}, {-13., 9.61955}, {-12., 9.38511}, {-11., 
  9.57344}, {-10., 9.54112}, {-9., 9.34235}, {-8., 9.17722}, {-7., 
  9.05361}, {-6., 8.90928}, {-5., 9.47539}, {-4., 9.52376}, {-3., 
  8.92076}, {-2., 8.10495}, {-1., 8.56173}, {0., 7.95154}, {0., 
  7.95154}, {1., 6.93814}, {2., 6.49608}, {3., 5.374}, {4., 
  2.30386}, {5., -1.51154}, {6., -4.77912}, {7., -7.21664}, {8., \
-8.28229}, {9., -9.17899}, {10., -9.38589}, {11., -9.55111}, {12., \
-9.86918}, {13., -9.60843}, {14., -9.82973}, {15., -9.8861}};    

ListPlot[Data]

This is my try

fitFunc = 
  Total@Table[
    a[i] Tanh[b[i] x + c[i]] + d[i] Tanh[e[i] x - f[i]], {i, 1, 1}];
params = Flatten@
   Table[{a[i], b[i], c[i], d[i], e[i], f[i]}, {i, 1, 1}];

nlm1 = NonlinearModelFit[Data, fitFunc, params, x];
Plot[{nlm1[t]}, {t, -15, 15}]

Define what you mean by 'closed shape profile'? Your code only has a curve, not a shape with area. — flinty, Commented Dec 20, 2022 at 14:21
This looks like a classical hysteresis curve. I would separately fit its upper and lower parts. — Alexei Boulbitch, Commented Dec 20, 2022 at 15:12
Consider ListCurvePathPlot and FindCurvePath. There is also the resource function CurveToBSpline. I would test them if you copy pasted the list of points in the question. I do not want to access a Dropbox account. — userrandrand, Commented Dec 20, 2022 at 19:10

azerbajdzan · Accepted Answer · 2022-12-20 22:13:31Z

Just do the same you did but twice for each half of the data.

dat = {{15., -9.8861}, {14., -9.80713}, {13., -9.60669}, {12., \
-9.37226}, {11., -9.56058}, {10., -9.52826}, {9., -9.3295}, {8., \
-9.16436}, {7., -9.04076}, {6., -8.89642}, {5., -9.46254}, {4., \
-9.5109}, {3., -8.90791}, {2., -8.09209}, {1., -8.54887}, {0., \
-7.93869}, {0., -7.93869}, {-1., -6.92529}, {-2., -6.48322}, {-3., \
-5.36115}, {-4., -2.29101}, {-5., 1.52439}, {-6., 4.79197}, {-7., 
    7.2295}, {-8., 8.29515}, {-9., 9.19184}, {-10., 9.39874}, {-11., 
    9.56396}, {-12., 9.88203}, {-13., 9.62128}, {-14., 
    9.84259}, {-15., 9.89896}, {-15., 9.89896}, {-14., 
    9.81998}, {-13., 9.61955}, {-12., 9.38511}, {-11., 
    9.57344}, {-10., 9.54112}, {-9., 9.34235}, {-8., 9.17722}, {-7., 
    9.05361}, {-6., 8.90928}, {-5., 9.47539}, {-4., 9.52376}, {-3., 
    8.92076}, {-2., 8.10495}, {-1., 8.56173}, {0., 7.95154}, {0., 
    7.95154}, {1., 6.93814}, {2., 6.49608}, {3., 5.374}, {4., 
    2.30386}, {5., -1.51154}, {6., -4.77912}, {7., -7.21664}, {8., \
-8.28229}, {9., -9.17899}, {10., -9.38589}, {11., -9.55111}, {12., \
-9.86918}, {13., -9.60843}, {14., -9.82973}, {15., -9.8861}};
Data = dat[[1 ;; 32]];
fitFunc = 
  Total@Table[
    a[i] Tanh[b[i] x + c[i]] + d[i] Tanh[e[i] x - f[i]], {i, 1, 1}];
params = Flatten@
   Table[{a[i], b[i], c[i], d[i], e[i], f[i]}, {i, 1, 1}];

nlm1 = NonlinearModelFit[Data, fitFunc, params, x];
Plot[{nlm1[t]}, {t, -15, 15}];

Data = dat[[33 ;; -1]];
fitFunc = 
  Total@Table[
    a[i] Tanh[b[i] x + c[i]] + d[i] Tanh[e[i] x - f[i]], {i, 1, 1}];
params = Flatten@
   Table[{a[i], b[i], c[i], d[i], e[i], f[i]}, {i, 1, 1}];

nlm1 = NonlinearModelFit[Data, fitFunc, params, x];
Plot[{nlm1[t]}, {t, -15, 15}];
Show[%%%%%%, %]

$\begingroup$ Nice I did not know the data was already sorted $\endgroup$
– userrandrand
Commented Dec 20, 2022 at 22:14 — userrandrand, Commented Dec 20, 2022 at 22:14

demm · Accepted Answer · 2022-12-23 11:10:42Z

We split the data set:

joinDS = Module[{a = SplitBy[Sort[Data], First]}, {a[[All, 1]], a[[All,2]]}];

We define the function:

branch[x_, a_, b_, c_] := a Tanh[b x + c]

and we fit the data using the MultiNonlinearModelFit

    fit = ResourceFunction["MultiNonlinearModelFit"][joinDS, {branch[x, a1, b1, c1], 
branch[x, a2, b2, c2]}, {a1, b1, c1,a2, b2, c2}, {x}, Method -> "NMinimize"]

(you can see that the output involves the Switch function)

Here are the results:

Show[{ListPlot[Data], Plot[{fit[1, x], fit[2, x]}, {x, -15, 15}]}]

However, since we are talking about hysteresis I would go for these functions:

B1[h_, b_, s_, x_] := -s  Tanh[1/h ArcTanh[b/s] (x + h)];
B2[h_, b_, s_, x_] := -s  Tanh[1/h ArcTanh[b/s] (x - h)];

and I would fit it with the same fitting parameters (h,b,s):

  fit2 = ResourceFunction["MultiNonlinearModelFit"][joinDS, 
<|"Expressions" -> {B2[h, b, s, x], B1[h, b, s, x]}, "Constraints" -> h > 4 && b > 8 && s > 9|>, {h, b, s}, {x}]

Which gives:

Show[{ListPlot[Data], Plot[{fit2[1, x], fit2[2, x]}, {x, -15, 15}]}]

This code is much simpler than mine. However, note that the parameters were not properly fit. In particular, with all parameters free, the fitted parameters to the upper branch are different from that of the lower branch, despite the data being the same. These are not real data, the points that make the lower branch are the reflection about the origin of those on the upper branch. — Vito Vanin, Commented Dec 23, 2022 at 12:34

userrandrand · Accepted Answer · 2022-12-20 22:11:29Z

The curve is difficult because of the noisy ends of the hysteresis. To see this consider ordering the points using FindCurvePath:

Note : data is Data in your list of points.

sorted = data[[#]] & /@ FindCurvePath[data];

We may plot these points:

The different colors represent different sections that FindCurvePath found. Notice multiple sections at the end suggesting that the noise there makes it difficult to find a path. You can remove those unwanted segmentations and study them separately. I will only focus on the large part to show some methods.

Note: You can directly use ListLinePlot[sorted[2]] and changed the InterpolationOrder option but because of the noise, interpolation should be avoided.

Method 1 : Find Splines

Here we use the resource function CurveToBSplineFunction :

spline = ResourceFunction["CurveToBSplineFunction"];

We can choose parameters using Manipulate although it can be a bit slow:

Manipulate[
 Show[ParametricPlot[
   spline[sorted[[2]], d, s, "CurveClosed" -> True][t], {t, 0, 1}], 
  ListPlot[data, PlotStyle -> Red]], {d, 1, 10, 2}, {s, 0.05, 1, 0.1}]

degree d=1 spline and scale s=0.05:

degree d=9 spline and scale s=0.85

Method 2: Using machine learning separately on the ordered x and y coordinates

Note: There is an issue at the boundary so you might not want to use this method although it has the advantage of not requiring parameter adjustments

Recalling that sorted[[2]] is the ordered data points in the middle, we can introduce a fictitious time with discrete values {1,2,...Length[sorted[[2]]]} and fit {x[t],t} and {y[t],t} and use a ParametricPlot. I imagine that the first method above does something similar internally.

With n the sized of the considered data (sorted[[2]]) we may create a list of data points {{x[1],1},{x[2],2}....{x[n],n}},{{y[1],1},{y[2],2},...,{y[n],n}} :

Note: ⎵=[UnderBracket]

parametric⎵data = 
  sorted[[2]] // Transpose // Map[Transpose[{Range@Length@#, #}] &];

In my limited experience fitting noisy data, Predict with the method GaussianProcess has given the nicest results. I learned about from an answer I found on stack exchange you may search for other examples if needed.

predict = 
  Predict[Rule @@@ #, Method -> "GaussianProcess", 
     PerformanceGoal -> "Quality"] & /@ parametric⎵data;

Plot and compare the fit with the data:

Show[ParametricPlot[Through[predict[x]], {x, 0, 35}], 
 ListPlot[sorted[[2]], PlotStyle -> Red]]

Notice that the curve is not closed I do not know why.

JimB · Accepted Answer · 2022-12-23 23:37:46Z

This is just an extended comment.

It appears that you really only have 32 data points rather than 64 data points so a description of how you constructed the data should be included in your question. Also, as I mentioned in a comment, if this is really a hysteresis loop, the third (but missing) variable that predicts both the "x" and "y" coordinate should be included.

Why do I think you only have 32 data points? Consider the first 32 responses and the negative of the last 32 responses:

d1 = Data[[1 ;; 32, 2]];
d2 = -Data[[33 ;; 64, 2]];
ListPlot[{d1, d2}, PlotStyle -> {PointSize[0.02], {PointSize[0.0075], White}}]

The differences in the two lists are essentially a constant:

d1 - d2
(* {0.01286, 0.01285, 0.01286, 0.01285, 0.01286, 0.01286, 0.01285, 
    0.01286, 0.01285, 0.01286, 0.01285, 0.01286, 0.01285, 0.01286, 
    0.01286, 0.01285, 0.01285, 0.01285, 0.01286, 0.01285, 0.01285, 
    0.01285, 0.01285, 0.01286, 0.01286, 0.01285, 0.01285, 0.01285, 
    0.01285, 0.01285, 0.01286, 0.01286} *)

That perfect fit doesn't happen with 64 real-life data points (even in physics). So I think you really only have 32 data points.

Vito Vanin · Accepted Answer · 2022-12-21 17:52:28Z

In order to warrant a closed curve, the parameters used in the two branches of the closed curve must be constrained. I suggest the inclusion of another predictive variable to distinguish the two branches, and fit the parameters to all data simultaneously. First, include the identification of the branch in the data:

dataC =
Join[Transpose[{data[[;;32,1]],ConstantArray["down",32], data[[;;32,2]]}],
Transpose[{data[[33;;,1]],ConstantArray["up",32], data[[33 ;; , 2]]} ] ]

Define the function

branch[x_, a_, b_, c_] := a Tanh[b  x + c]

fitFunc[x_, a_, b_, c_, d_, e_, f_, direction_] :=
    Switch[direction, "up", branch[x, a, b, c], "down", branch[x, d, e, f] ]

Fit:

result =  NonlinearModelFit[dataC, fitFunc[t, a, b, c, a, e, f, dir], 
 {{a, -10}, {b, 0.3}, {c, -1}, {e, 0.3}, {f, 1.2}}, {t, dir} ]

Notice that a is repeated in fitFunc[ ] call, which imposes the constraint of equal amplitudes for the two branches of the curve. Plot:

{af, bf, cf, ef, ff} = result["ParameterTableEntries"][[All, 1]]

Show[ ListPlot[{data[[33 ;;]], data[[;; 32]]}, PlotLegends ->{"up","down"}], 
      Plot[{fitFunc[x, af, bf, cf, af, ef, ff, "up"], 
         fitFunc[x, af, bf, cf, af, ef, ff, "down"]}, {x, -15, 15}, 
         PlotLegends -> {"up", "down"}]  ]

With real data, the parameters b and e will be different, as well as c and f. However, the constraint b == e can be imposed just changing e by b when calling fitFunc[ ]. If needed, the hypothesis that the curve is closed can be tested, just change one of the a's in fitFunc[ ] call by d, and providing an estimate to NonlinearModelFit.

Equivalent approach as the one presented by @Vito Vanin (i.e. using the Switch function) can be achieved using MultiNonlinearModelFit (resources.wolframcloud.com/FunctionRepository/resources/…) which uses the Switch command under the bonnet. — demm, Commented Dec 22, 2022 at 9:12
@demm You should post an answer based on MultiNonlinearModelFit[ ], which avoids all the chore of recasting the data, writing the Switch[ ], and constraining the parameters. It is the best answer for this question. From now on, I'll use it for constrained parameters fit. — Vito Vanin, Commented Dec 22, 2022 at 18:57

Stack Exchange Network

How can we fit a closed curve?

5 Answers 5

Your Answer

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

Hot Network Questions

How can we fit a closed curve?

5 Answers 5

Your Answer

Sign up or log in

Post as a guest

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

Related

Hot Network Questions

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