# Finding a polynomial representation for a sum function

I'm trying to find a polynomial representation for this horrendous function:

$$f(x)=\frac{\frac{1}{2}-\frac{4}{1+10,000x}\sum_{n=0}^{\infty}\frac{(-1)^n}{\big[\big(n+\frac{1}{2}\big)\pi\big]^4}\tanh\Bigg\{\frac{\big[\big(n+\frac{1}{2}\big)\pi\big](1+10,000x)}{2}\Bigg\}}{\frac{1}{3}-\frac{4}{1+10,000x}\sum_{n=0}^{\infty}\frac{1}{\big[\big(n+\frac{1}{2}\big)\pi\big]^5}\tanh\Bigg\{\frac{\big[\big(n+\frac{1}{2}\big)\pi\big] (1+10,000x)}{2}\Bigg\}}$$

I managed to plot it quite easily with this code:

Cs[x_] := 1/2 - (4/(1 + 10000*x))*Sum[((-1)^n/((n + 0.5)*Pi)^4)*Tanh[((n + 0.5)*(Pi/2))*(1 + 10000*x)], {n, 0, 100}];
Cp[x_] := 1/3 - (4/(1 + 10000*x))*Sum[(1/((n + 0.5)*Pi)^5)*Tanh[((n + 0.5)*(Pi/2))*(1 + 10000*x)], {n, 0, 100}];
Plot[Cs[x]/Cp[x], {x, 0, 1}]


The graph

Now I want to integrate this function, as part of an ODE. Since Wolfram is struggling with integrating a sum function, I want to find a polynomial for Wolfram to work with. I have tried different curve fitting function like FindFit and InterpolatingPolynomial, but haven't managed to get them to work and give me the coefficients of the polynomial.

Is there a simple way to extract a polynomial fit from the graph, perhaps?

Or a different way to get it to work?

Any help would be appreciated!

• There is a syntax error in your code; would you please correct it? Jul 22, 2017 at 12:54
• Sure, if you could just point me in the right direction. I can't find it. Jul 22, 2017 at 12:58
• I believe I corrected it. Jul 22, 2017 at 13:00
• Yes, you have. Thanks. Jul 22, 2017 at 13:01

## Fitting a polynomial: choose a good basis

Let's look at the function:

Plot[Cs[x]/Cp[x], {x, 0, 1}, PlotRange -> All]


It doesn't look much like a polynomial, except perhaps like the flipped-around graph of x^1000. Let's try something like that:

With[{r = 1.02},
lmf = LinearModelFit[
SetPrecision[Table[{r^p, Cs[r^p]/Cp[r^p]}, {p, -800, 0}],
MachinePrecision],
Table[(1 - x)^(200 n), {n, 26}],
x
]
]


Error is not bad, but one could hope for better:

