4.8 Numerical integration and differentiation

Two special functions, int_dx() and diff_dx(), may be used to integrate or differentiate algebraic expressions numerically. In each case, the letter x is the dummy variable which is to be used in the integration or differentiation and may be replaced by any valid variable name of up to 16 characters in length.

The function int_dx() takes three parameters – firstly the expression to be integrated, which may optionally be placed in quotes, followed by the minimum and maximum integration limits. These may have any physical dimensions, so long as they match, but must both be real numbers. For example, the following would plot the integral of the function $\sin (x)$ :

plot int_dt('sin(t)',0,x)

The function diff_dx() takes two obligatory parameters plus one further optional parameter. The first is the expression to be differentiated, which, as above, may optionally placed in quotes for clarity. This should be followed by the numerical value $x$ of the dummy variable at the point where the expression is to be differentiated. This value may have any physical dimensions, and may be a complex number if complex arithmetic is enabled. The final, optional, parameter to the diff_dx() function is an approximate step size, which indicates the range of argument values over which Pyxplot should take samples to determine the gradient. If no value is supplied, a value of $10^{-6}x$ is used, replaced by $10^{-6}$ if $x=0$ . The following example would evaluate the differential of the function $\cos (x)$ with respect to $x$ at $x=1.0$ :

$\includegraphics[width=0.9cm]{tick.eps}$

print diff_dx(’cos(x)’, 1.0)

When complex arithmetic is enabled, Pyxplot checks that the function being differentiated satisfies the Cauchy-Riemann equations, and returns an error if it does not, to indicate that it is not differentiable. The following is an example of a function which is not differentiable, and which throws an error because the Cauchy-Riemann equations are not satisfied:

$\includegraphics[width=0.9cm]{cross.eps}$

set num comp print diff_dx(Re(sin(x)),1)

Advanced users may be interested to know that int_dx() function is implemented using the gsl_integration_qags() function of the Gnu Scientific Library (GSL), and the diff_dx() function is implemented using the gsl_deriv_central() function of the same library. Any caveats which apply to the use of these routines also apply to Pyxplot’s numerical calculus.

Example: Integrating the function $\mathrm{sinc}(x)$

The function $\mathrm{sinc}(x)$ cannot be integrated analytically, but it can be shown that

$\int _0^{\pm \infty } \mathrm{sinc}(x)\, \mathrm{d}x = \pm \pi /2 .$

In the following script, we use Pyxplot’s facilities for numerical integration to produce a plot of

$y=\int _0^{x} \mathrm{sinc}(x)\, \mathrm{d}x .$

We reduce the number of samples taken along the abscissa axis to 80, as evaluation of the numerical integral may be time consuming on older computers. We use the set xformat command (see Section 8.8.8) to demark both the x- and y-axes in fractions of $\pi$ :
set samples 80
set key bottom right
set xformat r"%s$pi$"%(x/pi)
set yformat r"%s$pi$"%(y/pi)
set xrange [-5*pi:5*pi]
plot int_dz(sinc(z),0,x)
$\includegraphics[width=\textwidth ]{examples/eps/ex_integration}$