h=0.00001
h=0.00000001
h=0.000000000001
h=0.00000000000001
h=0.000000000000000001
question id: h-too-small-1
\[ \newcommand{\dnorm}{\text{dnorm}} \newcommand{\pnorm}{\text{pnorm}} \newcommand{\recip}{\text{recip}} \]
Exercise 1 In mathematical theory, we speak of \(\lim_{h \rightarrow 0}\). But in numerical calculations, such as those done by numD()
, you can set \(h\) too small in the slope function. This will become obvious in the graph when \(h\) is too small.
Active R chunk lst-h-too-small calculates and plots a numerical derivative of the sinusoid function, using a specified value of \(h\). (That’s the argument named .h
.) Remake the plot, using every smaller .h
, until you start to see imperfections in the calculated values.
Which of these values for h
is the smallest you can go before the problems of too-small h
start showing up?
h=0.00001
h=0.00000001
h=0.000000000001
h=0.00000000000001
h=0.000000000000000001
question id: h-too-small-1
Exercise 2
In Figure fig-reptile-put-kitchen, which of the line segments is tangent to the curve at the point marked with a dot?
question id: reptile-put-kitchen-1
Exercise 3 Sometimes a bit of algebra can help us see what’s going on with the instantaneous rate of change. Consider the exponential function \(e^x\).
Rather than writing the slope function definition with a 0.1, let’s substitute in the symbol \(h\). This gives \[{\cal D}\exp(x) = \frac{e^{x+h} - e^x}{h}\] Extracting out the common term \(e^x\) in the numerator, we arrive at \[{\cal D}\, e^x = \frac{e^x e^h - e^x}{h} = e^x \left[\frac{e^h - 1}{h}\right] \tag{1}\]
Math expression eq-Dexp-27 shows that the slope function of \(e^x\) is proportional to \(e^x\) itself. The constant of proportionality is \([e^h - 1]/h\).
question id: lamb-put-bulb-1
Keep in mind that our interest is in \({\cal D}\, e^x\). The role of \(h\) is analogous to scaffolding on a building site. We seek to remove it eventually. Ideally, if we make \(h\) small enough, it will vanish.
question id: lamb-put-bulb-2
The instantaneous rate of change involves making \(h\) very small, but not quite zero. If you make \(h\) exactly zero, the result will be ambiguous.
Inf
NaN
Bogus
An error message result.
question id: lamb-put-bulb-3
1
It varies with \(h\), but is always around 1.5 for \(h\) small enough
There is no particular result.
question id: lamb-put-bulb-4
`
Exercise 4 Suppose you want to borrow \(L=\) $10,000 to buy a car. You will have to make equal monthly payments \(P\) for five years (60 months) to pay off the loan. At an interest rate of \(r\) per month, your payment will be: \[P(L, r) = L r \frac{(1+r)^{60}}{(1+r)^{60} - 1}\]
For example, if the interest rate is 0.005 per month (roughly 6% = 0.06 per year), your monthly payment would be
You are a good negotiator and are trying to talk the car dealer down to a zero-percent loan. The dealer plugs \(r=0\) into the car-loan calculator. His response to you: “Sorry, but our system cannot handle zero-percent loans.”
You press the point. “Excuse me, Sir, but I’m reading MOSAIC Calculus. That NaN
appearing on your screen is the result of a division by zero in the formula. But there is no reason to do that division. At zero percent interest, the monthly payment will simply be the amount borrowed divided by the number of months of the loan, so $10,000 divided by 60, giving a monthly payment of $166.6666… .”
“I’m sure you’re very good at math,” the dealer says, “and I’m willing to agree to a monthly payment of $166.67, but I cannot process any loan without going through this computer system. And zero percent won’t fly.”
Having studied limits, you have an idea. “Suppose we agree on a non-zero interest rate—which your loan system can handle—that produces a monthly payment of $166.67? Can we write up the loan that way?”
Dealer: “I’d be happy to do that, but obviously it is impossible to find an interest rate greater than zero that gives the same result as your calculation for zero interest.”
You: “Let’s try. Instead of 0.005 for the monthly interest rate, put in \(r = 0.0005\). You’re computer can work with that, right?”
Dealer: “OK, let’s see. … Yes, the payment amount is $169.22. That’ll work.”
You: “We are making progress. But my sense of mathematical honor insists that we find an interest rate that gives $166.67, as close as we can practically get to the exact answer result $166.6666… .”
question id: low-low-rates-3
Later in the day, you tell the story to your roommate, who is a computer science major. She says that you were lucky. “Computer arithmetic isn’t the same as mathematical arithmetic. Computer calculations with very small numbers sometimes give different results than you would expect mathematically. I bet if you tried an even smaller \(r\), you would have gotten different answers or even just nonsense.”
The result steadily gets closer to 166.6666….
At first, the result gets closer to 166.6666, but then as \(r\) gets smaller the result stays the same.
The result falls below 166.666 and stays there.
The computer output shows a result of infinity.
question id: low-low-rates-4
See this news story from 1991 for a real-world example of a similar flaw.
Exercise 11 Use Active R chunk lst-rhinosaurus-toss-gloves1 for evaluating \(\sin(x)/x\) at tiny values of \(x\) or for plotting it on a domain that is somewhat wider. To answer the questions, you just have to change tinyx
.
tinyx
to make it even smaller. You can stop when you get tired. Does \(\sin(x)/x\) evaluate to something sensible for such tiny input values? If so, what value?question id: rhinosaurus-toss-gloves-1
Saying, “so small that I got tired typing the zeros” is not a convincing definition of “small” to a mathematician. For example, 0.0000000000000001 parsec (a unit of length) seems small but it is equivalent to about 10 feet—not so small. Mathematicians want you to take “small” to the limit, an arbitrarily large number of zeros, and when you’re done with that, even more zeros.
Fortunately, R and other computer languages have a scientific notation that allows you just to name the number of zeros you want after the decimal point. For instance 1e-2
is \(0.01\)—one zero. Similarly 1e-20
is \(0.00000000000000000001\), nineteen zeros.
x = 1e-31
, calculate sin(x)/x
. Then double the number of zeros, keep on doubling the number of zeros. The result will continue to be 1 … until it eventually becomes NaN
. How many zeros are there in the x
that produces NaN
as the answer to sin(x)/x
?question id: rhinosaurus-toss-gloves-2
What’s happening here has more to do with the nature of computers than the nature of numbers. Computers (in the manner they are ordinarily programmed) use packets of bits to represent numbers, and the chips have been engineered to make those bit packets respond to arithmetic operations as if they were the numbers they represent. A typical computer number, like 0.001, uses 64 bits in a special, standard format. Since there is a finite number of bits, there is a largest possible non-Inf
number and a smallest possible non-zero number. According to the IEEE standard for “floating-point” arithmetic the largest non-Inf
number is around 1e300
and the smallest non-zero number is around 1e-320
. This failure to behave like genuine mathematical numbers is called “overflow” (for large numbers which turn into Inf
) and “underflow” (for small numbers which turn into 0
).
1e300
, 1e301
and so on until you find the smallest 1e???
that prints as Inf
. Similarly, try numbers in the format 1e-320
and 1e-321
until you find the largest one that prints out as exactly zero. What are those two numbers?
1e305
and 1e-322
1e306
and 1e-323
1e308
and 1e-324
1e309
and 1e-327
question id: rhinosaurus-toss-gloves-3
The polynomial computer does not have any problem with overflow or underflow. The key to success is to write the Taylor polynomial for functions such as \(\sin(x)\) or \(x\) or \(x^2\) near \(x_0 = 0\). Such polynomials will always look like:
\[f(x) = a_1 x^1 + a_2 x^2 + a_3 x^3 + \cdots\]
What’s special here is that the a_0
term does not need to be included in the polynomial, since \(f(0) = 0\).
sin()
tan()
atan()
acos()
question id: rhinosaurus-toss-gloves-4
These zero divided by zero problems (like \(\sin(x) / x\)) always involve a ratio of two functions (\(\sin(x)\) and \(x\) here) that don’t need the \(a_0\) term in their Taylor series around \(x_0 = 0\). That makes them just a little bit simpler.
What’s more important than simpler is that, for the expansion of such functions to study the limit at \(x \rightarrow 0\), we only need the first term with a non-zero coefficient \(a_k\) to represent the function with complete accuracy.
Why? Consider the 2nd-order Taylor polynomial \(a_1 x + a_2 x^2\). If we are to be able to safely disregard the \(a_2\) term it is because that term, for small \(x\) is much, much smaller than the \(a_1 x\) term. And we can always choose non-zero \(x\) to make this so.
For instance, suppose our polynomial were \(x + 100 x^2\). For \(x=0.1\), the first and second terms are the same size; we need them both for accuracy. For \(x=0.01\), the second term is 1/100 the size of the first term, maybe we don’t need the second term so much. You can always make \(x\) so small that anyone will be satisfied that the second term is utterly negligible compared to the first.
Here’s the method:
Suppose you have a function \(f(x) \equiv u(x)/v(x)\) where \[\lim_{x\rightarrow 0} u(x) = 0\ \ \ \text{and}\ \ \ \lim_{x\rightarrow 0} v(x) = 0\] Given this, \(f(0)\) is not defined. But we can ask whether there is a sensible value that can be plugged in in place of \(f(0)\) that will cause the modified \(f()\) to be continuous at \(x=0\).
Step 1: Write the Taylor polynomial expansion around \(x-0 = 0\) for both \(u(x)\) and \(v(x)\). If both expansions have a non-zero first coefficient, you can stop there. Now we have:
\[u(x) \approx a_1 x\] \[v(x) \approx b_1 x\] where \(a_1 = \partial_x u(0)\) and \(b_1 = \partial_x v(0)\).
Step 2: Divide the polynomial (really just linear!) expansion of \(u()\) by the expansion of \(v()\) to get
\[\lim_{x\rightarrow 0}\frac{u(x)}{v(x)} = \lim_{x\rightarrow 0} \frac{a_1 x}{b_1 x} = \frac{a_1}{b_1}\]
that is the answer, \(a_1/b_1\), at least when \(b_1 \neq 0\). We will come back to that case later.
\(a_1 = 1\) and \(b_1 = 1\)
\(a_1 = \pi\) and \(b_1 = \pi\)
\(a_1 = -1\) and \(b_1 = -1\)
\(a_1 = 0\) and \(b_1 = -1\)
question id: rhinosaurus-toss-gloves-5
Sometimes the singularity is at some non-zero \(x\). For instance, \[h(x) \equiv \frac{x^2 - 16}{x - 4}\] The divide-by-zero comes into play when \(x=4\). So is there a sensible value to plug in for \(h(4)\) to replace the singularity.
Here, write your Taylor polynomials around \(x_0 = 4\), the location of the singularity. We will get:
\[x^2 - 16 = a_1 (x-4) + a_2 (x-4)^2 + \cdots\] \[x - 4 = b_1 (x-4)\] Using Taylor’s formula for coefficients we will get \[a_1 = \partial_x (x^2 - 16)\left.\right|_{x=4} = 2x\left.\right|_{x=4} = 8\] \[b_1 = \partial_x (x - 4) = 1\]
Consequently, \(\lim_{x\rightarrow 4} \frac{x^2 - 16}{x - 4} = 8\)
We’ve been discussing ratios of functions where the ratio cannot be calculated at the singularity using simply the limits of the functions approaching that singularity. (For instance \(\lim_{x\rightarrow 0} \sin(x) = 0\) and \(\lim_{x\rightarrow 0} x = 0\), but knowing this does not tell us what \(\lim_{x\rightarrow 0} \frac{\sin(x)}{x}\) will be. These are called “indeterminate forms.” As you’ve seen, if we know more about the functions than merely their individual limits, we can sometimes resolve the indeterminacy. Here we are doing that by writing each function as a low-order polynomial.
The indeterminate form \(\lim_{x\rightarrow 0} \frac{\sin(x)}{x}\) might be said to have the “shape” 0/0. But 0/0 is just a notation about the limits of the two individual functions. There are indeterminate forms with other shapes: \[\frac{0}{0}\ \ \ \ \ \ \frac{\infty}{\infty}\ \ \ \ \ \ 0\cdot {\infty\ \ \ \ \ 0^0\ \ \ \ \ \ \infty^0\ \ \ \ \ \ 1^\infty\ \ \ \ \ \ \infty - \infty}\] Keep in mind that something like \(0 \infty\) is not a multiplication problem but rather shorthand for \(u(x) v(x)\) where \(\lim_{x\rightarrow x_0} u(x) = 0\) and \(\lim_{x\rightarrow x_0} v(x) \rightarrow \infty\).
There is a variety of algebraic tricks to try to transform these different shapes of indeterminate forms into a ratio of functions, each of which goes to zero at the relevant \(x_0\). Once that is done, you can apply the method described above.
Indeterminate forms have been a bonanza for the author of calculus book exercises, who can write an large number of examples, many of which have never been seen in the wild.
One shortcut that works in practice is to make a graph of the indeterminate form near the singularity. If the limit as \(x\) approaches the singularity is a finite number, you can read off the result from the graph.
In Active R chunk lst-x-logx, the function \(g(x)\equiv x \ln(x)\) is set up. This function has a singularity at \(x=0\). Examine the plot and determine where the function value is going as you walk along the graph to the singularity.
Exercise 5 It is simple enough to say to an elementary-school pupil that “division by zero is an error.” But what is a computer to do when, for whatever reason, it divides by an amount that is zero generated by a not-erroneous early computation? In many computer languages (including R) the programmer can arrange that a unacceptable situation causes an error message to be “thrown.” This generally causes the program to stop. See this account of the USS Yorktown (under “Smart ship testbed”) for an example of why stopping the program is not always an acceptable practice.
Modern computer arithmetic systems, implemented in hardware, take a different approach to division by zero. Rather than an error, divide-by-zero calculations return a definite value called “not a number” and printed as NaN
. Rather than stopping the program, calculations can continue on. Any calculation involving NaN
, such as 3 + NaN
or NaN/NaN
produces a NaN
. This clever design allows the system-level programmer to determine where and whether a division by zero needs to be dealt with.
NaN
as the output.1/0
0/0
0*NaN
NaN/NaN
5e20 + NaN*1e-50
NaN - NaN
NaN * Inf
sin(NaN)
sqrt(NaN)
Inf/Inf
1/Inf
question id: NaN-numerics-0
The NaN
strategy is effective in plotting. Programs such as slice_plot()
simply ignore points whose value is NaN
. As a result, the shape of the graph is determined by those points nearby, not at, the singularity. The overall effect is that the graph indicates whether the singularity is removable without the user having to do any juggling with NaN
.
For each of the following functions of \(x\), use the following procedure:
NaN
. If so, the function has a singularity at \(x^\star\).Example: \(f(x) \equiv \frac{x^2 + 2x}{x^3 + 4x}\) at \(x^\star = 0\). The sandbox code has been set up for steps (1) and (2). You will have to do step (3) by hand, using algebra. Notice that f(xstar)
returns NaN
, indicating that there is a singularity at xstar
. The graph, however, is perfectly smooth and has a value of 0.5 at xstar
. So the limit \(\lim_{x\rightarrow 0} f(x) = 0.5\). Using algebra, we can rewrite \(f(x)\)
\[f(x) = \frac{x^2 + 2x}{x^3 + 4x} = \frac{x}{x} \frac{x + 2}{x^2 + 4}.\] Eliminating the \(\frac{x}{x}\), the remainder of the expression evaluates to 1/2 at \(x=0\) without any problematic division by zero. This agrees perfectly with the graph made by Active R chunk lst-NaN-numerics-1.
question id: NaN-numerics-1
question id: NaN-numerics-2
question id: NaN-numerics-3
question id: NaN-numerics-4
question id: NaN-numerics-5
Exercise 6 Consider the function \(f(x) \equiv \frac{x-3}{x-3}\). Simple algebra suggests that the output will be 1 for any \(x\) value input, but the rules are division mean that the output is undefined for input \(x=3\). Thus, the function has a “singularity” at \(x=3\). We want to test whether this singularity is removable, that is, whether the function has a limit as \(x \rightarrow 3\). The computer can help us out. The basic idea is to try inputs closer and closer to \(x=3\) and see if the output behaves nicely. To help you distinguish between non-removable and removable singularities (or, in the words of the song “Santa Claus is Coming to Town”, to know “who’s naughty or nice”), Active R chunk lst-singularity-numerics defines the function \(f()\) and a similar function \(g()\). Try out both \(f()\) and \(g()\) for a sequence of inputs that get closer and closer to \(x=3\) but are not exactly 3.
Essay 1: Explain what the options(digits=20)
statement is doing. You can do some experiments by changing the 20 to something much smaller.
Which of these statements is true about the successive values in x_sequence
?
They approach 3 from below.
They approach 3 from above.
They oscillate, getting nearer and farther from 3.
question id: singularity-numerics-1
Which of these statements is true about the successive values in f(x_sequence)
?
They approach 1 from below.
They approach 1 from above.
Their value is constant at 1.
question id: YFpUUK
Essay 2: Explain what about the output of the function \(g()\) qualifies it for being numerically ‘naughty’ near \(x=3\).
Exercise 7 For each exercise, you are given a series of intervals that get smaller and smaller. Your job is to calculate the average rate of change of the function \(f(x) \equiv x^2\) for each of the intervals. As the width of the intervals approach zero, our average rates of change become better approximations of the instantaneous rate of change. You should use the results you calculate to make an informed estimate of the instantaneous rate of change.
A. Use these three intervals to estimate the instantaneous rate of change \(\partial_x f(x=3)\) - [3, 3.1] - [3, 3.01] - [3, 3.001]
B. Use these three intervals to estimate the instantaneous rate of change \(\partial_x f(x=5)\) - [4.9, 5] - [4.99, 5] - [4.999, 5]
C. Use these three intervals to estimate the instantaneous rate of change \(\partial_x f(x=-2)\) - [-2, -1.9] - [-2, -1.99] - [-2, -1.999]
Exercise 8 Using numerical finite-difference differentiation, confirm the rules given in sec-d-pattern-book for the pattern-book functions.
Find a small-enough \(h\) so that the error over the domain 0 to 1 is less than \(10^{-6}\).
Find an even smaller \(h\) that produces a substantially larger error than in (a).
Exercise 9
Which of the line segments is tangent to the curve at the point marked with a dot?
question id: goat-pay-pot
Exercise 10 Let’s explore the derivative \(\partial_x \frac{1}{x}\) using the evanescent-h method.
The general definition of the derivative of a function \(f(x)\) is the limit as \(h\rightarrow 0\) of a ratio: \[\partial_x f(x) \equiv \underbrace{\lim_{h \rightarrow 0}}_\text{limit}\ \underbrace{\frac{f(x + h) - f(h)}{h}}_\text{ratio}\ .\] For an evanescent-h calculation, we put aside for a moment \(\lim_{h\rightarrow0}\) and operate algebraically on the ratio. In our operations, we treat \(h\) as non-zero.
Substitute in \(f(x) = 1/x\) and show step-by-step that the ratio (for non-zero \(h\)) is equivalent to \(- \frac{1}{x^2 + hx}\).
The domain of the function \(1/x\) is the entire number line, except zero. Examine the expression \(- \frac{1}{x^2 + hx}\) to see if setting \(h\) to zero would ever involve a division by zero for any \(x\) in the domain of \(1/x\).
If setting \(h=0\) cannot cause a division by zero, it is legitimate to do so. Set \(h=0\) and simplify the expression \(- \frac{1}{x^2 + hx}\). Compare this to the rule for \(\partial_x \frac{1}{x}\) given in sec-d-pattern-book.