How to plot a function in integral form with TikZ?

Question

how to plot the function $f:x\mapsto \int_x^{2x}\frac{4}{\sqrt{1+t^4}}\, \textrm{d}t$ with TikZ?

Answer 1

MWE using adaptive Simpson integration (Asymptote):

% s.tex:
\documentclass{article}
\usepackage[inline]{asymptote}
\usepackage{lmodern}
\begin{document}
\begin{asy}
size(300,200,IgnoreAspect);
import graph;
real F(real t){return 4/sqrt(1+t^4);}
real f(real x){return simpson(F,x,2x);}
pen axPen=darkblue;
pen fPen=red+1bp;
draw(graph(f,-7,7,n=200),fPen);
string noZero(real x) {return (x==0)?"":string(x);}
defaultpen(fontsize(10pt));
xaxis(axPen,LeftTicks(noZero,Step=2));
yaxis(axPen,RightTicks(noZero,Step=0.5));

label("$f:x\mapsto \displaystyle\int_x^{2x}"
     +"\frac{4}{\sqrt{1+t^4}}\, \textrm{d}t$"
  ,(1.7,f(1.7)),NE);
\end{asy}
\end{document}    

% To process it with `latexmk`, create file `latexmkrc`:
% 
%     sub asy {return system("asy '$_[0]'");}
%     add_cus_dep("asy","eps",0,"asy");
%     add_cus_dep("asy","pdf",0,"asy");
%     add_cus_dep("asy","tex",0,"asy");
% 
% and run `latexmk -pdf s.tex`.

Answer 2

Here is the PSTricks answer. I slightly changed the \psCumIntegral macro from pst-func to account for the different integration limits:

\documentclass[preview, varwidth, border=5pt]{standalone}
\usepackage{pst-func}

\makeatletter
\def\psMyIntegral{\pst@object{psMyIntegral}}
\def\psMyIntegral@i#1#2#3{%
  \begin@OpenObj%
  \addto@pscode{
    /xStart #1  def
    /dx #2 #1 sub \psk@plotpoints\space div def
    /a #1 def
    /b a 2 mul def
    /scx { \pst@number\psxunit mul } def
    /scy { \pst@number\psyunit mul } def
    tx@FuncDict begin /SFunc { #3 } def end
    \psk@plotpoints 1 add {
      a b \psk@Simpson
      tx@FuncDict begin Simpson I end
      scy a scx exch a xStart eq {moveto}{lineto}ifelse
      /a a dx add def
      /b a 2 mul def
    } repeat
  }%
  \end@OpenObj%
}
\makeatother

\begin{document}

\psset{xunit=0.8,yunit=1.5}
\begin{pspicture}(-7,-2)(7,2)
  \psMyIntegral[plotpoints=500, linecolor=red]{-7}{7}{4 exp 1 add sqrt 4 exch div}
  \psaxes[Dy=0.5, arrows=->](0,0)(-7,-2)(7,2)
  \rput[rt](7,2){$f:x\mapsto \displaystyle\int_x^{2x} \frac{4}{\sqrt{1+t^4}}\, \textrm{d}t$}
\end{pspicture}
\end{document}

That gives:

EDIT: Here comes a more general macro \psVarIntegral, which allows to specify both limits a(x) and b(x) in terms of functions operating on the x-value on the stack.

\documentclass[pstricks, border=5pt]{standalone}
\usepackage{pst-func}

\makeatletter
\def\psVarIntegral{\pst@object{psVarIntegral}}
\def\psVarIntegral@i#1#2#3#4#5{%
  \begin@OpenObj%
  \addto@pscode{
    /xStart #1 def
    /xCurr #1 def
    /dx #2 #1 sub \psk@plotpoints\space div def
    /a #1 #3 def
    /b #1 #4 def
    /scx { \pst@number\psxunit mul } def
    /scy { \pst@number\psyunit mul } def
    tx@FuncDict begin /SFunc { #5 } def end
    \psk@plotpoints 1 add {
      a b \psk@Simpson
      tx@FuncDict begin Simpson I end
      scy xCurr scx exch xCurr xStart eq {moveto}{lineto}ifelse
      /xCurr xCurr dx add def
      /a xCurr #3 def
      /b xCurr #4 def
    } repeat
  }%
  \end@OpenObj%
}
\makeatother

\begin{document}

\psset{xunit=0.8,yunit=1.5}
\begin{pspicture}(-7,-2)(7,2)
  \psVarIntegral[plotpoints=500, linecolor=red]{-7}{7}{}{2 mul}{4 exp 1 add sqrt 4 exch div}
  \psaxes[Dy=0.5, arrows=->](0,0)(-7,-2)(7,2)
  \rput[rt](7,2){$f:x\mapsto \displaystyle\int_x^{2x} \frac{4}{\sqrt{1+t^4}}\, \textrm{d}t$}
\end{pspicture}
\end{document}

Answer 3

Here is another, quite compact, PSTricks solution.

The TikZ solution using the same numerical approach is given below to satisfy the OP.

