2
$\begingroup$

I need to find a piecewise linear function that is a proper approximation for sunrise/sunset times of each day of year. Due to some limitations, I cannot use any trigonometric functions. So the sine-looking function has to be approximated by a piecewise linear one. Here is the code to generate data:

location = GeoPosition[{30, 60}];
year = DateRange["2022/3/21", "2023/3/21"];
sunrises = DateDifference[#, Sunrise[location, #], "Minute"] & /@ year;
sunsets = DateDifference[#, Sunset[location, #], "Minute"] & /@ year;

The sunrise data looks like this

sunrise

I want to approximate this function by $8$ straight lines. Thus the proposed function is defined as

piecewise linear

r[x_] := Piecewise[{
 {393 + ((-393 + b)*(-1 + x))/(-1 + x2), x <= x2},
 {b + ((-b + c)*(x - x2))/(-x2 + x3), x2 < x && x <= x3},
 {c + ((-c + d)*(x - x3))/(-x3 + x4), x3 < x && x <= x4},
 {d + ((-d + e)*(x - x4))/(-x4 + x5), x4 < x && x <= x5},
 {e + ((-e + f)*(x - x5))/(-x5 + x6), x5 < x && x <= x6},
 {f + ((-f + g)*(x - x6))/(-x6 + x7), x6 < x && x <= x7},
 {g + ((393 - g)*(x - x7))/(-x7 + x8), x7 < x && x <= x8}},
 393 + ((-393 + i)*(x - x8))/(366 - x8)]

Now the problem is, I cannot find such fit using any functions that I have tried. The closest thing that I have found was This question, whose proposed method failed with this error:

FindFit::nrlnum: The function value {0.,0.0204082,<<364>>} is not a list of real numbers with dimensions {366} at {b,c,d,e,f,g,h,x2,x3,x4,<<4>>} = {370.,350.,<<12>>}.

FWIW, I did set the initial values, as suggested there

init = {{b, 340}, {c, 320}, {d, 340}, {e, 390}, {f, 430}, {g, 450}, {h, 430},
 {x2, 50}, {x3, 80}, {x4, 120}, {x5, 200}, {x6, 270}, {x7, 300}, {x8, 330}};
NonlinearModelFit[First /@ sunrises, r[x], init, x]

Is there any way to do this?

$\endgroup$
7
  • $\begingroup$ Would ListInterpolation[sunrises, InterpolationOrder -> 1] work? Also, Extracting the function from InterpolatingFunction object $\endgroup$
    – Roman
    Apr 6, 2022 at 13:03
  • $\begingroup$ @Roman no, it didn't work. I want the intervals to be exactly calculated. Since the function will be implemented in an embedded device with C $\endgroup$
    – polfosol
    Apr 6, 2022 at 13:21
  • $\begingroup$ Can you afford more multiplications on the embedded CPU? If so, you can try fitting a cubic or a quartic. Now you have one multiplication, one addition and 3 comparisons. A quartic needs 4 multiplications and 4 additions, if you rewrite it into the form a0+x(a1+x(a2+x(a3+a4 x))). Too expensive? $\endgroup$ Apr 6, 2022 at 14:00
  • $\begingroup$ This is not what you asked for, but once you have the Chebyshev series to approximate the function, the function can be evaluated using only arithmetic via Clenshaw's algorithm. See mathematica.stackexchange.com/a/115621/4999 $\endgroup$
    – Michael E2
    Apr 6, 2022 at 14:05
  • 1
    $\begingroup$ You're missing the variable i in init, but it doesn't really fix things, although NonLinearModelFit will give an answer if you include it. $\endgroup$
    – Michael E2
    Apr 6, 2022 at 14:25

2 Answers 2

6
$\begingroup$

You were on the right track.

I fitted the simular function data = Table[-70 Sin[x] + 400 // N, {x, 0, 2 Pi, 2 Pi/365}]; because didn't get sunrise data.

data = Table[-70 Sin[x] + 400 // N, {x, 0, 2 Pi, 2 Pi/365}]; 

ListPlot[data]

r[x_] = Piecewise[{{a + ((-a + b)*(-1 + x))/(-1 + x2), 
0 < x <= x2}, {b + ((-b + c)*(x - x2))/(-x2 + x3), 
x2 < x && x <= x3}, {c + ((-c + d)*(x - x3))/(-x3 + x4), 
x3 < x && x <= x4}, {d + ((-d + e)*(x - x4))/(-x4 + x5), 
x4 < x && x <= x5}, {e + ((-e + f)*(x - x5))/(-x5 + x6), 
x5 < x && x <= x6}, {f + ((-f + g)*(x - x6))/(-x6 + x7), 
x6 < x && x <= x7}, {g + ((h - g)*(x - x7))/(-x7 + x8), 
x7 < x && x <= x8}}, h + ((-h + i)*(x - x8))/(366 - x8)]

init = {{a, 380}, {b, 340}, {c, 320}, {d, 340}, {e, 390}, {f, 
   430}, {g, 450}, {h, 390}, {i, 390}, {x2, 50}, {x3, 80}, {x4, 
   120}, {x5, 200}, {x6, 270}, {x7, 300}, {x8, 330}}

nlfit = NonlinearModelFit[data, r[x], init, x]

nlfit["BestFitParameters"]

(*   {a -> 398.58, b -> 343.842, c -> 327.174, d -> 344.802, e -> 451.384, 
 f -> 470.609, g -> 467.262, h -> 445.643, i -> 400.682, 
 x2 -> 51.9572, x3 -> 91.9915, x4 -> 134.41, x5 -> 229.068, 
 x6 -> 264.087, x7 -> 295.086, x8 -> 326.01}   *)

Plot[nlfit[x], {x, 0, 366}, 
 Epilog -> {PointSize[.005], Point@Transpose[{Range[366], data}]}, 
 PlotStyle -> {Red, Thickness[.015]}]

enter image description here

It is interesting to observe, you do not get absoute symmetry here, although this data are symmetric. Seems Fit doesn't find the global minimum of leastsquares. Think, it depends on initial values you give.

$\endgroup$
1
$\begingroup$

Here is an example with 8 straight lines. It still needs some fine tuning, but this you can do yourself:

x[0] = 1; x[1] = 60; x[2] = 80; x[3] = 110; x[4] = 200; x[5] = 280; 
x[6] = 300; x[7] = 320;x[8] = 360;
fit[x_] = 
 Piecewise[
  Transpose[{Fit[
       Transpose[{Table[x, {x, x[#[[1]]], x[#[[2]]]}], 
         sunrises[[x[#[[1]]] ;; x[#[[2]]]]]}], {1, x}, x], 
      x[#[[1]]] <= x <= x[#[[2]]]}] & /@ Partition[Range[0, 8], 2, 1]]
Show[Plot[fit[x], {x, 1, 350}, PlotStyle -> Red], ListPlot[sunrises]]

enter image description here

$\endgroup$
2
  • 1
    $\begingroup$ The discretization points are fixed here. Is there any way to optimize them too? $\endgroup$
    – polfosol
    Apr 6, 2022 at 13:29
  • $\begingroup$ You could do this, but I think it is not worth the effort, you can do it by hand. Otherwise you would define a function of the points that makes the fits and returns the sum of squared deviations. Then you would use a minimization function to minimize the sum of squares. $\endgroup$ Apr 6, 2022 at 15:47

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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