# Approximation of the Fabius function with a quotient of exponentials

Approximation of the Fabius function $$f(x) = \text{FabiusF}[x+1]\cdot \text{HeavisideTheta}[1-x^2]$$ - FabiusF[x] doesn't work in Wolfram-Alpha

I am looking to figure out how well the displaced version of the Fabius function $$f(x) = \text{FabiusF}[x+1],\quad |x|\leq 1$$ could being approximated by the function: $$g(x) = \frac{1+\exp\left(\frac{-1}{|x|^{\frac{3}{4}\sqrt{2\pi}}}\right)}{1+\exp\left(\frac{1-2\ |x|}{x^2-|x|}\right)}$$ but I don't have Mathematica and I use instead Wolfram-Alpha where the function FabiusF[x] is not working as you could see here, so I would like to ask for help in order to compare in a plot $$f(x)$$ and $$g(x)$$ for $$x\in [-1,\ 1]$$ (just a picture will work), and also if possible the plot of $$f'(x)$$ and $$g'(x)$$.

This is related to this other question, so if you want to improve it by changing some constants, you are also welcome. "Supposedly" $$f(x)$$ should fulfill $$f'(x) = 2f(2x+1)-2f(2x-1)$$ which is what I am trying to match in the approximation (or something "similar" like $$f'(x) = k\left(f(2x+1)-f(2x-1)\right)$$ for some $$k \in \mathbb{R}$$).

Thanks a lot beforehand.

I have realized later that the following function: $$q(x)=\frac{1}{1+\exp\left(\frac{1-2|x|}{x^2-|x|}\right)}$$ Matches "almost perfectly" the comstruction of a "smooth" transition curve show in Wikipedia Non-analytic smooth function:

1. Define $$f(x)=\begin{cases}e^{-\frac{1}{x}},\quad x>0\\ 0,\quad x\leq 0\end{cases}$$
2. Define $$p(x)=\frac{f(x)}{f(x)+f(1-x)}$$
3. Then for some real valued constants $$a one could built a smooth bump function as: $$r(x)=p\left(\frac{x-a}{b-a}\right)\cdot p\left(\frac{d-x}{d-c}\right)$$ which is non-zero in $$[a,\ d]$$ such is flat on $$[b,\ c]$$

One can "experimentally" see that for $$[a,\ d]\equiv [-1,\ 1]$$, by chosing $$b\to 0,\ c\to 0$$ then $$q(x)$$ matches $$r(x)$$ on $$(-1,\ 1)\setminus\{0\}$$.

Then, by playing with $$q(x)$$ I realized two interesting things about the function: $$q(x,a)=\frac{1}{1+\exp\left(\frac{a(1-2|x|)}{x^2-|x|}\right)}$$

• The function $$q(x,4)$$ fits "really good" the function $$r(x)$$ if I change the definition of $$f(x)$$ to $$f(x)=\begin{cases}e^{-\frac{1}{x^2}},\quad x>0\\ 0,\quad x\leq 0\end{cases}$$
• The function $$q(x,\frac{\sqrt{3}}{2})$$ "looks like" a perfect smooth transition function with straight-line edges, fulfilling having a derivative with flat tops - but unfortunately it don't fulfill solving a DDE like $$y'(x)=k\left(y(2x+1)-y(2x-1)\right)$$ for some constant $$k$$.

You could see them on Desmos:

The lines that fit the slope of $$q(x,\frac{\sqrt{3}}{2})$$ are: $$L^\pm=\pm \sqrt{3}\ (x\pm 1)-\frac12(\sqrt{3}-1)$$

• Comments have been moved to chat; please do not continue the discussion here. Before posting a comment below this one, please review the purposes of comments. Comments that do not request clarification or suggest improvements usually belong as an answer, on Mathematica Meta, or in Mathematica Chat. Comments continuing discussion may be removed.
– Kuba
Commented Jun 2, 2023 at 10:49
• @VladimirReshetnikov In this question I explained in detail how I got the approximated function. Commented Jun 30, 2023 at 12:17

This question is very interesting since the Fabius function is a little mathematical treasure and many of its properties can be successfully demonstrated with Mathematica. It has a potential of opening new and surprising perspectives on "well-known" theories of mathematics and physics. FabiusF can be still improved, nonetheless it already works well.

