Avoiding Numeric Instability

I. Cnop
Vrije Universiteit Brussel
icnop@vub.ac.be

A classic limit

Plot[x Log[x],  {x, 0, 3}, AspectRatio -> Automatic]

Plot[x Log[x],  {x, 0, 0.1}, 
      AspectRatio -> Automatic]

What is the slope near the origin?

(∂_x) (x Log[x])

Limit[%, x -> 0]

Did we see this by magnifying the last graph?

Plot[x Log[x],  {x, 0, 0.001}, 
      AspectRatio -> Automatic]

Well to me it looks more like -10, not really making contact with the vertical axis.

The number   e   as a limit

Recall the classic def of  e .

A plot

Plot[(((x + 1))^(1/x)), {x, 0, 4}]

Plot[(((x + 1))^(1/x)), {x, 0, 0.001}]

Very close to 0  we get numeric instability

Because no CAS can compute millionth powers of reals accurately. 
If in the future numeric packages can do such computations, the same problem will still arise by the hundred millionth or many billionth powers.

Plot[(((x + 1))^(1/x)), {x, 0, 
        0.0000001}]

But Mathematica  knows the value of this limit:

Some limits

(lim_(x -> ∞)) (((1 + (1/x)))^x)

and we can even obtain information about how this limit is obtained:

((∂(((1 + (1/x)))^x))/(∂x))

(lim_(x -> ∞)) %

It is a nice asymptote.
Inverting near zero we get

(lim_(x -> 0)) (((x + 1))^(1/x))

Direct computation of the slope near zero  is not easy:

((∂(((x + 1))^(1/x)))/(∂x))

Limit[%, x -> 0]

Here is why:

((∂log((((x + 1))^(1/x))))/(∂x))

We recognise factors and do the factoring

Apart[%]

Limit[%, x -> 0]

Therefore the slope is

((-1)/(2 E))

The problem is that we cannot get this information in a plot near zero.

A simple reasoning

Taking a long polynomial

Use Mc Laurin development to approximate the previous plots by some polynomials. 
We do not have to go all the way to

Series[log(x + 1), {x, 0, 50}];

exponent = Normal[%]

Series[(E^(exponent/x)), {x, 0, 
        50}];

or we do not have to go even to 20

approxim = Normal[Series[%, {x, 0, 20}]]

which is unreadable and of no intrest! 
In fact, a quadratic approximation is sufficient for our purposes:

A smaller degree  polynomial is sufficient

It is sufficient to write an approxim of second order,which we could do even by hand :

Series[log(x + 1), {x, 0, 2}];

exponent2 = Normal[%]

Series[(E^(exponent2/x)), {x, 0, 
        1}]

tangent  = Normal[%]

Plot[tangent, {x, 0, 1}, PlotRange -> {2, 3}, 
      AspectRatio -> Automatic]

\(Show[
	
      Plot[{tangent, approxim}, {x, 0,  .8}, AspectRatio -> Automatic]
, 
      PlotRange -> {{0, 0.2}, {2.6, 2.8}}]\)

The drawing of the tangent avoids the numeric instability in the graph near 0:

Show[{Plot[(((x + 1))^(1/x)),{x, 0.000001, 0.00001},
PlotRange->(({{0, 0.00001`}, {E - 
                              0.00002`, E + 2.`*^-6}}))],
Plot[tangent,{x, 0, 0.00001},
PlotRange->(({{0, 0.00001`}, {E - 
                              0.00002`, E + 2.`*^-6}}))]}
,AspectRatio -> Automatic]

Comment

This notebook shows how simple mathematical reasoning can avoid serious pitfalls and inaccuracies in ICT.

Extensions

Compute the following limits and try to check these by plotting.

(lim_(x -> 0)) ((((x^3) + 3 (x^2) + 1))^((1/(x^2))))

Plot[((((x^3) + 3 (x^2) + 1))^((1/(x^2)))), {x, 
        0,  .1}]

(lim_(x -> 0)) ((((x^3) + 1))^((1/(x^2))))

Plot[((((x^3) + 1))^((1/(x^2)))), {x, 
        0,  .1}]

The following shows bad behavior:

(lim_(x -> 0)) (((x + 1))^((1/(x^2))))

Plot[(((x + 1))^((1/(x^2)))), {x, 
        0.01,  .1}]

These are all situations of the form

((f(x))^(g(x)))

If   f(x)   goes to   1   as   x^m   and the exponent   g(x)    goes to   0   as   x^-m  , then we have the same conclusion, the limit is finite
and we can control the way the expression goes to the finite lim it.
Different behavior make the expressions  tend to 1 or  infintiy.