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.
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))
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;
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.
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.