Lower and upper Riemann sums

Question

I am trying to create a few images for educational purposes. Therefore I need to be able to illustrate the area under various functions, using lower and upper riemann sums (rectangles)

Reason for not including an MWE is that I do not know which package or tools are best suited for the task, wheter it is tikz, pgf-plots, or something else. Although prefferably not gnuplot, as it is not working for me.

So my goal is to type in the function, specify the start and end value, upper or lower sum, number of rectangles, and produce something like below.

From another question, I was able to discover how to make middle riemann sums. But I absolutely must have lower and upper =(

\documentclass[10pt,a4paper]{minimal}
\usepackage{pgfplots}

\pgfplotsset{
    integral segments/.code={\pgfmathsetmacro\integralsegments{#1}},
    integral segments=3,
    integral/.style args={#1:#2}{
        ybar interval,
        domain=#1+((#2-#1)/\integralsegments)/2:#2+((#2-#1)/\integralsegments)/2,
        samples=\integralsegments+1,
        x filter/.code=\pgfmathparse{\pgfmathresult-((#2-#1)/\integralsegments)/2}
    }
}

\begin{document}

\pagecolor{black}

\begin{center}
\begin{tikzpicture}[color=white,/pgf/declare function={f=3*e^(-x)*x^3;}]
\begin{axis}[
    domain=0:8.1,
    samples=100,
    axis lines=middle
]
\addplot [ultra thick] {f};
\addplot [
    white,
    integral segments=4,
    integral=0:8
] {f};
\end{axis}
\end{tikzpicture}
\end{center}

\end{document}

Answer 1

A solution with PSTricks. Run it with xelatex or `latex->dvips->ps2pdf

\documentclass{article}
\usepackage{pstricks-add}
\begin{document}

\psset{plotpoints=200,algebraic}
\begin{pspicture}(-0.5,-2.5)(10,3)
\psStep[linecolor=magenta,StepType=upper,
  fillstyle=solid,fillcolor=cyan!50](0,9){20}{sqrt(x)*sin(x)}
\psStep[linecolor=blue,fillstyle=solid,fillcolor=blue!30,
  opacity=0.4](0,9){20}{sqrt(x)*sin(x)}
\psaxes[labelFontSize=\scriptstyle]{->}(0,0)(0,-2.25)(10,3)
\psplot[linewidth=1.5pt]{0}{10}{sqrt(x)*sin(x)}
\end{pspicture}

\psset{yunit=1.25cm}
\begin{pspicture}(-0.5,-1.75)(10,1.5)
\psaxes[labelFontSize=\scriptstyle]{->}(0,0)(0,-1.5)(10,1.5)
\psStep[StepType=Riemann,fillstyle=solid,
  fillcolor=magenta!20](0,10){50}{sqrt(x)*cos(x)*sin(x)}
\psplot[linewidth=1.5pt]{0}{10}{sqrt(x)*cos(x)*sin(x)}
\end{pspicture}
\end{document}

Answer 2

I don't claim that this is a particularly elegant solution, but this is how I did a small illustration of Riemann sums and upper and lower Darboux sums:

\documentclass[10pt]{article}
\usepackage{mathtools,relsize,tikz}
\usetikzlibrary{calc,intersections,through} % I'm not 100% on which of these are needed
\begin{figure}[h]%\centering
  % f(x) = x^3/20 - x^2/5 - x/5 + 3 from -3 to 3 (drawing -3.1 to 2.8).
  % f(-3) = 0.45,  f(3) = 1.95
  % f has max at approx. -0.4305 and f(-0.4305) ~= 3.0450.
  % We're computing the integral from -2.8 to 2.5
  % (x0=a,x1,x2,x3,x4=b) = (-2.8,     -1.3,      0.2,      1.3,      2.5)
  % (x0*,x1*,x2*,x3*)    = (     -2.1,     -0.8,      0.7,      1.9)
  \def\GraphAndAxis{%
    \draw[->] (-3.0,0) -- (2.7,0) coordinate (x axis);
    \draw[smooth,semithick,fill=none,name path=graph,domain=-2.9:2.6]% 
    plot[id=poly] (\x,0.05*\x*\x*\x - 0.2*\x*\x - 0.2*\x + 3);} % end GraphAndAxis
  \begin{minipage}{0.5\linewidth}
    \centering
    \begin{tikzpicture}
      \def\TagBox#1#2#3#4{% 1=x_i, 2=x*_i, 3=x_(i+1), 4=i
        % \useasboundingbox (0,0) rectangle (0,0);
        \path[name path=upward at tag] (#2,0.0) -- (#2,3.1);
        % In general the meaning of (p) |- (q) is "the intersection of a vertical line through p
        % and a horizontal line through q."
        \path[name intersections={of=graph and upward at tag, by=fofxstar}];
        \draw[fill=gray!40] (#3,0.0 |- fofxstar) rectangle (#1,0.0) 
        node[below,opacity=1.0] {\smaller\ensuremath{x^{\phantom{*}}_{#4}}};
        \draw[thin,dashed] (#2,0.0) node[below] {\smaller\ensuremath{x^*_{#4}}}-- (fofxstar);
      } % end TagBox
      \GraphAndAxis %draw it so \TagBox can find intersection with it
      \TagBox{-2.8}{-2.1}{-1.3}{0}
      \TagBox{-1.3}{-0.8}{0.2}{1}
      \TagBox{0.2}{0.7}{1.3}{2}
      \TagBox{1.3}{1.9}{2.5}{3}
      \GraphAndAxis
      \begin{scope}
        \path (2.5,1.0) -- (2.5,0.0) node [below] {\smaller\ensuremath{x^{\phantom{*}}_{4}}}; %ugly manual hack
      \end{scope}
    \end{tikzpicture}
    \caption{Riemann sum.}\label{Fig::riemann sum}
  \end{minipage}
  \begin{minipage}{0.5\linewidth}
    \centering
    \begin{tikzpicture}
      \def\UpperAndLowerBoxes#1#2#3#4#5{% #1=x_i, #2=x_(i+1),
        % #3=x s.t. f(x) = min_x f(x), 
        % #4=x s.t. f(x) = max_x f(x),
        % #5= index 
        \path[name path=upward at max x] (#4,0.0) -- (#4,3.1);
        \path[name intersections={of=graph and upward at max x, by=maxf}];
        \draw[fill=gray!30] (#1,0.0) rectangle (#2,0.0 |- maxf);
        \path[name path=upward at min x] (#3,0.0) -- (#3,3.1);
        \path[name intersections={of=graph and upward at min x, by=minf}];
        \draw[fill=gray!50] (#2,0.0 |- minf) rectangle (#1,0.0)
        node[below,opacity=1.0] {\smaller{\ensuremath{x^{\phantom{*}}_{#5}}}};
      } % end UpperAndLowerBoxes
      \GraphAndAxis
      \UpperAndLowerBoxes{-2.8}{-1.3}{-2.8}{-1.3}{0}
      \UpperAndLowerBoxes{-1.3}{0.2}{-1.3}{-0.4305}{1}
      \UpperAndLowerBoxes{0.2}{1.3}{1.3}{0.2}{2}
      \UpperAndLowerBoxes{1.3}{2.5}{2.5}{1.3}{3}
      \GraphAndAxis
      \begin{scope}
        \path (2.5,1.0) -- (2.5,0.0) node [below] {\smaller\ensuremath{x^{\phantom{*}}_{4}}}; %ugly manual hack
      \end{scope}
    \end{tikzpicture}
    \caption{Upper and lower Darboux sums.}\label{Fig::upper and lower sums}
  \end{minipage}
\end{figure}

Answer 3

Now it works as expected. To demonstrate the behavior I plotted two different functions.

How it works:

You have the following keys:

• integral segments

  number of the equal parts which the domain of the function is divided into;

• integral samples

  number of computation cycles for getting the maximum value in the given interval defined by integral segments

• integral min

  draw the lower Riemann sum

• integral max

  draw the upper Riemann sum

I tested the example with tikz 2.10 CSV. There is a known bug and I added a fix to the code. For more details have a look at this question. Why am I getting the PGF Math Error: Unknown function `getargs'?

Another relevant answer was provided by @Christian Feuersänger: How do I use pgfmathdeclarefunction to create define a new pgf function?

\documentclass[10pt,a4paper]{article}
\usepackage{pgfplots}
\pgfplotsset{compat=newest}
\pgfplotsset{
    integral segments/.code={\pgfmathsetmacro\integralsegments{#1}},
    integral segments=3,
    integral samples/.code={\edef\integralsamples{#1}},
    integral samples = 10,
    integral min/.style args={#1:#2}{
        ybar interval,
        domain=#1:#2,
        samples=\integralsegments+1,
        x filter/.code={
                               \edef\lastx{\pgfmathresult}
                                \pgfmathresult
                              },%
        y filter/.code={%
                          \pgfmathparse{(#2/(\integralsegments))/\integralsamples}%
                           \edef\tempstep{\pgfmathresult}%
                            \pgfmathparse{f(\lastx)}%
                            \edef\tempa{\pgfmathresult}%
                            \edef\tempb{\pgfmathresult}%
                            \foreach \x in {0,1,...,\integralsamples}%
                               {%
                                \pgfmathparse{f(\lastx+\x*\tempstep)}%
                                 \xdef\tempb{\tempb,\pgfmathresult}%
                                 }%
                               \pgfmathmin{\tempb}{\tempb}
                             },
    },
    integral max/.style args={#1:#2}{
        ybar interval,
        domain=#1:#2,
        samples=\integralsegments+1,
        x filter/.code={
                               \edef\lastx{\pgfmathresult}
                                \pgfmathresult
                              },%
        y filter/.code={%
                          \pgfmathparse{(#2/(\integralsegments))/\integralsamples}%
                           \edef\tempstep{\pgfmathresult}%
                            \pgfmathparse{f(\lastx)}%
                            \edef\tempa{\pgfmathresult}%
                            \edef\tempb{\pgfmathresult}%
                            \foreach \x in {0,1,...,\integralsamples}%
                               {%
                                \pgfmathparse{f(\lastx+\x*\tempstep)}%
                                 \xdef\tempb{\tempb,\pgfmathresult}%
                                 }%
                               \pgfmathmax{\tempb}{\tempb}
                             },
    },   
}
\makeatletter
%see http://tex.stackexchange.com/questions/9722/why-am-i-getting-the-pgf-math-error-unknown-function-getargs
\def\pgfmathmax#1#2{%
%   \pgfmathparse{getargs(#1,#2)}%
     \pgfmathparse{#1,#2}%
    \expandafter\pgfmathmax@\expandafter{\pgfmathresult}%
}
\def\pgfmathmin#1#2{%
    \pgfmathparse{#1,#2}%
    \expandafter\pgfmathmin@\expandafter{\pgfmathresult}%
}
%see http://tex.stackexchange.com/questions/15435/how-do-i-use-pgfmathdeclarefunction-to-create-define-a-new-pgf-function
\makeatother
\begin{document}
\begin{tikzpicture}
\pgfset{declare function={f(\x)=3*exp(-(\x))*(\x)^3+1;}}
\begin{axis}[
    domain=0:8.1,
    samples=100,
    axis lines=middle
]
\addplot [ultra thick] {f(x)};

\addplot [
    red,
    integral segments=4,
    integral min=0:8
] {f(x)};
\end{axis}
\end{tikzpicture}
\hfill
\begin{tikzpicture} 
\pgfset{declare function={f(\x)=3*exp(-(\x))*(\x)^3+1;}}
\begin{axis}[
    domain=0:8.1,
    samples=100,
    axis lines=middle
]
\addplot [ultra thick] {f(x)};

\addplot [
    blue,
    integral segments=5,
    integral max=0:8
] {f(x)};
\end{axis}
\end{tikzpicture}

\begin{tikzpicture}
\pgfset{declare function={f(\x)=sin(2*deg(\x))*exp(0.1*\x)+2;}}
\begin{axis}[
    domain=0:8.1,
    samples=100,
    axis lines=middle
]
\addplot [ultra thick] {f(x)};

\addplot [
    red,
    integral segments=4,
    integral min=0:8
] {f(x)};
\end{axis}
\end{tikzpicture}
\hfill
\begin{tikzpicture} 
\pgfset{declare function={f(\x)=sin(2*deg(\x))*exp(0.1*\x)+2;}}
\begin{axis}[
    domain=0:8.1,
    samples=100,
    axis lines=middle
]
\addplot [ultra thick] {f(x)};

\addplot [
    blue,
    integral segments=5,
    integral max=0:8
] {f(x)};
\end{axis}
\end{tikzpicture}

\end{document}

TEST CODE to handle functions with f(x)<0:

\documentclass[10pt,a4paper]{article}
\usepackage{pgfplots}
\pgfplotsset{compat=newest}
\pgfplotsset{
    integral segments/.code={\pgfmathsetmacro\integralsegments{#1}},
    integral segments=3,
    integral samples/.code={\edef\integralsamples{#1}},
    integral samples = 10,
    integral min/.style args={#1:#2}{
        ybar interval,
        domain=#1:#2,
        samples=\integralsegments+1,
        x filter/.code={
                               \edef\lastx{\pgfmathresult}
                                \pgfmathresult
                              },%
        y filter/.code={%
                          \pgfmathparse{(#2/(\integralsegments))/\integralsamples}%
                           \edef\tempstep{\pgfmathresult}%
                            \pgfmathparse{f(\lastx)}%
                            \edef\tempa{\pgfmathresult}%
                            \edef\tempb{\pgfmathresult}%
                            \foreach \x in {0,1,...,\integralsamples}%
                               {%
                                \pgfmathparse{f(\lastx+\x*\tempstep)}%
                                 \xdef\tempb{\tempb,\pgfmathresult}%
                                 }%
                               \pgfmathmin{\tempb}{\tempb}
                               \let\savepgfmathresult\pgfmathresult
                               \pgfmathgreater{\pgfmathresult}{0}
                               \ifdim\pgfmathresult pt> 0 pt \relax
                                         \let\pgfmathresult\savepgfmathresult 
                               \else
                                    \pgfmathparse{f(\lastx)}%
                                    \edef\tempa{\pgfmathresult}%
                                    \edef\tempb{\pgfmathresult}%
                                      \foreach \x in {0,1,...,\integralsamples}%
                                         {%
                                          \pgfmathparse{f(\lastx+\x*\tempstep)}%
                                           \xdef\tempb{\tempb,\pgfmathresult}%
                                           }%
                                         \pgfmathmax{\tempb}{\tempb}
                               \fi 
                             },
    },
    integral max/.style args={#1:#2}{
        ybar interval,
        domain=#1:#2,
        samples=\integralsegments+1,
        x filter/.code={
                               \edef\lastx{\pgfmathresult}
                                \pgfmathresult
                              },%
        y filter/.code={%
                          \pgfmathparse{(#2/(\integralsegments))/\integralsamples}%
                           \edef\tempstep{\pgfmathresult}%
                            \pgfmathparse{f(\lastx)}%
                            \edef\tempa{\pgfmathresult}%
                            \edef\tempb{\pgfmathresult}%
                            \foreach \x in {0,1,...,\integralsamples}%
                               {%
                                \pgfmathparse{f(\lastx+\x*\tempstep)}%
                                 \xdef\tempb{\tempb,\pgfmathresult}%
                                 }%
                               \pgfmathmax{\tempb}{\tempb,0}
                               \let\savepgfmathresult\pgfmathresult
%                               \pgfmathparse{ifthenelse(\pgfmathresult>=0,1,0)} 
                               \pgfmathgreater{\pgfmathresult}{0}
                               \ifdim\pgfmathresult pt> 0 pt \relax
                                         \let\pgfmathresult\savepgfmathresult 
                               \else
                                    \pgfmathparse{f(\lastx)}%
                                    \edef\tempa{\pgfmathresult}%
                                    \edef\tempb{\pgfmathresult}%
                                      \foreach \x in {0,1,...,\integralsamples}%
                                         {%
                                          \pgfmathparse{f(\lastx+\x*\tempstep)}%
                                           \xdef\tempb{\tempb,\pgfmathresult}%
                                           }%
                                         \pgfmathmin{\tempb}{\tempb}
                               \fi 
                             },
    },   
}
\makeatletter
%see http://tex.stackexchange.com/questions/9722/why-am-i-getting-the-pgf-math-error-unknown-function-getargs
\def\pgfmathmax#1#2{%
%   \pgfmathparse{getargs(#1,#2)}%
     \pgfmathparse{#1,#2}%
    \expandafter\pgfmathmax@\expandafter{\pgfmathresult}%
}
\def\pgfmathmin#1#2{%
    \pgfmathparse{#1,#2}%
    \expandafter\pgfmathmin@\expandafter{\pgfmathresult}%
}
%see http://tex.stackexchange.com/questions/15435/how-do-i-use-pgfmathdeclarefunction-to-create-define-a-new-pgf-function
\makeatother
\begin{document}
\begin{tikzpicture}
\pgfset{declare function={f(\x)=3*exp(-(\x))*(\x)^3+1;}}
\begin{axis}[
    domain=0:8.1,
    samples=100,
    axis lines=middle
]
\addplot [ultra thick] {f(x)};

\addplot [
    red,
    integral segments=4,
    integral min=0:8
] {f(x)};
\end{axis}
\end{tikzpicture}
\hfill
\begin{tikzpicture} 
\pgfset{declare function={f(\x)=3*exp(-(\x))*(\x)^3+1;}}
\begin{axis}[
    domain=0:8.1,
    samples=100,
    axis lines=middle
]
\addplot [ultra thick] {f(x)};

\addplot [
    blue,
    integral segments=5,
    integral max=0:8
] {f(x)};
\end{axis}
\end{tikzpicture}

\begin{tikzpicture}
\pgfset{declare function={f(\x)=sin(2*deg(\x))*exp(0.1*\x)+2;}}
\begin{axis}[
    domain=0:8.1,
    samples=100,
    axis lines=middle
]
\addplot [ultra thick] {f(x)};

\addplot [
    red,
    integral segments=4,
    integral min=0:8
] {f(x)};
\end{axis}
\end{tikzpicture}
\hfill
\begin{tikzpicture} 
\pgfset{declare function={f(\x)=sin(2*deg(\x))*exp(0.1*\x)+2;}}
\begin{axis}[
    domain=0:8.1,
    samples=100,
    axis lines=middle
]
\addplot [ultra thick] {f(x)};

\addplot [
    blue,
    integral segments=5,
    integral max=0:8,
] {f(x)};
\end{axis}
\end{tikzpicture}


\begin{tikzpicture}
\pgfset{declare function={f(\x)=sqrt(\x)*cos(deg(\x))*sin(deg(\x));}}
\begin{axis}[
    domain=0:8.1,
    samples=100,
    axis lines=middle
]
\addplot [ultra thick] {f(x)};

\addplot [
    red,
    integral segments=30,
    integral min=0:8
] {f(x)};
\end{axis}
\end{tikzpicture}
\hfill
\begin{tikzpicture} 
\pgfset{declare function={f(\x)=sqrt(\x)*cos(deg(\x))*sin(deg(\x));}}
\begin{axis}[
    domain=0:8.1,
    samples=100,
    axis lines=middle
]
\addplot [,ultra thick] {f(x)};

\addplot [
    blue,
    integral segments=15,
    integral max=0:8
] {1};
\end{axis}
\end{tikzpicture}

\end{document}

Answer 4

This seems to work but sort of hackish.

\documentclass{article}
\usepackage{pgfplots}
\begin{document}
\begin{tikzpicture}
\begin{axis}
\addplot[const plot mark right,fill=yellow!20]
coordinates
{(0,0.1) (0.1,0.15) (0.2,0.5) (0.3,0.62)
(0.4,0.56) (0.5,0.58) (0.6,0.65) (0.7,0.6)
(0.8,0.58) (0.9,0.55) (1,0.52) }\closedcycle;

\addplot[const plot,fill=blue!20,opacity=0.5]
coordinates
{(0,0.1) (0.1,0.15) (0.2,0.5) (0.3,0.62)
(0.4,0.56) (0.5,0.58) (0.6,0.65) (0.7,0.6)
(0.8,0.58) (0.9,0.55) (1,0.52) }\closedcycle;

\addplot+[smooth]
coordinates
{(0,0.1) (0.1,0.15) (0.2,0.5) (0.3,0.62)
(0.4,0.56) (0.5,0.58) (0.6,0.65) (0.7,0.6)
(0.8,0.58) (0.9,0.55) (1,0.52)};
\end{axis}
\end{tikzpicture}
\end{document}

See the manual about const plot mark <right/left/mid>

Answer 5

As is mentioned above, TeX is not a good CAS, so I would suggest the following (you'll need 3 tabs open in your browser): First, staying on this page, create another tab, go here and find "Numerical integrals with various rules". Highlight the code from that section and copy it (Control-C). Second, in the third tab go here. Paste the code in over the existing code (Control-V). Press 'Evaluate' on the screen to make sure the code is working properly (i.e. you've copied everything properly [spacing/indentation is important). You should now have a manipulative that does what you've asked for, except there's no PDF output. We'll change that. On the Sage Cell Server page, go to the line of code that says, " show(plot(func,a,b) + rects, xmin = a, xmax = b, ymin = min_y, ymax = max_y)". That 1 line will be replaced by 3 lines.

C= plot(func,a,b) + rects 
C.show(xmin = a, xmax = b, ymin = min_y, ymax = max_y)
C.save("MyPicture.pdf")

Again, spacing/indentation is critical. Press 'Evaluate' and you'll get output like what is shown below. I couldn't get everything on the screen, but under the display (below 'Session Notes' there is a link to "MyPicture.pdf.". Right click on it and save it to your computer. Note, it only shows the graph, not the HTML Riemann sum calculation or sliders that you see in the screenshot.

It shouldn't take more than several minutes to run through all the steps. Here is a screenshot which includes the altered code section and image.