MATLAB Examples

# How does limest work?

How might one best compute the limit of a function at a specific point using numerical techniques? Since we need to compute a limit, the function will generally be singular in some fashion at that point. Either a discontinuity exists, or if we tried to evaluate the function, we might encounter a 0/0 problem. In either case, we need to infer the limit purely from function evaluations taken in the vicinity of the limit point. In some cases, we might even find that if we try to evaluate the function too closely to the point in question, there are numerical problems that occur.

Likewise, if we sample our function too far away from the point, the function values will not be a good indicator of the limiting value.

So an adaptive numerical tool must be careful, like Goldilocks in the fairy tale, it must find the middle ground.

The approach taken by limest is to sample the function at a geometric sequence of points approaching the limit point. (See also my derivest tools and residueEst.)

For example, suppose we wish to compute the limit of the function f(x) = sin(x)./x, at x = 0. (The limit is 1.) We cannot evaluate this function at 0, since MATLAB will return a NaN there. If we evaluate the function too far away from 0, then it will not be an accurate indicator of the true limit.

```fun = @(x) sin(x)./x; ezplot(fun) format long g fun(0) fun(.000001) fun(.001) ```
```ans = NaN ans = 0.999999999999833 ans = 0.999999833333342 ```

A somewhat worse test case is the function f(x) = exp(x) - 1 - x. Ezplot suggests the limit is 0.5 at x = 0.

```fun = @(x) (exp(x) - 1 - x)./x.^2; ezplot(fun,[-1,1]) grid on ```

At x = 0, we get NaN.

```fun(0) ```
```ans = NaN ```

If we evaluate the function near the limit point, there are serious numerical problems. Below approximately sqrt(eps), this function returns pure garbage.

```x = 10.^(-1:-1:-15)'; [x,fun(x)] ```
```ans = 0.1 0.517091807564771 0.01 0.501670841679489 0.001 0.500166708384577 0.0001 0.500016671409274 1e-05 0.500000696499348 1e-06 0.499962183698476 1e-07 0.494336803055876 1e-08 -0.607747099184471 1e-09 82.7403709368088 1e-10 827.403709626582 1e-11 8274.03708980347 1e-12 88900582.3410317 1e-13 -7992778373.59144 1e-14 -79927783735.9112 1e-15 110223024625156 ```

Limest solves this problem by choosing some offset, dx, from x0. Pick some small offset from the limit point, say 1e-8. Now, evaluate the function at a sequence of points:

`x0 + dx, x0 + k*dx, x0 + k^2*dx, x0 + k^3*dx, ...`

The default in limest is to use the geometric factor of k = 4.

```x0 = 0; dx = 1e-8; k = 4; delta = dx*k.^[-10:15]' ```
```delta = 9.5367431640625e-15 3.814697265625e-14 1.52587890625e-13 6.103515625e-13 2.44140625e-12 9.765625e-12 3.90625e-11 1.5625e-10 6.25e-10 2.5e-09 1e-08 4e-08 1.6e-07 6.4e-07 2.56e-06 1.024e-05 4.096e-05 0.00016384 0.00065536 0.00262144 0.01048576 0.04194304 0.16777216 0.67108864 2.68435456 10.73741824 ```

Now, evaluate the function at each of these points, x0 + delta

```fun = @(x) (exp(x) - 1 - x)./x.^2; fun_of_x = fun(x0+delta); [delta,fun_of_x] ```
```ans = 9.5367431640625e-15 122868749999.998 3.814697265625e-14 30717187499.9994 1.52587890625e-13 -1857446289.06264 6.103515625e-13 131684875.488247 2.44140625e-12 -4331684.1125574 9.765625e-12 -1082921.02813935 3.90625e-11 20308.0475324997 1.5625e-10 -4017.93513460437 6.25e-10 132.384593565068 2.5e-09 -2.43098839673788 1e-08 -0.607747099184471 4e-08 0.541952615594605 1.6e-07 0.499780083853781 6.4e-07 0.50007897264343 2.56e-06 0.499984288251391 1.024e-05 0.500000851218376 4.096e-05 0.500006844844694 0.00016384 0.50002731050272 0.00065536 0.500109244595775 0.00262144 0.500437193156052 0.01048576 0.501752217589173 0.04194304 0.507064426657616 0.16777216 0.529175319451948 0.67108864 0.63344285370923 2.68435456 1.52161462365596 10.73741824 399.292392213823 ```

Use a polynomial to fit a subsequence of points. Polyfit would do in a pinch, although more efficient methods are employed in practice. Note that if we look at the points that are very close to z0, then the polynomial may have strange coefficients.

