This is documentation for Mathematica 5, which was
based on an earlier version of the Wolfram Language.
View current documentation (Version 11.1)

Documentation / Mathematica / Add-ons & Links / Standard Packages / Introduction /

Numerical Mathematics Packages

The standard add-on packages for numerical mathematics extend the set of built-in numerical functions, and take advantage of the arbitrary-precision arithmetic that is the foundation of numerical computations in Mathematica. The NumericalMath packages provide fitting functions (polynomial, spline, trigonometric), numerical versions of some of the kernel functions (ND, NLimit, NResidue, NSeries), numerical integration functions (CauchyPrincipalValue, ListIntegrate, NIntegrateInterpolationFunction), support for numerical solution of differential equations (BesselZeros, Butcher, OrderStar), alternatives to FindRoot using interpolation or interval methods, functions for approximating by a ratio of polynomials, and pedagogical functions for exploring floating-point arithmetic and numerical quadrature.

This initializes the numerical mathematics packages so that they load as needed.

In[1]:= <<NumericalMath`

Here is a set of points that lie approximately on a circle.

In[2]:= (data = Table[(x = Random[Real, {-1, 1}];
Table[Random[Real, {-.1, .1}], {2}] +
{x, Sign[.5 - Random[]] Sqrt[1 - x^2]}), {15}];
ListPlot[data, AspectRatio -> 1,
Ticks -> {{-1., -.5, .5, 1.}, Automatic}])


The points given to SplineFit are first sorted according to polar angle. The resulting spline interpolates the points.

In[3]:= (data = Map[Last,
Sort[Map[{Apply[ArcTan, #], #}&, data]]];
spline = SplineFit[Join[data, {First[data]}], Cubic];
ParametricPlot[spline[u], {u, 0, 15},
AspectRatio -> 1, Compiled -> False,
Ticks -> {{-1., -.5, .5, 1.}, Automatic}])


You can use GeneralMiniMaxApproximation to approximate the inverse to the Euler gamma function on an interval. This shows the gamma function on .

In[4]:= Plot[Gamma[t], {t, 1.5, 6},
AxesOrigin -> {1.5, 0}, AspectRatio ->1,
AxesLabel -> {"t", "Gamma[t]"}]


This gives a degree rational function that approximates the inverse on .

In[5]:= (gmma = GeneralMiniMaxApproximation[{Gamma[t], t},
{t, {1.5, 6}, 2, 4}, s];
inverseGamma = gmma[[2, 1]])


Here is a plot of the approximation. For a minimax approximation, the error is distributed over the interval.

In[6]:= (Plot[inverseGamma, {s, Gamma[1.5], Gamma[6]},
AxesOrigin -> {0, 1.5},
AspectRatio ->1, AxesLabel -> {"s", "inverseGamma[s]"}])


When you want to solve differential equations approximately, there are many considerations to take into account. You must ensure that the numerical solution process is sufficiently accurate (it is consistent with the differential system), that it actually solves the problem in which you are interested (it is convergent), and that any errors do not cause your approximate solution to drift far from the true solution (it is stable). Further, you may wish to ensure that your approximate solution samples (interpolates) the exact solution at specific points. Order stars provide a framework for the study and construction of numerical methods with these desirable properties.

This loads the function Pade for calculating Padé approximations.

In[7]:= <<Calculus`Pade`

Here the Padé approximation to represents the approximate solution of a differential equation, while represents the exact solution.

In[8]:= (f = (z^2*(1 + z^2) + (1 + z/2)*(1 - z - z^3))/
(z^2*(1 + z^2) + (1 - z/2)*(1 - z - z^3));
r = Pade[f, {z, 0, 1, 1}])


This is an order star plot for analyzing the numerical method represented by the Padé approximation.

In[9]:= OrderStar[r, f, OrderStarInterpolation -> True,
OrderStarLegend -> True, ImageSize -> {500, Automatic}]


The order star reflects the fact that the order of the approximant at the origin is , which is one less than the number of adjoining sectors there. The shaded region gives the region of growth of the approximant over the function and yields information on the accumulation of numerical errors.