☰

基于的精确解决方案expl3

tag-icon • tag-icon tikz-pgf pgfplots pgfmath
基于的精确解决方案expl3

(1 - x^2)^n我有以下代码，用于绘制的图形n=1,3,5,7。有没有办法直接在 pgfplots 中计算从 -1 到 1 的积分？目标是拥有规范化函数。

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

在此处输入图片描述

答案1

基于的精确解决方案expl3

在这个解决方案中，每个需要的积分都使用非常精确的l3fp引擎和 40 个矩形（矩形的数量只是一个参数，可以自由更改）预先计算一次。

\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}%
    \edef\theKnIntegral{%
      \csname knIntegral\romannumeral\pgfmathresult\space\endcsname}%
    \pgfmathparse{ (1 - (#2)^2)^(#1) / \theKnIntegral }%
    \pgfmathsmuggle\pgfmathresult
  \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}

在此处输入图片描述

解决方案基于pgfmath

以下解决方案使用pgfmath其fpu库来计算积分（每个积分恰好一次）。我在这里只使用了 20 个矩形，不是因为速度慢，而是因为引擎fpu不是pgfmath很准确，而且我不想由于操作次数过多而积累太多错误（使用此引擎，有效数字的数量与l3fp第一个解决方案中用于计算积分的引擎所提供的相比非常少）。

有一个注释掉的代码路径，它提供了一种解决方法，以防您遇到pgfmath提及的错误消息@@str@@:。几天前我需要这个解决方法，但今天升级我的 TeX Live 软件包（从 Debian 不稳定版升级）后，它似乎不再需要（甚至会导致错误）。因此，只有在收到错误时才启用此解决方法。

\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
%     \let\pgfmath@basic@stack@push@operand\@firstofone
%     \edef\my@tmp{%
%       \endgroup\def\noexpand#2%
%         {\unexpanded\expandafter\expandafter\expandafter{%
%            \pgfmathfloat@stack@push@operand@single@str #1\relax}}}%
%   \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 '\def\my@funcname{#1}' line:
    % \expandafter\my@decode@fpu@string@argument\expandafter{#1}{\my@funcname}%
    \def\my@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}%
    \pgfmathsmuggle\pgfmathresult
  \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{knIntegral\romannumeral #1\space}{\pgfmathresult}%
  %
  \pgfmathdeclarefunction{kn#1}{1}{%
    \pgfmathparse{ knbase#1(##1) / \csuse{knIntegral\romannumeral #1\space} }%
  }%
}
\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}

在此处输入图片描述

答案2

一个sagetex解决方案。我从之前的答案中获取了代码这里并对其进行了修改以适合您的问题。代码当然可以通过删除为您提供图表外观选项的行来简化。

\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= %f\textwidth,"%(yscale)
output += r"width = %f\textwidth,"%(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}

在 Cocalc 中运行我们得到： 在此处输入图片描述

该sagetex软件包需要计算机代数系统 SAGE 才能工作。正如 @Benjamin McKay 所说，在 Windows 计算机上安装它并使其与 LaTeX 完美配合有时会出现问题。免费可钙帐户可以避免这些问题，因为您的工作是在云端完成的。Cocalc 的性能在过去几个月里有所下降，但对于像这样的轻量级工作来说应该足够好了。

注意：输出结果与您的图片不同。我对 n=1 进行了完整性检查，得到 1-x^2 的积分为 4/3，除以 -1 到 1。1-x^2 在 0 处的高度为 1，1/(4/3) 为 3/4。

CTANsagetex文档这里. SAGE 的文档是这里。

答案3

另一个准确的解决方案是使用 PSTricks。它（滥用）使用\pstODEsolve（RKF45) 来计算定积分。 在此处输入图片描述

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)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\def\I(#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}

相关内容