To compute the limit, we are really only interested in the value of the constant term for this polynomial model. Essentially, this is the extrapolated prediction of the limiting value of our product extrapolated down to delta = 0. In effect, this process is highly related to the idea of Richardson extrapolation , except that we do not use a polynomial interpolant to derive the necessary coefficients. The use of a regression polynomial provides the error estimates from limest.

```P = polyfit(delta(1:4),fun_of_x(1:4),2) ```
```P = 1.26830908463711e+36 -9.41559254138715e+23 102784215715.195 ```

Our estimate of the limit at x = 0 is just P1(end), the constant term in the polynomial. Recall that the true limit for ths function was 0.5. This estimate is a poor one.

```P(end) ```
```ans = 102784215715.195 ```

A nice feature of this sequence of points, is we can re-use the function values for the points, adding one more to the end of the sequence, and dropping the first one in our call to polyfit. (By the way, limest carefully scales its problems to avoid the problems that polyfit would have seen. The actual code does not call polyfit anyway. For this example, I'll just turn off those warning messages.)

```warning('off','MATLAB:polyfit:RepeatedPointsOrRescale') P = zeros(23,3); for i = 2:23 P(i,:) = polyfit(delta(i:i+3),fun_of_x(i:i+3),2); end P ```
```P = 0 0 0 1.73964236149361e+34 -5.09621201795639e+22 20907400210.698 -6.57631553730906e+31 7.69529260258972e+20 -1255568987.79351 3.00505692815173e+29 -1.40675165291052e+19 91753708.1986385 -6.71190236062665e+26 1.27875217994924e+17 -3614068.6586673 -9.4990833950136e+24 7.14220588600797e+15 -759299.095081526 1.04188714046047e+22 -30971161588071.8 12448.979532125 -1.37785872691531e+20 1654680714319.14 -2791.97050517086 2.88122013179984e+17 -13819800588.0643 93.085996068978 -456180813797903 89403838.2743891 -2.13468277075317 -9264964688531.65 7118379.81187322 -0.267173798691355 22100106854.8417 -68096.1960063015 0.529708609951636 -5517277.86172314 67.5795567367358 0.499885550415498 155267.272474548 -7.36018231272769 0.500048373377069 -2126.20836704672 0.587163326012658 0.499988101779456 -6.8596909037101 0.172097785181083 0.49999939995763 0.0512616218594046 0.166636188126643 0.500000014200183 0.04186392442111 0.166665324579313 0.500000002022359 0.0421156754180716 0.166662478241848 0.500000004949546 0.0435030810610506 0.166597904348441 0.500000323876313 0.0497270418428591 0.165439254718907 0.500023205364048 0.0913430121766019 0.134458010136504 0.502467086279486 4.68115184308082 -13.5267488831478 4.8013380041388 ```

See how, as we move along this sequence using a sliding window of 4 points, the constant terms will approach 0.5. Then at some point, we move just too far away from the pole, and our extrapolated estimate of the limit becomes poor again.

```[delta(1:23),P(:,end)] ```
```ans = 9.5367431640625e-15 0 3.814697265625e-14 20907400210.698 1.52587890625e-13 -1255568987.79351 6.103515625e-13 91753708.1986385 2.44140625e-12 -3614068.6586673 9.765625e-12 -759299.095081526 3.90625e-11 12448.979532125 1.5625e-10 -2791.97050517086 6.25e-10 93.085996068978 2.5e-09 -2.13468277075317 1e-08 -0.267173798691355 4e-08 0.529708609951636 1.6e-07 0.499885550415498 6.4e-07 0.500048373377069 2.56e-06 0.499988101779456 1.024e-05 0.49999939995763 4.096e-05 0.500000014200183 0.00016384 0.500000002022359 0.00065536 0.500000004949546 0.00262144 0.500000323876313 0.01048576 0.500023205364048 0.04194304 0.502467086279486 0.16777216 4.8013380041388 ```

The trick is to learn from Goldilocks. Choose a prediction of the limit for some model that is not too close to the limit point, nor too far away. The choice is made by a simple set of rules. First, discard any predictions of the limit that are either NaN or inf. Then trim off a few more predictions, leaving only those predictions in the middle. Next, each prediction is made from a polynomial model with ONE more data point than coefficients to estimate. This yields an estimate of the uncertainty in the constant term using standard statistical methodologies. While a 95% confidence limit has no true statistical meaning, since the data is not truly random, we can still use the information.

Limest chooses the model that had the narrowest confidence bounds around the constant term.

```[res,errest] = limest(fun,0) ```
```res = 0.499999999681485 errest = 2.20308196660258e-09 ```