# Is there a way to fit a piecewise linear function to a list?

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 I want to approximate this function by $$8$$ straight lines. Thus the proposed function is defined as 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?

• Would ListInterpolation[sunrises, InterpolationOrder -> 1] work? Also, Extracting the function from InterpolatingFunction object Apr 6, 2022 at 13:03
• @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 Apr 6, 2022 at 13:21
• 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? Apr 6, 2022 at 14:00
• 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 Apr 6, 2022 at 14:05
• 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. Apr 6, 2022 at 14:25

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, data}]},
PlotStyle -> {Red, Thickness[.015]}] 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.

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

x = 1; x = 60; x = 80; x = 110; x = 200; x = 280;
x = 300; x = 320;x = 360;
fit[x_] =
Piecewise[
Transpose[{Fit[
Transpose[{Table[x, {x, x[#[]], x[#[]]}],
sunrises[[x[#[]] ;; x[#[]]]]}], {1, x}, x],
x[#[]] <= x <= x[#[]]}] & /@ Partition[Range[0, 8], 2, 1]]
Show[Plot[fit[x], {x, 1, 350}, PlotStyle -> Red], ListPlot[sunrises]] • The discretization points are fixed here. Is there any way to optimize them too? Apr 6, 2022 at 13:29
• 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. Apr 6, 2022 at 15:47