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:

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.

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:

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:

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:

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 .

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.