Plot[Cs[x]/Cp[x] - lmf[x] // Evaluate, {x, 0, 1}, PlotRange -> All]


The coefficients are given by

lmf@"BestFitParameters"


The coefficients are quite frightening numerically, many between 10^6 and upwards of 10^8 and alternating signs. Note that the basis functions (1 - x)^k should never be expanded, since k ranges up to over 5000.

## Chebyshev approximation

According to Weierstrass, there should be no limit to how well we can approximate (at least with exact coefficients and inputs). Chebyshev interpolations typically are robust, but we will see that does not make it easy to approximate the OP's funcion. The function chebSeries below computes the degree n series of Chebyshev coefficients of the Chebyshev polynomials of the interpolation through the Chebyshev points (x in the code):

chebSeries[f_, a_, b_, n_, prec_ : MachinePrecision] :=
Module[{x, y, cc},
x = Rescale[Sin[N[(π Range[n, -n, -2])/(2 n), MachinePrecision]], {-1, 1}, {a, b}];
y = f /@ x;
cc = Sqrt[2/n] FourierDCT[y, 1];
cc[[{1, -1}]] /= 2; cc]


Let's do an interpolation of degree 2048, which is not a whole lot more than is needed in this case:

cc = chebSeries[x \[Function] Cs[x]/Cp[x], 0, 1, 2^11];


As the order of the interpolation increases, the coefficients of the interpolation approach the (infinite) Chebyshev series expansion of the function. The error of a truncated series is at most the sum of the absolute values of the lopped-off coefficients. For a nice function, these coefficients eventually converge exponentially to zero, which is the case here:

ListPlot[RealExponent@cc,
GridLines -> {None, Log10[Max[Abs@cc] $MachineEpsilon] {1/2, 1}}]  We can get an numerical estimate of the error of truncation from the coefficients themselves: cc[[510 ;; 513]] cc[[-4 ;;]] (* {-9.30101*10^-8, 9.18019*10^-8, -9.05997*10^-8, 8.94039*10^-8} {-5.89806*10^-17, 5.26922*10^-17, -5.37764*10^-17, 1.11022*10^-16} *)  The corresponding polynomial (truncated series of degree n) is given by With[{coeffs = cc[[ ;; n+1]]}, coeffs.ChebyshevT[Range[0, Length@ coeffs - 1], 1 - 2 x] ]  We can get almost machine precision this way. Let's begin by estimating the error by summing the tail, for both single and double precision accuracy, and calculate the number of terms of the Chebyshev series needed to obtain the desired precisions. single = Length@cc - Module[{sum = 0.}, LengthWhile[Reverse@cc, (sum += #) < Sqrt@$MachineEpsilon/2 &]]
double = Length@cc -
Module[{sum = 0.}, LengthWhile[Reverse@cc, (sum += #) < $MachineEpsilon &]] (* 617 1535 *)  The function cheval evaluates a given Chebyshev series with coefficients cc over the interval {a, b}. The argument to ChebyshevT needs to be rescaled to run from -1 to 1 as x runs from a to b. It is also important that ChebyshevT[] not evaluate on symbolic x, because it will evaluate to a polynomial in powers of x, which will be numerically unstable to evaluate, especially at such extremely high degrees as we have. Simplify@Rescale[x, {a, b}, {-1, 1}] (* (a + b - 2 x)/(a - b) *) ClearAll[cheval]; cheval[x_?NumericQ, a_?NumericQ, b_?NumericQ, cc_] := cc.ChebyshevT[Range[0, Length@cc - 1], (a + b - 2 x)/(a - b)];  Here is a plot of the logarithm of the absolute error of the "single" and "double" precision series. We see that it is possible to achieve single precision accuracy (when computing with double precision floats), but the double precision series sometimes has small but noticeable rounding errors. Plot[Cs[x]/Cp[x] - {cheval[x, 0, 1, cc[[;; single]]], cheval[x, 0, 1, cc[[;; double]]]} // RealExponent // Evaluate, {x, 0, 1}, PlotLabel -> "Log Absolute Error", MaxRecursion -> 2, PlotLegends -> {Row[{"Order ", single - 1}], Row[{"Order ", double - 1}]}, PlotRange -> {0, -19}, GridLines -> {None, Log10[Max[Abs@cc]$MachineEpsilon] {1/2, 1}}]


The OP stated than a polynomial approximation was desired. Chebyshev series are probably the best way to do it easily, but their degrees might be prohibitive.

## Piecewise Chebyshev approximation

Finally, piecewise interpolation will make it easier to have small degrees, and in Mathematica the inconvenience of piecewise functions is small.

The following divides the interval {0, 1} at 2^n for n = -1, -2,..., -14 and forms the degree-16 Chebyshev interpolation over each subinterval.

pwf = Piecewise@
Join[{{cheval[x, 0, 2^(-14),
chebSeries[x \[Function] Cs[x]/Cp[x], {0, 2^(-14)}, 16]],
0 <= x <= 2^(-14)}},
Table[{cheval[x, 2^(n - 1), 2^(n),
chebSeries[x \[Function] Cs[x]/Cp[x], {2^(n - 1), 2^(n)}, 16]],
2^(n - 1) <= x <= 2^(n)}, {n, -13, 0}]
];


Here is a plot of the logarithm of the absolute error. We quite easily get near machine precision accuracy.

Plot[Cs[x]/Cp[x] - pwf // RealExponent // Evaluate,
{x, 0, 1},
PlotLabel -> "Log Absolute Error", PlotRange -> {-13, -17},
GridLines -> {None, Log10[Max[Abs@cc] \$MachineEpsilon] {1/2, 1}}]


One could also use chebInterpolation to construct an InterpolatingFunction instead of a Piecewise function. Both can be used inside NDSolve.

• Wow, bravo man! The Chebyshev approximation gives marvelous results. Even the LinearFitModel yielded a decent result. Thanks a lot for your time, this was very helpful! Jul 24, 2017 at 13:55
• How would I extend the interval of interpolation from [0,1] to [0,10], for the Chebyshev approximation? Jul 25, 2017 at 13:02
• cc = chebSeries[x \[Function] Cs[x]/Cp[x], 0, 10, 2^11];? You might need to increase the 2^11 to 2^12 (the Fourier DCT is more efficient if you use a power of two). Then evaluate with cheval[x, 0, 10, cc] (or use something like cc[[;; double]] for cc to reduce the degree). That is, in the code the variables {a, b} represent the interval. Jul 25, 2017 at 13:55
• @VortexSheet I forgot to do the @ thing on my mobile, but I've been without M & a computer for a week. Jul 31, 2017 at 15:45
Cs[x_] :=
1/2 - (4/(1 + 10000*x))*
Sum[((-1)^n/((n + 1/2)*Pi)^4)*
Tanh[((n + 1/2)*(Pi/2))*(1 + 10000*x)], {n, 0, 100}];
Cp[x_] :=
1/3 - (4/(1 + 10000*x))*
Sum[(1/((n + 1/2)*Pi)^5)*
Tanh[((n + 1/2)*(Pi/2))*(1 + 10000*x)], {n, 0, 100}];

pts = Table[{x, Cs[x]/Cp[x]}, {x, 0, 1, 0.01}];

Manipulate[
intp = Interpolation[pts, InterpolationOrder -> order];
{Plot[{Cs[x]/Cp[x], intp[x]}, {x, 0, 1.},
PlotStyle ->
{{Blue, Thick}, {Red, AbsoluteDashing[{5, 10}]}},
PlotLegends -> Placed["Expressions", {0.8, 0.25}],
ImageSize -> 360], "",
StringForm["Integral = ",
NIntegrate[intp[x], {x, 0, 1}]]} //
Column,
{{order, 1, "InterpolationOrder"}, Range[1, 9, 2]},
SynchronousUpdating -> False]


EDIT:

For InterpolationOrder->1 the interpolation is equivalent to a Piecewise function consisting of a straight line drawn between adjoining pairs of points

Clear[line];
line[{{x1_, y1_}, {x2_, y2_}}] :=
Evaluate[a*x + b /. Solve[{y1 == a*x1 + b, y2 == a*x2 + b}, {a, b}][[1]]]

intp1[x_] =
Piecewise[{line@#, #[[1, 1]] <= x < #[[2, 1]]} & /@ Partition[pts, 2, 1]]


Piecewise functions can be integrated directly

Integrate[intp1[x], {x, 0, 1}]

(*  1.49914  *)

• Excellent, thank you. Is there any way to extract the coefficients of this polynomial? Jul 23, 2017 at 11:36
• Then it's a good method for plotting purposes, but not so much for simplifying unpleasant functions to polynomials. Thank you for your effort. Jul 24, 2017 at 11:36

I'm not quite sure what you're looking for, and you said you already tried InterpolatingPolynomial, but maybe you struggled with its input form, so I offer:

InterpolatingPolynomial[
Table[{x, Cs[x]/Cp[x]}, {x, 0, 1, 0.01}],
x
] // HornerForm

1.3337 + x (-3.4577*10^13 +
x (1.78876*10^16 +
x (-4.34422*10^18 +
x (6.655*10^20 +
x (-7.2754*10^22 +
x (6.08033*10^24 +
x (-4.06021*10^26 +
x (2.23423*10^28 +
x (-1.03653*10^30 +
x (4.12621*10^31 +
x (-1.42925*10^33 +
x (4.35704*10^34 +
x (-1.18005*10^36 +
x (2.86216*10^37 +
x (-6.25917*10^38 +
x (1.24143*10^40 +
x (-2.24454*10^41 +
x (3.7161*10^42 +
x (-5.65629*10^43 +
x (7.94324*10^44 +
x (-1.03244*10^46 +
x (1.24559*10^47 +
x (-1.39845*10^48 +
x (1.46453*10^49 +
x (-1.43368*10^50 +
x (1.31449*10^51 +
x (-1.13079*10^52 + x (9.14192*10^52 +
x (-6.95619*10^53 + x (4.98866*10^54 +
x (-3.37618*10^55 + x (2.15876*10^56 +
x (-1.30554*10^57 + x (7.47502*10^57 +
x (-4.05578*10^58 + x (2.0871*10^59 +
x (-1.01943*10^60 + x (4.72968*10^60 +
x (-2.08571*10^61 + x (8.74757*10^61 + x (-3.49121*10^62 +
x (1.3266*10^63 + x (-4.80151*10^63 + x (1.65605*10^64 +
x (-5.44493*10^64 + x (1.70717*10^65 + x (-5.1057*10^65 +
x (1.45694*10^66 + x (-3.9676*10^66 + x (1.03133*10^67 +
x (-2.55922*10^67 + x (6.06335*10^67 + x (-1.37165*10^68 +
x (2.96287*10^68 + x (-6.11121*10^68 + x (1.20356*10^69 +
x (-2.26313*10^69 + x (4.06259*10^69 + x (-6.96121*10^69 +
x (1.13835*10^70 + x (-1.77614*10^70 + x (2.64347*10^70 +
x (-3.75173*10^70 + x (5.07569*10^70 + x (-6.54318*10^70 +
x (8.03368*10^70 + x (-9.38969*10^70 + x (1.04412*10^71 +
x (-1.10394*10^71 + x (1.10898*10^71 + x (-1.0577*10^71 +
x (9.56945*10^70 + x (-8.20533*10^70 + x (6.66106*10^70 +
x (-5.11376*10^70 + x (3.70807*10^70 + x (-2.53613*10^70 +
x (1.63363*10^70 + x (-9.89394*10^69 + x (5.6236*10^69 +
x (-2.99361*10^69 + x (1.48906*10^69 + x (-6.90316*10^68 +
x (2.97403*10^68 + x (-1.18681*10^68 + x (4.37059*10^67 +
x (-1.47898*10^67 + x (4.57624*10^66 + x (-1.28728*10^66 +
x (3.26969*10^65 + x (-7.43831*10^64 + x (1.50067*10^64 +
x (-2.65224*10^63 + x (4.04283*10^62 + x (-5.20733*10^61 +
x (5.51104*10^60 + x (-4.60123*10^59 + x (2.84166*10^58 +
x (-1.15408*10^57 +
2.31199*10^55 \
x)))))))))))))))))))))))))))))))))))))))))))))))))))))))))))))))))))))))))))))))))))\
))))))))))))))))


The method above has serious oscillation problems. VortexSheet proposed nonuniform sampling with more points at the ends. This helps with the documentation example but I was not able to solve the case above. For what it's worth here is the code I applied to the documentation example:

f[x_] := 1/(1 + 25 x^2);

ptfn[n_Integer?Positive] := -# ⋃ # & @ Array[N@Log[n + 1, #] &, n + 1]

points = Table[{x, f[x]}, {x, -1, 1, .1}];

points2 = {#, f@#} & /@ ptfn[15];

Plot[Evaluate[InterpolatingPolynomial[#, x]], {x, -1, 1}, PlotRange -> All,
Epilog -> {Red, Map[Point, #]}] & /@ {points, points2} // Column


• Excellent, this will work. I guess I'm just a noob in Mathematica. Jul 22, 2017 at 13:17
• When I try to plot this polynomial, I get an oscillating graph that does not behave similarly to the original sum function... What's going on? Jul 22, 2017 at 13:41
• @VortexSheet Okay, I didn't check the output of this before posting, as I didn't know if you had tried this syntax or not. Looking at not now I see that it is not good, but also that you Accepted this answer, which confuses me. Were you able to make this work for you, and if so, how? If not why did you Accept it? Jul 23, 2017 at 13:08
• Your answer helped me get an understanding, but you're right, I shouldn't have accepted it as the polynomial is wrong. What's wrong in the code then, do you know? Jul 23, 2017 at 13:14
• @VortexSheet I believe it is the problem described under the Possible issues section of the InterpolatingPolynomial documentation. I don't know how to create a single polynomial expression that fits your function properly; as mentioned in that section "Interpolation uses a lower-order piecewise polynomial that does not have this problem." Jul 23, 2017 at 13:19