\pstODEsolve (RKF45 method) from the pst-ode package is used to solve the definite integral between x and 2 x at each of the 501 plot points in the interval [-7,7]. The initial value for each \pstODEsolve invocation is set to zero to immediately get the definite integral at 2 x.

\documentclass[pstricks,border=5pt]{standalone}
\usepackage{pst-ode,pst-plot}

\begin{document}
\pstVerb{/result {} def} %initialise empty result list
\multido{\nX=-7.00+0.028}{501}{% 501 plotpoints
  %integral = [x 0 2x F(2x)] (two output points)-------------v   v----initial value
  \pstODEsolve[algebraicAll]{integral}{t | y[0]}{\nX}{2*\nX}{2}{0.0}{4/sqrt(1+t^4)}
  %append [x F(2x)] to results list
  \pstVerb{/result [result integral exch pop exch pop] cvx def}
}
%plot result
\psset{xunit=0.8,yunit=1.5}
\begin{pspicture}(-7,-2)(7,2)
  \psaxes[Dy=0.5, arrows=->](0,0)(-7,-2)(7,2)
  \listplot[linecolor=red]{result}
  \rput[rt](7,2){$f:x\mapsto \displaystyle\int_x^{2x} \frac{4}{\sqrt{1+t^4}}\, \textrm{d}t$}
\end{pspicture}
\end{document}

---

TikZ/PGFPlots solution, requires pdflatex --shell-escape:

\documentclass[tikz,border=5pt]{standalone}
\usepackage{pgfplots} \pgfplotsset{width=\linewidth,compat=1.9}
\usepackage{filecontents}

\begin{filecontents}{xyz.tex}
  \input pst-ode \input multido
  \pstVerb{/statefile (result.dat) (w) file def}
  \multido{\nX=-7.00+0.028}{501}{% 501 plotpoints
    \pstODEsolve[algebraicAll]{integral}{t | y[0]}{\nX}{2*\nX}{2}{0}{4/sqrt(1+t^4)}
    \pstVerb{[integral exch pop exch pop] tx@odeDict begin writeresult end}
  }
  \pstVerb{statefile closefile} \bye
\end{filecontents}
\immediate\write18{tex xyz}\immediate\write18{dvips xyz}
\immediate\write18{ps2pdf -dNOSAFER xyz.ps}

\begin{document}
\begin{tikzpicture}
  \begin{axis}[
      axis x line=center, axis y line=center, unit vector ratio=0.8 1.5,
      ymin=-2, ymax=2, xtick={-7,...,6}, ytick={-2,-1.5,...,1.5},
      y tick label style={/pgf/number format/.cd, fixed, fixed zerofill, precision=1},
    ]
    \addplot[red] table {result.dat};
    \node [anchor=north east] at (axis cs:7,2)
      {$f:x\mapsto \displaystyle\int_x^{2x} \frac{4}{\sqrt{1+t^4}}\, \textrm{d}t$};
  \end{axis}
\end{tikzpicture}
\end{document}

Answer 4

You can use the GNU Scientific Library (GSL) via the FFI of LuaJITTeX (and LuaTeX ≥ 1.0.3). Needs --shell-escape.

\documentclass{article}

\usepackage{pgfplots}
\pgfplotsset{compat=newest}

\usepackage{luacode}
\begin{luacode*}
local ffi = require("ffi")

ffi.cdef[[
typedef double (*gsl_callback) (double x, void * params);

typedef struct {
  gsl_callback F;
  void * params;
} gsl_function;

typedef void gsl_integration_workspace;

gsl_integration_workspace * gsl_integration_workspace_alloc (size_t n);

void gsl_integration_workspace_free (gsl_integration_workspace * w);

int gsl_integration_qags(gsl_function * f, double a, double b, double epsabs, double epsrel, size_t limit,
                         gsl_integration_workspace * workspace, double * result, double * abserr);
]]

local gsl = ffi.load("gsl")

function gsl_qags(f, a, b, epsabs, epsrel, limit)
   local limit = limit or 50
   local epsabs = epsabs or 1e-8
   local epsrel = epsabs or 1e-8

   local gsl_function = ffi.new("gsl_function")
   gsl_function.F = ffi.cast("gsl_callback", function(x, params) return f(x) end)
   gsl_function.params = nil

   local result = ffi.new('double[1]')
   local abserr = ffi.new('double[1]')

   local workspace = gsl.gsl_integration_workspace_alloc(limit)
   gsl.gsl_integration_qags(gsl_function, a, b, epsabs, epsrel, limit, workspace, result, abserr)
   gsl.gsl_integration_workspace_free(workspace)

   gsl_function.F:free()

   return result[0]
end

function f(x)
    tex.sprint(gsl_qags(function(t) return 4/math.sqrt(1+t^4) end, x, 2*x))
end
\end{luacode*}

\begin{document}

\begin{tikzpicture}[
  declare function={f(\x) = \directlua{f(\x)};}
  ]
  \begin{axis}[
    use fpu=false, % very important!
    no marks, samples=101,
    ]
    \addplot {f(x)};
  \end{axis}
\end{tikzpicture}

\end{document}