Let's define

fb = ResourceFunction["FabiusF"];


We can evaluate the Fabius function exactly at half integer arguments:

Table[fb[x/2], {x, 0, 8}]

  {0, 1/2, 1, 1/2, 0, -(1/2), -1, -(1/2), 0}


We could expect that more symbolic properties can be exploited in next versions of this function, especially when it is built-in. We introduce

g[x_] /; -1 < x < 1 := (1 + Exp[-(1/RealAbs[x]^(3/4 Sqrt[2 Pi]))])/(
1 + Exp[(1 - 2 RealAbs[x])/(x^2 - RealAbs[x])])
g[x_] /; x <= -1 || x >= 1 := 0


The above definition of g function makes it exactly zero outside the range $$(-1,1)$$. We have used RealAbs instead of Abs to evaluate properly derivatives. Now plotting the both functions we find that g is a a quite good approximation of fb:

Plot[{fb[x + 1], g[x]}, {x, -1.02, 1.02}, AspectRatio -> Automatic,
PlotStyle -> {{Thickness[0.003]}, {Thick, Dashing[0.01]}}]


Here we find maximal deviations between fb and g searching for numerical extremal values, symbolic counterparts are not able to evaluate them:

Through @ { NMaximize, NMinimize}[{fb[x + 1] - g[x], 0 <= x <= 1}, x]

  {{0.00336294, {x -> 0.837602}}, {-0.0254221, {x -> 0.632589}}}

Plot[fb[x + 1] - g[x], {x, -1, 1}, PlotStyle -> Thick, AspectRatio -> 1/4] // Quiet


Here one can wonder whether there might be better approximations and Weierstrass Approximation Theorem says there are better ones even among polynomials since we are interested in approximations on a compact set. See e.g. Generating a polynomial that's accurate to within an error of no more than 1/10^5 and Get polynomial interpolation formula.

Let's compare functions suggested in the commets:

