# Lesson 2: Floating point computations
# ===========================
#
# As we saw in Lesson 1, Maple will write decimals or floating-point numbers
# with a certain number of significant digits, specified by the variable "Digits".
# Here are some floating-point numbers, with the default Digits = 10.
> evalf(exp(20)), evalf(exp(-7)), evalf(exp(-20));
9 -8
.4851651954*10 , .0009118819656, .2061153622*10
# In principle a floating point number should be written as
# (possibly -) .digits * 10^exponent, but Maple doesn't use the exponent if it is
# small. The exponent can be anywhere from -9989 to 10009 for a floating
# point number to be written this way.
> bignumber:=evalf(.10*10^10009);
10009
bignumber := .1000000000*10
> smallnumber:= evalf(.10*10^(-9989));
-9989
smallnumber := .1000000000*10
# That should be a big enough range for most purposes, but actually there is
# another notation that is used for exponents outside that range, and these can be
# used for some REALLY big or small numbers.
> bignumber*10, bignumber^1000, smallnumber^10000;
10010 10008001 -99899999
.1000000000*10 , .1000000000*10 , .1000000000*10
# Many implementations of floating point actually store and calculate numbers
# in base 2 (converting to base 10 for output), but Maple actually uses base 10.
# So this calculation has no roundoff error, because 1/5 can be represented
# exactly in base 10:
> evalf(1/5)*5;
1.000000000
# This one, on the other hand, does have roundoff error, because 1/3 can't be
# represented exactly.
> evalf(1/3)*3;
.9999999999
# You can think of a floating point number as representing a whole interval of
# numbers, all of which have the same floating point representation. Thus
# everything from 1.9999999995 to 2.0000000005 is represented as 2.000000000.
# When you add a small floating point number to a larger number, some
# significant digits of the small number are lost.
> .1234567891*10^(-5)+1;
1.000001235
# In an extreme case, adding two numbers that differ by at least "Digits" orders
# of magnitude, the smaller number is lost altogether.
> 1.+ 10^(-10);
1.000000000
# This can have serious consequences if you are trying to add up many numbers
# that are each very small (e.g. in numerical integration). Suppose we start with
# 0 and add 10000 numbers, each equal to 4.444444444 * 10^(-4). The correct
# answer would be 4.444444444.
# I'll use a "do loop": the command total:= total + 4.444444444 * 10^(-4) will be
# executed
# 10000 times, once for each integer value of j from 1 to 10000. Any sequence of
# commands
# could go between the "do" and the "od" (which is "do" backwards). Note the
# ":" instead
# of ";", because we don't want to see 10000 outputs.
> total:= 0:
> for j from 1 to 10000 do
> total:= total + 4.444444444 * 10^(-4) od:
> total;
4.444440910
# The way to fix this would be to increase Digits: 14 should be enough.
> Digits:= 14:
> total:= 0:
> for j from 1 to 10000 do
> total:= total + 4.444444444 * 10^(-4) od:
> total;
4.4444444440000
> Digits:= 10:
# Another case where roundoff error can be severe is when we subtract two
# numbers that are almost equal. This is called "catastrophic cancellation".
> 1.000000001 - 1.000000000;
-8
.1*10
# We subtracted two numbers with 10 significant digits and the result has only
# one significant digit.
# One good way to think about this is in terms of absolute and relative errors. In
# approximating a "true value" T by an approximate value A, the absolute error
# is |T-A|. The relative error is |T-A|/|T|. In representing a number in floating
# point the relative error should be at most 5 * 10^(-Digits).
# When two numbers are multiplied or divided, the relative error of the result is
# at most (approximately) the sum of the two relative errors. This usually
# doesn't cause any problems.
# When two numbers are added or subtracted, the absolute error of the result is
# at most the sum of the two absolute errors. If the result is close to 0, this
# means the relative error can be much larger than the relative errors of the
# original numbers.
# Here is an example of catastrophic cancellation affecting the solution of a
# quadratic equation. Here's the solution formula from high school.
> r:= solve(a*x^2+b*x+c=0, x);
2 1/2 2 1/2
- b + (b - 4 a c) - b - (b - 4 a c)
r := 1/2 ---------------------, 1/2 ---------------------
a a
# Let's try it when b is large compared to a and c.
> s:=solve(x^2 + 9751*x + 1 = 0, x);
1/2 1/2
s := - 9751/2 + 1/2 95081997 , - 9751/2 - 1/2 95081997
# Let's use "evalf" to get floating point answers.
# There's no accuracy problem with the second one. We get 10 significant digits.
> evalf(s[2]);
-9750.999898
# The first, however, involves subtracting two numbers that are very close
# together. The result has about the same absolute error as for the second
# solution, but only 3 significant digits.
> evalf(s[1]);
-.000102
# Using a higher setting of "Digits" gives us more accuracy.
> evalf(s[1], 20); evalf(s[2], 20);
-.0001025535853263
-9750.9998974464146737
# But we can actually get better relative error at the same setting of "Digits" by
# using a different formula that won't involve this catastrophic cancellation.
# Recall that the product of the roots is the constant term c.
> expand(r[1]*r[2]);
c/a
# So the first root should be
> r[1] = c/r[2];
2 1/2
- b + (b - 4 a c) c a
1/2 --------------------- = 2 ---------------------
a 2 1/2
- b - (b - 4 a c)
# which won't produce a catastrophic cancellation (when b is positive). In our
# case:
> evalf(1/s[2]);
-.0001025535853
# In addition to "solve" which finds exact solutions (when possible), Maple has
# "fsolve" which finds floating-point approximations to solutions of equations,
# even those that "solve" can't handle. The syntax is similar to that of "solve".
# Let's see how it handles our quadratic:
> fsolve(x^2 + 9751*x + 1 = 0);
-9750.999897, -.0001025535853
# It won't work, of course, on equations that have symbolic parameters.
> fsolve(a*x^2 + b*x + c=0, x);
Error, (in fsolve) should use exactly all the indeterminates
# How about our fifth-degree polynomial from last time?
> fsolve(x^5 - x^3 + 2*x + 1);
-.5603499931
# For a polynomial, "fsolve" finds all the real roots. This polynomial has only
# one. There should be four other complex roots. To tell "fsolve" to find both
# real and complex roots, you give it a third input, the word "complex".
> s:=fsolve(x^5 - x^3 + 2*x + 1, x, complex);
s := - .7961643382 - .6551644381 I, - .7961643382 + .6551644381 I,
-.5603499931, 1.076339335 - .7212070670 I, 1.076339335 + .7212070670 I
# Should we trust "fsolve" here? One way to check a solution is to substitute the
# answer into the original problem. The "subs" command does this.
> subs(x = s[1], x^5 - x^3 + 2*x + 1);
-9 -9
.6*10 + .4*10 I
# Not bad. You can't expect more accuracy than that at Digits=10. To get more
# accuracy, of course, you can use a higher "Digits" setting. Now let's try a
# non-polynomial equation that stumped "solve".
> fsolve(sin(x)=1-x, x);
.5109734294
# This one has only one real solution (since x + sin(x) is an increasing function).
# But even if there were more, "fsolve" would only give us one solution for a
# non-polynomial equation.
> fsolve((x^2-1)/(x+2), x);
-1.000000000
# For the polynomial x^2-1 itself, "fsolve" would give us both solutions.
> fsolve(x^2-1, x);
-1.000000000, 1.
# Of course the solutions of (x^2-1)/(x+2)=0 are the two solutions of x^2-1=0,
# but since the expression was not a polynomial "fsolve" only gave us one
# solution.
# Let's try something more complicated.
> f:= x - 10*sin(x);
f := x - 10 sin(x)
# Before trying "fsolve", let's get a picture of what happens here by plotting a
# graph.
# The "plot" command, in its simplest form, takes two inputs, an expression f(v)
# depending
# on a variable v, and a range for the v axis in the form v = a .. b. This plots the
# curve
# y = f(v) for v from a to b. The scale for the y axis is chosen automatically. If
# you prefer,
# you can set it yourself by including another input y = c .. d.
> plot(f, x=-15 .. 15);
# The plot comes up in a separate window, which can be copied (in the Edit
# menu of the plot window) and pasted into the worksheet (in the Edit menu of
# the worksheet).
#
** Maple V Graphics **
# There appear to be 7 real solutions of f=0, at about x=-8.5, -7, -3, 0, 3, 7 and
# 8.5. In general you'd have to worry about solutions outside the interval we
# plotted, but here It's not hard to show there are no solutions outside the
# interval [-10, 10], because -1 <= sin(x) <= 1 means x-10 <= f <= x+10.
# By the way, there's a very useful feature of these plot windows: if you click the
# left mouse button with the mouse pointer in the window, Maple shows you the
# X and Y coordinates of the tip of the arrow. This lets you locate any point on
# the plot with better precision than you would have guessing "by eye", especially
# for points not on the axes.
> fsolve(f);
2.852341894
# This is the solution near 3. To get other solutions, you can use an option in
# "fsolve" to find roots in a given interval.
> fsolve(f, x=6 .. 8);
7.068174358
> fsolve(f, x=8 .. 9);
8.423203932
# The other roots can be obtained by symmetry (f is an odd function).
#
#
