MATLAB Examples

## How does residueEst work?

```%{ How might one best compute the residue of a complex function at a known pole of low order using numerical techniques? Since the function is not finite at the pole, we cannot evaluate it directly there. The approach taken by residueEst is to sample the function at a geometric sequence of points approaching the pole. This is not unlike the approach I took with my derivest tools. For example, suppose we wish to compute the residue of the first order pole of 1./sin(z), at z = pi. (The residue is -1.) For a first order pole, the residue is given by the limit: lim (f(z)*(z-z0)) z --> z0 ResidueEst solves this problem by choosing some offset, dz, from z0. Pick some small offset, say 1e-16. Now, evaluate the function at a sequence of points z0 + dz, z0 + k*dz, z0 + k^2*dz, ..., z0 + k^n*dz. The default in derivest is to use the geometric factor of k = 4. %} format long g z0 = pi; dz = 1e-16; k = 4; delta = dz*k.^(0:30) % Now, evaluate the function at each of these points, z0 + delta, and % multiply by (z0 + delta - z0) == delta. fun = @(z) 1./sin(z); fun_of_z = fun(z0+delta).*delta ```
```delta = Columns 1 through 3 1e-16 4e-16 1.6e-15 Columns 4 through 6 6.4e-15 2.56e-14 1.024e-13 Columns 7 through 9 4.096e-13 1.6384e-12 6.5536e-12 Columns 10 through 12 2.62144e-11 1.048576e-10 4.194304e-10 Columns 13 through 15 1.6777216e-09 6.7108864e-09 2.68435456e-08 Columns 16 through 18 1.073741824e-07 4.294967296e-07 1.7179869184e-06 Columns 19 through 21 6.8719476736e-06 2.74877906944e-05 0.0001099511627776 Columns 22 through 24 0.0004398046511104 0.0017592186044416 0.0070368744177664 Columns 25 through 27 0.0281474976710656 0.112589990684262 0.45035996273705 Columns 28 through 30 1.8014398509482 7.20575940379279 28.8230376151712 Column 31 115.292150460685 fun_of_z = Columns 1 through 3 0.816561967659769 -1.24368623276475 -0.96741494953197 Columns 4 through 6 -1.05007818637943 -0.998645996304325 -0.999393504822813 Columns 7 through 9 -1.0006650249676 -1.00016932159893 -1.0000454724843 Columns 10 through 12 -0.999997574134254 -1.0000025407217 -1.00000060598878 Columns 13 through 15 -1.00000012230672 -1.00000000138628 -1.00000000424339 Columns 16 through 18 -1.00000000082177 -0.999999999966395 -1.0000000000115 Columns 19 through 21 -1.00000000003004 -1.00000000013473 -1.0000000020163 Columns 22 through 24 -1.00000003223861 -1.00000051580866 -1.00000825298128 Columns 25 through 27 -1.00013205914403 -1.00211587978403 -1.03462138368425 Columns 28 through 30 -1.85044046871399 -9.03930662202117 55.26083092194 Column 31 -142.066423215372 ```

## Use a polynomial to fit a subsequence of points. Polyfit would do.

Note that if we look at the points that are very close to z0, then the polynomial may have strange coefficients.

```% For a first order pole, we are really only interested in the value % of the constant terms for this polynomial model. Essentially, this % is the extrapolated prediction of the limiting value of our product % extrapolated down to delta == 0. P0 = polyfit(delta(1:4),fun_of_z(1:4),2) % A nice feature of this sequence of points, is we can re-use the % function values, building a sliding sequence of models. P0 = polyfit(delta(2:5),fun_of_z(2:5),2) P0 = polyfit(delta(3:6),fun_of_z(3:6),2) P0 = polyfit(delta(4:7),fun_of_z(4:7),2) P0 = polyfit(delta(5:8),fun_of_z(5:8),2) % See how, as we move along this sequence using a sliding window of 4 % points, the constant terms is approaching -1. P0 = polyfit(delta(5:8),fun_of_z(5:8),2) P0 = polyfit(delta(12:15),fun_of_z(12:15),2) % At some point, we move just too far away from the pole, and our % extrapolated estimate at the pole becomes poor again. P0 = polyfit(delta(26:29),fun_of_z(26:29),2) P0 = polyfit(delta(27:30),fun_of_z(27:30),2) ```
```Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT. P0 = 1.40722697171131e+29 -1.08866090694287e+15 0.165206417147826 Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT. P0 = -7.74680625721916e+26 25552618526783.7 -1.14677041257162 Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT. P0 = 4.67339267331006e+24 -451780266912.317 -1.00159018384202 Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT. P0 = -1.05363184798865e+24 514903808559.053 -1.03508462443342 Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT. P0 = 3.39997173195685e+21 -6525897035.24524 -0.9986021313655 Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT. P0 = 3.39997173195685e+21 -6525897035.24524 -0.9986021313655 Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT. P0 = -3107361376.59391 101.434422959505 -1.00000049047424 P0 = -0.114884717006806 -0.297892275588167 -0.928218832274888 P0 = 0.149132944160628 -2.4122832973759 0.881916750679264 ```

## The trick is learn from Goldilocks.

Choose a prediction of the pole for some model that is not too close to the pole, nor to 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.

ResidueEst choose the model that had the narrowest confidence bounds around the constant term.

```[res,errest] = residueEst(@(z) 1./sin(z),pi) % See that if the true residue is really is -1, then we should have % (approximately) % % abs(res - (-1)) < errest [abs(res - (-1)),errest] ```
```res = -1.00000000000222 errest = 1.23720470391972e-11 ans = 2.22333262911434e-12 1.23720470391972e-11 ```

## Residue of exp(z)/z, around z0 == 0

Here the residue should be exp(0) == 1

```[res,errest] = residueEst(@(z) exp(z)./z,0) ```
```res = 1 errest = 0 ```

## A second order pole, exp(3*z)./z.^2

The residue should be 3.

```[res,errest] = residueEst(@(z) exp(3*z)./z.^2,0,'poleorder',2) % The error estimate is a reaonable one here too. [abs(res - 3),errest] ```
```res = 3.00000000019214 errest = 2.26398426905353e-10 ans = 1.92144078425827e-10 2.26398426905353e-10 ```

## A second order pole, 1./z.^2, at z0 = 0

The residue should be 0

```[res,errest] = residueEst(@(z) 1./z.^2,0,'poleorder',2) ```
```res = 1.77121275308468e-15 errest = 1.83193682355051e-14 ```

## Another 2nd order pole with a zero residue, 1./sin(z).^2, z0 = pi

The residue here is zero

```[res,errest] = residueEst(@(z) 1./(sin(z).^2),pi,'poleorder',2) ```
```res = -6.43430530792682e-09 errest = 3.61963048877858e-08 ```