u[x_] /; -1 <= x <= 1 := fb[x + 1]
u[x_] /; x < -1 || x > 1 := 0
Plot[{u'[x], 2 (u[2 x + 1] - u[2 x - 1])}, {x, -1, 1},
PlotStyle -> {Dashing[0.008], Dashing[0.004]}]


These functions clearly overlap, however currently we cannot demonstrate symbolically equality of these functions, e.g. PossibleZeroQ[u'[x] - 2 (u[2 x + 1] - u[2 x - 1]), Assumptions -> -1 < x < 1] yields False. FullSimplify doesn't work either.

Next we re-evaluate first derivatives with appropriate values of fb yielding a kind of differential equation for fb:

Plot[{fb'[x], 2 (fb[2 x] - fb[2 (x - 1)])}, {x, 0, 1},
PlotStyle -> {Dashing[0.008], Dashing[0.004]}]


We have numerical deviations

Plot[fb'[x] - 2 (fb[2 x] - fb[2 (x - 1)]), {x, 0, 1},
PerformanceGoal -> "Quality", WorkingPrecision -> 20]


however we can significantly reduce numerical error setting larger WorkingPrecision.

• @Joako You are welcome, I'm glad you find it helpful. I compared numerically the both functions in my answer, however they don't overlap, changing the arguments they do. Symbolically we cannot prove it if not digging deeper. Commented May 31, 2023 at 17:36
• @Joako Your observation can be regarded as a straightforward consequence of the above calculations. I'll update my answer a bit later when I find time enough. Commented May 31, 2023 at 18:19
• @Joako Yes, we have shown numerically that equation. Commented Jun 1, 2023 at 1:53
• @Joako For any $x<0$ we haxe $F(x)=0$ so it is satisfied trivially. For $x>0$ it should be clear if you evaluate e.g. Plot[fb[x], {x, 0, 36}, PlotRange -> All]. I guess you should invest a bit of your resources and buy Mathematica (home edition or student edition or just test the system without paying for it (for a month or two weeks only). I can say it deserves its cost. It can make many things easier. Nonetheless in general one should beware of addiction to computer algebra systems. Commented Jun 1, 2023 at 23:14
• @Joako In case you haven't yet seen these papers by Juan Arias de Reyna, I recommend them very much: arxiv.org/abs/1702.05442, arxiv.org/abs/1702.06487 Commented Jul 20, 2023 at 15:38

The Fabius function is approximated well by this piecewise polynomial function:

Sum[(-1)^ThueMorse[n] (2 n + 1 - x 2^m)^(m - 1) Sign[2 n + 1 - x 2^m], {n, 0, 2^m - 1}] *
(-1)^m/(2 * 2^Binomial[m, 2] (m - 1)!)


It gives the exact Fabius function in the limit m -> ∞. It's derived from the Fourier transform of the Fabius function, which has a representation as an infinite product of Sinc functions, by truncating that product after a finite number of terms, and converting it back using the inverse Fourier transform.

You can use this expression with WolframAlpha if you substitute an integer for m.

There is also an exact formula for the values of the Fabius function at dyadic rationals using finite summation, although I don't see how to easily use it to plot the function:

FabiusF[m/2^n] == Sum[QBinomial[n, k, 1/2]/(2^(k (k - 1)) (n + k)!)
Sum[(-1)^ThueMorse[j] (j - m*2^k + 1/2)^(n + k), {j, 0, m*2^k - 1}],
{k, 0, n}]/(2^(n^2) QPochhammer[1/2, 1/2, n])]


The Fabius function has a series expansion in Legendre polynomials that is valid for 0 < x < 2:

FabiusF[x] == Sum[LegendreP[2 n, x - 1] 4^(-n) (4 n + 1)
Sum[(-1)^(n + k) Binomial[2 n, n + k] Binomial[2 n + 2 k, 2 n] (2 k)!
2^Binomial[2 k + 1, 2] FabiusF[2^(-2 k - 1)], {k, 0, n}], {n, 0, ∞}]


You can get a good polynomial approximation by truncating the outer series and replacing FabiusF[2^(-2 k - 1)] factors with corresponding rational numbers.

A very good simple polynomial approximation on 0 < x < 1/2:

FabiusF[x] ≈ 32 x^5 (7 - 28 x + 40 x^2 - 20 x^3)


You can hardly see any difference between their graphs. In a certain sense, this approximation is the best among approximations by interpolating polynomials.

In case you need an exact asymptotic of the Fabius function for x -> 0 in terms of elementary function and the Lambert 𝑊-function:

FabiusF[x] ∝ Exp[(EulerGamma^2/2 - π^2/12 - ProductLog[-1, -x Log[2]] -
1/2 ProductLog[-1, -x Log[2]]^2 + StieltjesGamma[1])/Log[2]]/(2^(7/12) √π √x)


• I fixed the typo, thanks Commented May 31, 2023 at 21:37
• It's a piecewise polynomial, it's stitched from 32 quartic polynomial pieces. You can see it as the result of the iterated refinement process described here: math.stackexchange.com/q/240687/19661 Commented May 31, 2023 at 22:02
• I also added a formula for the exact values of the Fabius function at dyadic rationals, but I'm not sure if it's useful for plotting. Theoretically, because the dyadic rationals are dense in ℝ, and the Fabius function is continuous, this formula determines its values at all points. Commented May 31, 2023 at 22:28
• @Joako Yes, because of the symmetry. Also, you can verify it using its antiderivative. Commented Jun 19, 2023 at 20:06
• @Joako There is a representation of Thue–Morse via binomial coefficients: math.stackexchange.com/q/1891648/19661. Also, you can express it via a hypergeometric function: mathworld.wolfram.com/… Commented Jun 22, 2023 at 23:43

Extended comment:

Optimal exponent a=? for the approximation

If we take the general ansatz FabiusF[x+1]~(1 + Exp[-1/x^(a Sqrt[2 Pi])])/(1 + Exp[(1 - 2 x)/(x^2 - x)]), 0<x<1, a>0 the approximation error depends on a:

err[a_?NumericQ] :=NIntegrate[(ResourceFunction["FabiusF"][x + 1] - (1 + Exp[-1/x^(a Sqrt[2 Pi])])/(1 + Exp[(1 - 2 x)/(x^2 - x)]))^2, {x, 0, 1}]

pic = LogPlot[err[a], {a, 0, 2.5}, AxesLabel -> {"a", "err[a]"} ](* takes some time...*)


In the plot we see a minimum at a~1.64 !

pi = Cases[pic, Line[p_] :> p, -1][[1]] ;
opta = Select[pi, Min[pi] == #[[2]] &][[1]] (*{1.64427, -11.8112}*)


Comparing for different a-values

Plot[Evaluate[
Map[ResourceFunction["FabiusF"][x + 1] - (1 + Exp[-1/x^(# Sqrt[2 Pi])])/(1 + Exp[(1 - 2 x)/(x^2 - x)]) &, { 3/4, opta[[1]], GoldenRatio}]], {x, 0, 1},PlotStyle -> {Blue, {Red, Dashed}, Gray},
PlotLegends ->LineLegend[{ 3/4, 1.64427, GoldenRatio // N}, LegendLabel -> "parameter a:"],PlotRange->All]


it looks like parameter value a=1.64427 ~GoldenRatio  gives a better approximation for the fabius-function!

comparison of the derivatives of the approximation error:

Plot[Evaluate[
Map[Derivative[1, 0][
Function[{x, a},
ResourceFunction["FabiusF"][x + 1] - (
1 + Exp[-1/x^(# Sqrt[2 Pi])])/(
1 + Exp[(1 - 2 x)/(x^2 - x)])]][u, a] &, { 3/4, opta[[1]],
GoldenRatio}]], {u, 0, 1},
PlotStyle -> {Blue, {Red, Dashed}, Gray},
PlotLegends ->
LineLegend[{ 3/4, 1.64427, GoldenRatio // N},
LegendLabel -> "parameter a:"], PlotRange -> All]


Hope it helps!

• The problem with your value is that we were "optimizing" different things: my objective where to have "flat" lobes in the derivative, which happen near the value $1.87$, if you change the exponent with something bigger than $2.5$ its become rounded, as you could see here in Desmos, at your value the derivative already looks rounded as the derivative $g'(x)$ of $g(x)=\left(1+\exp\left(\frac{1+2|x|}{x^2-|x|}\right)\right)$ which is different of the Fabius function $\text{FabiusF}[x+1]$ that do have flat lobes in its derivatives Commented Jun 1, 2023 at 15:39
• My idea was founding at least a function that was such as its derivatives where proportional to displaced versions of itself $u'(x)=k\left(u(2x+1)-u(2x-1)\right)$ for some $k>0$, and since the original functions have a flat top, I was trying to find also a function which derivatives have also flat top lobes Commented Jun 1, 2023 at 15:43
• @Joako My approximation gives smaller slopes too, I think! Or did I misunderstand "flat lobes in the derivative"? Commented Jun 2, 2023 at 9:07
• @Joako The quadratic error err[a_?NumericQ] := NIntegrate[(ResourceFunction["FabiusF"][x + 1] - 1 /( 1 + Exp[a (1 - 2 x)/(x^2 - x)]))^2 , {x, 0, 1}] shows (LogPlot[err[a], {a, 0, 2}, AxesLabel -> {"a", "err[a]"} ] )a minimum for a==1. Not really surprising... Commented Jun 8, 2023 at 15:58
• @Joako NMinimize[{Log[err2[a]], 0 < a < 2}, a, Method -> "NelderMead"] evaluates {-11.6135, {a -> 0.997659}} Commented Jun 8, 2023 at 16:05

As I mentioned in my comment, I didn't take limits, I just used Table in MMA to plug numbers in and only did this as a very quick very simple attempt to get you something.

Please, be very cautious with this and check all the details of this before you even think of beginning to trust this. I didn't study the links that others supplied to explore the depths of evaluating the Fabius function. I just did the three lines of MMA code to let you peek at the results and hoped that would be enough to help you. If you study the depths of the steps of precise evaluation of your approximation then you may be able to understand why MMA is concerned about the value of this at 0.

Plot will sometimes have problems just plugging in numbers.

In the next step I used /.x->0 to just plug in zero.

In the final step I used Limit to check that result at zero.

If you find any mistakes in what I did then I am sorry and if you try to simply explain how I need to fix those then I will try to do what I can. Hope this helps you get what you need done.

• thanks you very much for taking the time to answer. I don't know why your functions don't match the full plot, so I recomend you to check the answer by @Artes where the graphs do match. Commented Jun 1, 2023 at 1:12