# Find good polynomial fit of square root

I would like to find a good polynomial fit of the square root. I tried the following:

m = Fit[Table[Sqrt[x], {x, 0, Pi}], {x, x^2}, x]
Plot[{Sqrt[x], m}, {x, 0, Pi}]


But the results are not convincing:

Any help would be much appreciated !

• Square root is a classic example of a function that polynomials don't fit well. It has a branch point at zero, something no polynomial can reproduce. Sep 3, 2018 at 22:12
• @JohnDoty is quite right. The gradient at zero is infinite, which cannot be matched with a polynomial. Can you explain why you want to do this? A polynomial approximation of Sqrt[1+x] could work well. Sep 3, 2018 at 22:21
• @JohnDoty I agreee, but there ought to be a better solution than mine, right ? What other non polynomial could be used ? Sep 3, 2018 at 22:25
• @james try m = Fit[Table[{x, Sqrt[x]}, {x, 0, Pi, .01}], {1, x, x^2}, x] Plot[{Sqrt[x], m}, {x, 0, Pi}] Sep 3, 2018 at 22:28
• @RohitNamjoshi This is actually pretty good ! Thanks ! Sep 3, 2018 at 22:32

Here's a Gauss-Legendre interpolation, adapted from FunctionInterpolation over an open interval:

ClearAll[gaussLegendreInterpolation];
(*construct barycentric interpolation at the Gauss-Legendre nodes*)
gaussLegendreInterpolation[deg_Integer?Positive, f_, {a_, b_}] :=
Module[{xj, fj, lj, wj},
{xj, wj} = Most@NIntegrateGaussRuleData[deg + 1, 18 + 1.5 deg];
wj = 2 wj;
lj = (-1)^Range[0, deg] Sqrt[(1 - Rescale[xj, {0, 1}, {-1, 1}]^2) wj];
xj = Rescale[xj, {0, 1}, {a, b}];
fj = f /@ xj;
StatisticsLibraryBarycentricInterpolation[xj, fj, Weights -> lj]];


Degree 2 interpolant:

sqrt = gaussLegendreInterpolation[2, Sqrt, {0, Pi}][x] // Simplify // N
(*  0.358018 + 0.698352 x - 0.0817347 x^2  *)

Plot[{Sqrt[x], sqrt}, {x, 0, Pi}, WorkingPrecision -> Precision@sqrt]


Degree 20 interpolant:

sqrt = gaussLegendreInterpolation[20, Sqrt, {0, Pi}][x] // Simplify // N
(*
0.0547549 + 4.91465 x - 43.8281 x^2 + 333.68 x^3 - 1747.49 x^4 +
6419.58 x^5 - 17104.3 x^6 + 34004. x^7 - 51553.3 x^8 + 60568.2 x^9 -
55757.4 x^10 + 40489.7 x^11 - 23253.5 x^12 + 10544. x^13 -
3750.35 x^14 + 1033.25 x^15 - 215.885 x^16 + 33.0421 x^17 -
3.49217 x^18 + 0.227653 x^19 - 0.00689549 x^20
*)

Plot[{Sqrt[x], sqrt}, {x, 0, Pi}, WorkingPrecision -> Precision@sqrt]


Note that as you go higher in degree, if the polynomial is put into power-basis form as above, the polynomial suffers from catastrophic numerical errors, unless high precision is used. Omit the // N to use high precision. Using barycentric interpolation is more stable than converting to the power basis, but I chose the above form because the OP seemed interested in getting the formula. The barycentric formula is a bit of a mess and looks like a rational function. Indeed you get divide-by-zero errors unless you use baryInterp[].

• It occurs to me one can use Interpolation on the Gauss-Legendre points, instead of using BarycentricInterpolation[], if the power-basis form is desired. (I started with the BC interp. in mind for numerical reasons, but as I thought through my answer, I clung to BC even after it need evaporated.) Sep 6, 2018 at 12:39

Padé (rational) approximations are much better.

asqrt[x] = PadeApproximant[Sqrt[x], {x, 1, 2}]
(* (1 + 5/4 (-1 + x) + 5/16 (-1 + x)^2)/(1 + 3/4 (-1 + x) + 1/16 (-1 + x)^2) *)
Plot[{Sqrt[x], asqrt[x]}, {x, 0, Pi}]


• Even more snug rational fit: Module[{f}, f[a_?NumericQ, b_, c_, d_] := Quiet@NIntegrate[(Sqrt[x] - x (a - x)/(b x^2 + c x + d))^2, {x, 0, Pi}]; Plot[{Sqrt[x], x (a - x)/(b x^2 + c x + d) /. Last@NMinimize[f[a, b, c, d], {a, b, c, d}, Reals]}, {x, 0, Pi}, Evaluated -> True]] - here I do assume that value of the function is zero at zero, of course. Sep 4, 2018 at 6:59
• @kirma This is great ! How do I get the resulting function ? Sep 4, 2018 at 8:29
• Thanks a lot for your answer ! Sep 4, 2018 at 8:29
• @james Module[{f}, f[a_?NumericQ, b_, c_, d_] := Quiet@NIntegrate[(Sqrt[x] - x (a - x)/(b x^2 + c x + d))^2, {x, 0, Pi}]; Plot[{Sqrt[x], x (a - x)/(b x^2 + c x + d) /. Last@NMinimize[f[a, b, c, d], {a, b, c, d}, Reals] evaluates to ((-0.6450108179399511 - x)*x)/(-0.1008496912186526 - 1.2894051380615201*x - 0.2625364225597352*x^2) Sep 4, 2018 at 8:34
• And, if you want to keep going down this rabbit hole, there's the FunctionApproximations package. Sep 4, 2018 at 11:15

You can use a Quantile Regression fit of Chebyshev polynomials through the package "MonadicQuantileRegression.m". See this blog post for details.

sqPoints = Table[{x, Sqrt[x]}, {x, 1/2, 8, 0.1}];
ListPlot[sqPoints]

p =
QRMonUnit[sqPoints]⟹
QRMonQuantileRegressionFit[20, 0.5]⟹
QRMonPlot[]⟹
QRMonErrorPlots⟹
QRMonErrors⟹QRMonEchoFunctionValue["Relative errors summary:", RecordsSummary /@ # &];


qFunc = (p⟹QRMonTakeRegressionFunctions)[0.5];

Simplify[qFunc[x]]

(* 27.05 - 245.125 x + 1011.11 x^2 - 2479.16 x^3 + 4080.66 x^4 - 4812.07 x^5 +
4233.13 x^6 - 2853.77 x^7 + 1501.85 x^8 - 624.791 x^9 + 207.112 x^10 -
54.9287 x^11 + 11.6595 x^12 - 1.97386 x^13 + 0.264376 x^14 -
0.0276295 x^15 + 0.00220351 x^16 - 0.000129425 x^17 + 5.27363*10^-6 x^18 -
1.3307*10^-7 x^19 + 1.56554*10^-9 x^20 *)
`