How to plot this function?

Question

I'd like to plot the function $f$ in this answer https://math.stackexchange.com/a/788818/66096

Let $h(x) = \begin{cases} e^{ -\frac{1}{1 - (1-2x)^2}} & \mbox{ for } 0<x < 1\\ 0 & \mbox{ otherwise} \end{cases}$ Then let $$g(x)=\sum_{n=1}^\infty h(n^2(x-n))=\begin{cases}h(n^2(x-n))&\text{if }n\le x<n+1, n\in\mathbb N\\0&\text{otherwise}\end{cases}$$ And let $f(x)=\int_0^xg(t)\,\mathrm dt$.

I only know the basics of Mathematica, so I've been trying to plot piecewise functions without success.

Here is what I got

h[x_] = Piecewise[{{E^(-(1 - (1 - 2 x)^2)^(-1)), 0 < x < 1}, {0, 
    x > 1}, {0, x < 0}}]

But how do I deal with $g$ ? I could input partial sums with the series expression, but at the same time, I could save computations by just defining it piecewise.

Answer 1

Pretty cool example! I'm sure the following could be sped up a bit.

h[x_] = Piecewise[{{Exp[-1/(1 - (1 - 2 x)^2)], 0 < x < 1}}];
g[x_] := With[{n = Floor[x]}, h[n^2 (x - n)]];
f[x_?NumericQ] := NIntegrate[g[t], {t, 0, x}];
Plot[f[x], {x, 0, 8}]

A plot of $g(x)$ helps illustrate what's going on here as well.

Plot[g[x], {x, 0, 9.1},
  PlotPoints -> 2000,
  PlotRange -> All]

The humps are basically terms in the definition of $g$. The fact that they are non-overlapping allows us to define $g$ without resorting to a sum. $f$ is just an accumulated area function for $g$.

Answer 2

Here are two ways to speed up the computation. Both use NDSolve

(1) The simplest, for a given interval:

Clear[f, g, h];
h[x_?NumericQ] := Piecewise[{{Exp[-1/(1 - (1 - 2 x)^2)], 0 < x < 1}}];
g[x_] := With[{n = Floor[x]}, h[n^2 (x - n)]];
f = NDSolveValue[{f0'[t] == g[t], f0[0] == 0}, f0, {t, 0, 8},
   Method -> "StiffnessSwitching"];

(2) Alternative, fast and valid for all (nonnegative) real numbers: It takes advantage of the facts that the infinite sum is really finite because $h$ is zero for arguments greater than $1$ and that the integral of $g$ is just a sum of the integral of $h$, under a substitution $u = n^2(x-n)$.

Clear[f];
if = NDSolveValue[{f'[x] == Exp[-1/(1 - (1 - 2 x)^2)], 
       f[1/2] == NIntegrate[Exp[-1/(1 - (1 - 2 K[1])^2)], {K[1], 0, 1/2}, 
         AccuracyGoal -> 15, PrecisionGoal -> 15, WorkingPrecision -> 30]},
       f, {x, 0 + 10^-5, 1 - 10^-5}, 
       AccuracyGoal -> 15, PrecisionGoal -> 15, WorkingPrecision -> 30];
if1 = NIntegrate[Exp[-1/(1 - (1 - 2 K[1])^2)], {K[1], 0, 1}, 
       AccuracyGoal -> 15, PrecisionGoal -> 15, WorkingPrecision -> 30];
f[x_] := Piecewise[{
          {Sum[1/n^2, {n, Floor[x] - 1}] if1 + 
            1/Floor[x]^2 *
             Piecewise[{{Quiet@if[Floor[x]^2 (x - Floor[x])], x - Floor[x] < Floor[x]^2}}],
           x > 1}}];

Both produce:

Plot[f[x], {x, 0, 8}]