# How to plot this function?

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.

• Let us see what you tried first :) So we don't have to start typing every single letter but only fix your code :) – Öskå May 10 '14 at 9:46
• @Öskå See my edit please. – Gabriel Romon May 10 '14 at 9:50
• @Öskå Yes I checked it. I really don't know how to proceed in the same manner with $g$ since it has infinitely many pieces. – Gabriel Romon May 10 '14 at 9:54
• Do you mean g[x_]=h[Floor[x]^2 (x-Floor[x])] ? – b.gates.you.know.what May 10 '14 at 10:58

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$.

• The output looks exactly like what I had in mind. Could you elaborate on the "With[{n=Floor[x]}]" and NumericQ ? It's the first time I see them. Thanks. – Gabriel Romon May 10 '14 at 11:09
• @GabrielR. _?NumericQ is useful to prevent evaluation prior to an actual number being plugged in and I use it out of habit when defining functions using numerical procedures. Strictly speaking, it's not necessary here. If you were to try to plot $f'$, however, some errors would be issued. The With construct is simply a scoping mechanism to declare n as local. Module could have been used just as well. – Mark McClure May 10 '14 at 11:18

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}, 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)^2)], {K, 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)^2)], {K, 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}] 