Approximation of the Fabius function with a quotient of exponentials

Question

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.

---

Added after the question was answered

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<b<c<d$ 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)$$

Answer 1

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.

Answer 2

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.

Example:

---

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)

---

See also: How do I numerically evaluate and plot the Fabius function?

Answer 3

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!

addendum

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!

Answer 4

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.