Floating-Point Arguments and Function Sensitivity

Particular choices of parameters can reduce some special functions to simpler special functions, elementary functions, or numbers. Nevertheless, for most parameters MuPAD® returns the symbolic notation of a special function. In such cases, you can approximate the value of a special function numerically. To approximate a special function numerically, use the float command or call the special function with floating-point arguments.

When approximating the value of a special function numerically, remember that floating-point results can be extremely sensitive to numeric precision. Also, floating-point results are prone to roundoff errors. The following approaches can help you recognize and avoid incorrect results:

  • When possible, use symbolic computations. Switch to floating-point arithmetic only if you cannot obtain symbolic results. See Use Symbolic Computations When Possible.

  • Numeric computations are sensitive to the DIGITS environment variable that determines the numeric working precision. Increase the precision of numeric computations, and check if the result changes significantly. See Increase Precision.

  • Compute the value of a special function symbolically, and then approximate the result numerically. Also, compute the value of a special function using the floating-point parameters. Significant difference in these two results indicates that one or both approximations are incorrect. See Approximate Parameters and Approximate Results.

  • Plot the function. See Plot Special Functions.

Use Symbolic Computations When Possible

By default, MuPAD performs computations in exact symbolic form. For example, standard mathematical constants have their own symbolic representations in MuPAD. Using these representations, you can keep the exact value of the constant throughout your computations. You always can find a numeric approximation of a constant by using the float function:

pi := float(PI)

Avoid unnecessary conversions to floating-point numbers. A floating-point number approximates a constant; it is not the constant itself. Using this approximation, you can get incorrect results. For example, the heaviside special function returns different results for the sine of π and the sine of 10-digit floating-point approximation of π:

heaviside(sin(PI)), heaviside(sin(pi))

Increase Precision

The Riemann hypothesis states that all nontrivial zeros of the Riemann Zeta function have the same real part . To locate possible zeros of the Zeta function, plot its absolute value . The following plot shows the first three nontrivial roots of the Zeta function :

plot(abs(zeta(1/2 + I*y)), y = 0..30,
        AxesTitles = ["y", "|zeta|"])

Use the numeric solver to approximate the first three zeros of this Zeta function:

numeric::solve(zeta(1/2 + I*y), y = 13..15),
numeric::solve(zeta(1/2 + I*y), y = 20..22),
numeric::solve(zeta(1/2 + I*y), y = 24..26)

Now, consider the same function, but slightly increase the real part: . According to the Riemann hypothesis, this function does not have a zero for any real value y. By default, MuPAD uses 10 significant decimal digits for computations that involve floating-point numbers. When you use the numeric::solve solver with the default number of digits, the solver finds the following (nonexisting) zero of the Zeta function:

numeric::solve(zeta(1000000001/2000000000 + I*y), y = 14..15)

Increasing the numbers of digits shows that the result is incorrect. The Zeta function does not have a zero at 14 < y < 15:

DIGITS:=15:
numeric::solve(zeta(1000000001/2000000000 + I*y), y = 14..15)

delete DIGITS;

Approximate Parameters and Approximate Results

Bessel functions with half integer indices return exact symbolic expressions. Approximating these expressions by floating-point numbers, you can get very unstable results. For example, the exact symbolic expression for the following Bessel function is:

B := besselJ(53/2, PI)

Use the float command to approximate this expression numerically:

float(B)

Now, call the Bessel function with the floating-point parameter. Significant difference in these two approximations indicates that one or both results are incorrect:

besselJ(53/2, float(PI))

Increase the numeric working precision to obtain more accurate approximations:

DIGITS:= 45: float(B); besselJ(53/2, float(PI))

delete DIGITS;

Now you can see that using the floating-point parameter to compute the Bessel function produces the correct result (within working precision). Approximation of the exact symbolic expression for that Bessel function returns the wrong result because of numerical instability.

Plot Special Functions

Plotting the function can help you recognize incorrect floating-point approximations. For example, the numeric approximation of the following Bessel function returns:

B := besselJ(53/2, PI):
float(B)

Plot the function for the values of x around 53/2. The function plot shows that the floating-point approximation is incorrect:

plot(besselJ(x, PI), x = 26..27)

Sometimes, to see that the floating-point approximation is incorrect, you must zoom the particular parts of the function plot. For example, the numeric solver finds the unexpected zero of the Zeta function:

numeric::solve(zeta(1000000001/2000000000 + I*y), y = 14..15)

To investigate whether the Zeta function actually has a zero at that point or whether the result appears because of the roundoff error, plot the absolute value of Zeta function .

plot(abs(zeta(1000000001/2000000000 + I*y)), y = 0..30,
                          AxesTitles = ["y", "|zeta|"])

To see more details of the function plot near the possible zero, zoom the plot. To see that the numeric result is incorrect, enlarge that part of the function plot beyond the numeric working precision, and then reevaluate the plot. When you zoom and reevaluate, MuPAD recalculates the part of the function plot with the increased numeric precision.

    Note:   When zooming, MuPAD does not automatically reevaluate the function plot.

To get accurate results after zooming the plot, use the Recalculate button . After zooming and reevaluating the plot, you can see that the function does not have a zero at that interval.

Was this topic helpful?