TransWikia.com

Calculating definite integral of smooth function in pgfplots

TeX - LaTeX Asked on December 12, 2020

I have the following code which plots the graph of (1 - x^2)^n for n=1,3,5,7. Is there a way to calculate the integrals from -1 to 1 directly in pgfplots? The goal is to have normalized functions.

documentclass{standalone}
usepackage{pgfplots}
pgfplotsset{
  compat=1.17,
  compat/show suggested version=false,
}

pgfmathdeclarefunction{kn}{1}{%
  % should be normalized dividing it by its integral from -1 to 1
  pgfmathparse{(1 - x^2)^#1}%
}

begin{document}
begin{tikzpicture}
  begin{axis}[
      axis lines=center,
      xlabel={$x$},
      ylabel={$y$},
      xmin=-1.2, xmax=1.2,
      ymin=-0.2, ymax=5,
      xtick={-1,1},
      ytick={1},
      every axis plot/.append style={
        smooth,
        domain=-1:1,
      },
    ]

    addplot [red] {kn(1)};
    addplot [blue] {kn(3)};
    addplot [yellow] {kn(5)};
    addplot [green] {kn(7)};

  end{axis}
end{tikzpicture}
end{document}

enter image description here

3 Answers

Accurate solution based on expl3

In this solution, each of the needed integrals is precomputed exactly once using the very accurate l3fp engine and 40 rectangles for the midpoint rule method (the number of rectangles is just a parameter and can be freely changed).

documentclass[tikz, border=2mm]{standalone}
usepackage{xparse}
usepackage{pgfplots}
pgfplotsset{
  compat=1.17,
  compat/show suggested version=false,
}

ExplSyntaxOn
fp_new:N l__noibe_result_fp
fp_new:N l__noibe_currentx_fp
fp_new:N l__noibe_deltax_fp

% Compute an approximation of the integral of a function over an interval
% using the midpoint rule.
%
% Arguments: macro or tl var for storing the result, unary function, interval
% start, interval end, number of rectangles
cs_new_protected:Npn noibe_set_to_midpoint_rule_riemann_sum:NNnnn #1#2#3#4#5
  {
    fp_zero:N l__noibe_result_fp
    fp_set:Nn l__noibe_deltax_fp { (#4 - #3) / (#5) }
    fp_set:Nn l__noibe_currentx_fp { #3 + 0.5*l__noibe_deltax_fp }

    int_step_inline:nn {#5}
      {
        fp_add:Nn l__noibe_result_fp { #2 { l__noibe_currentx_fp } }
        fp_add:Nn l__noibe_currentx_fp { l__noibe_deltax_fp }
      }

    tl_set:Nx #1 { fp_eval:n { l__noibe_deltax_fp * l__noibe_result_fp } }
  }

cs_generate_variant:Nn noibe_set_to_midpoint_rule_riemann_sum:NNnnn { c }

% Macro name stem for results, parameter, nb rectangles
cs_new_protected:Npn noibe_compute_kn_integral:nnn #1#2#3
  {
    cs_set:Npn noibe__tmp_function:n ##1 { (1 - (##1)^2)^(#2) }
    noibe_set_to_midpoint_rule_riemann_sum:cNnnn { #1 int_to_roman:n {#2} }
      noibe__tmp_function:n { -1 } { 1 } {#3}
  }

% Document-level interface
NewDocumentCommand computeKnIntegral { m m m }
  {
    noibe_compute_kn_integral:nnn {#1} {#2} {#3}
  }
ExplSyntaxOff

% Compute the integrals for parameters 1, 3, 5, 7
pgfplotsinvokeforeach{1, 3, 5, 7}{%
  computeKnIntegral{knIntegral}{#1}{40}% 40 is the number of rectangles
}

% Declare a kn function with two arguments: the parameter and the variable ('x')
pgfmathdeclarefunction{kn}{2}{%
  begingroup
    pgfmathfloatparsenumber{#1}%
    pgfmathfloattoint{pgfmathresult}%
    edeftheKnIntegral{%
      csname knIntegralromannumeralpgfmathresultspaceendcsname}%
    pgfmathparse{ (1 - (#2)^2)^(#1) / theKnIntegral }%
    pgfmathsmugglepgfmathresult
  endgroup
}

begin{document}
begin{tikzpicture}
  begin{axis}[
      axis lines=center,
      xlabel={$x$},
      ylabel={$y$},
      enlarge x limits=0.1,
      enlarge y limits=auto,
      every axis plot/.append style={smooth, domain=-1:1},
    ]

    addplot [red] {kn(1, x)};
    addplot [blue] {kn(3, x)};
    addplot [yellow] {kn(5, x)};
    addplot [green] {kn(7, x)};

  end{axis}
end{tikzpicture}
end{document}

enter image description here

Solution based on pgfmath

The following solution uses pgfmath with its fpu library to compute the integrals (exactly once each). I use only 20 rectangles here not because of slowness, but because the fpu engine of pgfmath is not very accurate and I don't want to accumulate too many errors due to a high number of operations (with this engine, the number of significant digits is quite small compared to what l3fp—the engine used to compute the integrals in the first solution—offers).

There is a commented-out code path that provides a workaround in case you encounter a pgfmath error message mentioning @@str@@:. I needed this workaround a couple of days ago, but after todays's upgrade of my TeX Live packages (from Debian unstable), it appears not to be necessary anymore (and even causes an error). So, only enable this workaround if you get the error.

documentclass[tikz, border=2mm]{standalone}
usepackage{etoolbox}
usepackage{pgfplots}
pgfplotsset{
  compat=1.17,
  compat/show suggested version=false,
}
usepgflibrary{fpu}

makeatletter

% Workaround for a problem I had before the last update of my TeX Live
% packages (Debian unstable). Uncomment the definition if you have an error
% message mentioning '@@str@@:'
% newcommand*{my@decode@fpu@string@argument}[2]{%
%   begingroup
%     letpgfmath@basic@stack@push@operand@firstofone
%     edefmy@tmp{%
%       endgroupdefnoexpand#2%
%         {unexpandedexpandafterexpandafterexpandafter{%
%            pgfmathfloat@stack@push@operand@single@str #1relax}}}%
%   my@tmp
% }

% Compute an approximation of the integral of a function over an interval
% using the midpoint rule.
%
% Arguments: function (prefixed with pgfmath@fpu@stringmarker), x_min, x_max,
%            number of rectangles.
pgfmathdeclarefunction{midrule}{4}{%
  begingroup
    pgfset{fpu=true}%
    pgfmathsetmacro{my@result}{0}%
    pgfmathsetmacro{my@delta@x}{((#3) - (#2)) / (#4)}%
    pgfmathsetmacro{my@x}{(#2) + 0.5*my@delta@x}%
    % If you have an error message mentioning '@@str@@:', uncomment this line
    % and comment out the following 'defmy@funcname{#1}' line:
    % expandaftermy@decode@fpu@string@argumentexpandafter{#1}{my@funcname}%
    defmy@funcname{#1}%
    %
    pgfplotsforeachungrouped x in {1,...,#4}{%
      pgfmathsetmacro{my@result}{my@result + my@funcname(my@x)}%
      pgfmathsetmacro{my@x}{my@x + my@delta@x}%
    }%
    %
    pgfmathparse{my@delta@x * my@result}%
    pgfset{fpu=false}%
    pgfmathfloattofixed{pgfmathresult}%
    pgfmathsmugglepgfmathresult
  endgroup
}

newcommand*{defineknForParam}[2]{%
  pgfmathdeclarefunction{#1#2}{1}{%
    pgfmathparse{(1 - (##1)^2)^(#2)}%
  }%
}

% Define functions kn1, kn3, kn5 and kn7.
pgfplotsinvokeforeach{1, 3, 5, 7}{%
  defineknForParam{knbase}{#1}%
  % Compute and store the integral corresponding to parameter #1. 20 is the
  % number of rectangles used for the midpoint rule.
  pgfmathmidrule{"knbase#1"}{-1}{1}{20}%
  csedef{knIntegralromannumeral #1space}{pgfmathresult}%
  %
  pgfmathdeclarefunction{kn#1}{1}{%
    pgfmathparse{ knbase#1(##1) / csuse{knIntegralromannumeral #1space} }%
  }%
}
makeatother

begin{document}
begin{tikzpicture}
  begin{axis}[
      axis lines=center,
      xlabel={$x$},
      ylabel={$y$},
      enlarge x limits=0.1,
      enlarge y limits=auto,
      every axis plot/.append style={smooth, domain=-1:1},
    ]

    addplot [red] {kn1(x)};
    addplot [blue] {kn3(x)};
    addplot [yellow] {kn5(x)};
    addplot [green] {kn7(x)};

  end{axis}
end{tikzpicture}
end{document}

enter image description here

Correct answer by frougon on December 12, 2020

A sagetex solution. I've grabbed the code from a previous answer here and modified it to fit your question. The code could certainly be simplified by removing lines giving you options on what the graph should look like.

documentclass{standalone}
usepackage[usenames,dvipsnames]{xcolor}
usepackage{pgfplots}
usepackage{sagetex}
usetikzlibrary{spy}
usetikzlibrary{backgrounds}
usetikzlibrary{decorations}
pgfplotsset{compat=newest}% use newest version
begin{document}
begin{sagesilent}
####### SCREEN SETUP #####################
LowerX = -1.0
UpperX = 1.0
LowerY = -0.2
UpperY = 1.8
step = .01
Scale = 1.0
xscale=1.0
yscale=1.0
#####################TIKZ PICTURE SET UP ###########
output = r""
output += r"begin{tikzpicture}"
output += r"[line cap=round,line join=round,x=8.75cm,y=8cm]"
output += r"begin{axis}["
output += r"grid = none,"
#Change "both" to "none" in above line to remove graph paper
output += r"minor tick num=4,"
output += r"every major grid/.style={Red!30, opacity=1.0},"
output += r"every minor grid/.style={ForestGreen!30, opacity=1.0},"
output += r"height= %ftextwidth,"%(yscale)
output += r"width = %ftextwidth,"%(xscale)
output += r"thick,"
output += r"black,"
output += r"axis lines=center,"
#Comment out above line to have graph in a boxed frame (no axes)
output += r"domain=%f:%f,"%(LowerX,UpperX)
output += r"line join=bevel,"
output += r"xmin=%f,xmax=%f,ymin= %f,ymax=%f,"%(LowerX,UpperX,LowerY, UpperY)
#output += r"xticklabels=empty,"
#output += r"yticklabels=empty,"
output += r"major tick length=5pt,"
output += r"minor tick length=0pt,"
output += r"major x tick style={black,very thick},"
output += r"major y tick style={black,very thick},"
output += r"minor x tick style={black,thin},"
output += r"minor y tick style={black,thin},"
#output += r"xtick=empty,"
#output += r"ytick=empty"
output += r"]"
##############FUNCTIONS#################################
##FUNCTION 1
t1 =  var('t1')
const1 = numerical_integral(1-x^2, -1, 1, max_points=100)
x1_coords = srange(LowerX,UpperX,step)
y1_coords = [((1-t1^2)/const1[0]).n(digits=6) for t1 in x1_coords]
output += r"addplot[thin, NavyBlue, unbounded coords=jump] coordinates {"
for i in range(0,len(x1_coords)):
    if (y1_coords[i])<LowerY or (y1_coords[i])>UpperY:
        output += r"(%f,inf) "%(x1_coords[i])
    else:
        output += r"(%f,%f) "%(x1_coords[i],y1_coords[i])
output += r"};"
##FUNCTION 2 #########################################
t2 =  var('t2')
const2 = numerical_integral((1-x^2)^3, -1, 1, max_points=100)
x2_coords = srange(LowerX,UpperX,step)
y2_coords = [((1-t2^2)^3/const2[0]).n(digits=6) for t2 in x2_coords]
output += r"addplot[thin, Orchid, unbounded coords=jump] coordinates {"
for i in range(0,len(x2_coords)):
    if (y2_coords[i])<LowerY or (y2_coords[i])>UpperY:
        output += r"(%f,inf) "%(x2_coords[i])
    else:
        output += r"(%f,%f) "%(x2_coords[i],y2_coords[i])
output += r"};"
##FUNCTION 3 ##############################################
t3 =  var('t3')
const3 = numerical_integral((1-x^2)^5, -1, 1, max_points=100)
x3_coords = srange(LowerX,UpperX,step)
y3_coords = [((1-t3^2)^5/const3[0]).n(digits=6) for t3 in x3_coords]
output += r"addplot[thin, Peach, unbounded coords=jump] coordinates {"
for i in range(0,len(x3_coords)):
    if (y3_coords[i])<LowerY or (y3_coords[i])>UpperY:
        output += r"(%f, inf) "%(x3_coords[i])
    else:
        output += r"(%f, %f) "%(x3_coords[i],y3_coords[i])
output += r"};"
##FUNCTION 3 ##############################################
t4 =  var('t4')
const4 = numerical_integral((1-x^2)^7, -1, 1, max_points=100)
x4_coords = srange(LowerX,UpperX,step)
y4_coords = [((1-t4^2)^5/const4[0]).n(digits=6) for t4 in x4_coords]
output += r"addplot[thin, ForestGreen, unbounded coords=jump] coordinates {"
for i in range(0,len(x3_coords)):
    if (y4_coords[i])<LowerY or (y4_coords[i])>UpperY:
        output += r"(%f, inf) "%(x4_coords[i])
    else:
        output += r"(%f, %f) "%(x4_coords[i],y4_coords[i])
output += r"};"
##### COMMENT OUT A LINE OF SAGESILENT BY STARTING WITH #
output += r"end{axis}"
output += r"end{tikzpicture}"
end{sagesilent}
sagestr{output}
end{document}

Running in Cocalc we get: enter image description here

The sagetex package requires a computer algebra system, SAGE, in order to work. Installing it and getting it to work nicely with LaTeX on Windows computers can be problematic at times, as @Benjamin McKay comments. A free Cocalc account avoids these problems as your work is done in the cloud. Cocalc performance has deteriorated a bit over the last few months but should be good enough for light work like this.

NOTE: The output looks different than your picture. I did a sanity check for n=1 getting an integral of 1-x^2 to be 4/3 over -1 to 1. The height of 1-x^2 at 0 is 1 and 1/(4/3) is 3/4.

The CTAN documentation on sagetex is here. The documentation for SAGE is here.

Answered by DJP on December 12, 2020

Another, accurate solution using PSTricks. It (ab)uses pstODEsolve (RKF45) for computing the definite integrals. enter image description here

latex+dvips+ps2pdf

documentclass[pstricks]{standalone}
usepackage{pst-ode,pst-plot,pstricks-add}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% I(n)=int_{-1}^1 (1-t^2)^n dt
% #1: n
% #2: PS variable for result I(n)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
defI(#1)#2{%     two output points are enough---v   v---y[0](-1) (initial value)
  pstODEsolve[algebraicAll]{#2}{y[0]}{-0.999}{1}{2}{0.0}{(1-t^2)^#1}
  %            integration interval t_a---^    ^---t_b
  %  From ret value `#2', we throw away initial value y(n,-1)
  pstVerb{/#2 #2 exch pop def}
}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

begin{document}
% compute and save the definite integrals to int*
I(1){int1}% n=1
I(3){int3}% n=3
I(5){int5}% n=5
I(7){int7}% n=7
%
begin{pspicture}(-0.4,-0.7)(0.5,5)
begin{psgraph}[xAxisLabel={$x$},yAxisLabel={$y$},linewidth=0.5pt,
    Dx=0.5,Dy=0.5, arrows=->](0,0)(-1.2,0)(1.2,1.7){6cm}{!} % x-y-axis with same unit
  psplot[linecolor=red,plotpoints=100,algebraic]{-1}{1}{ (1-x^2)^1 / int1 }
  psplot[linecolor=blue,plotpoints=100,algebraic]{-1}{1}{ (1-x^2)^3 / int3 }
  psplot[linecolor=yellow,plotpoints=100,algebraic]{-1}{1}{ (1-x^2)^5 / int5 }
  psplot[linecolor=green,plotpoints=100,algebraic]{-1}{1}{ (1-x^2)^7 / int7 }
end{psgraph}
end{pspicture}
end{document}

Answered by AlexG on December 12, 2020

Add your own answers!

Ask a Question

Get help from others!

© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP