Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I'd like to plot the function $f$ in this answer http://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.

share|improve this question
1  
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 at 9:46
    
@Öskå See my edit please. –  G.T.R May 10 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. –  G.T.R May 10 at 9:54
1  
Do you mean g[x_]=h[Floor[x]^2 (x-Floor[x])] ? –  b.gatessucks May 10 at 10:58

2 Answers 2

up vote 8 down vote accepted

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}]

enter image description here

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]

enter image description here

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

share|improve this answer
    
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. –  G.T.R May 10 at 11:09
2  
@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 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] == 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}]

Mathematica graphics

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.