(This is a reprint of the second ever Cleve's Corner from the Winter 1990 MathWorks Newsletter).... read more >>

]]>(This is a reprint of the second ever Cleve's Corner from the Winter 1990 *MathWorks Newsletter*).

The other day at lunch with a couple of other MathWorks people, I posed the following problem:

“I'm thinking of two numbers. Their average is 3. What are the numbers?”

Now stop a minute and, before you read any further, provide your own answer to my question.

Of course, my problem doesn't have a “right” answer. The solution is not unique. The problem is “underdetermined” or “ill posed.”

Most of the people I was having lunch with wouldn't give me an answer. When I pressed them, one person finally said, “Both numbers could be 3." Another person said, “Yes, but one could be 6 and the other 0. That's what I hoped they would say. In some sense, these are both “nice” answers. Nobody said “23 and minus 17,” or “2.71828 and 3.28172.” Those would also have been been correct, but not as “nice.”

What does this have to do with MATLAB? Well, MATLAB doesn't balk at impossible problems. It will solve this one. And it doesn't without complaining like my buddies at lunch.

Naturally the problem is a matrix problem. It is a "system"

A*x = b

of one "simultaneous" linear equation in two unknowns. The matrix is

A = [1/2 1/2]

and the right-hand side is

b = 3

MATLAB’s backslash solves such equations. Typing

x = A\b

tell's me

x = 6 0

The solution is a 2-by-1 matrix representation of one of the “nice” answers I expected. The second half of the help entry for “\” gives some indication of where this solution came from.

If A is an m-by-n matrix with m < or > n and B is a column vector with m components, or a matrix with several such columns, then X = A\B is the solution in the least squares sense to the under- or overdetermined system of equations A*X = B. The effective rank, k, of A is determined from the QR decomposition with pivoting. A solution X is computed which has at most k nonzero components per column. If k < n this will usually not be the same solution as PINV(A)*B.

In our case, we have `m = 1` and `n = 2`. My matrix has full rank, but the most that can be is `k = 1`. So, we get a solution with at most one nonzero component. Even that isn't unique. There are exactly two solutions with one nonzero component, `[6 0]’` and `[0 6]’`. We get the one where the nonzero component has the smallest index.

What about my other “nice” solution? That comes from the pseudoinverse.

x = pinv(A)*b

produces

x = 3.0000 3.0000

(The fact that I got that I get 3.000 instead of the integer 3 without the trailing zeros indicates that there has been some roundoff error and I don't have the other “nice” solution exactly, but you can pursue that yourself if you really want to.)

Of all the possible solutions to `A*x = b`, this one, `x = [3 3]’`, is the “shortest.” It minimizes `norm(x)`. In fact `x = A\b` has

norm(x) = 6.0000

while `x = pinv(A)*b` has

norm(x) = 4.2426

Now, finally, we have uniqueness. Of all the possible solutions to `A*x = b`, the one which also minimizes `norm(x)` is unique.

So, MATLAB not only solves the problem, it gives us a choice between two different solutions, `x = A\b` and `x = pinv(A)*b`. I think the first solution, `A\b = [6 0]'`, is "nice", because it is "simple"; i.e, it has the fewest possible nonzero components. I think the second solution, `pinv(A)*b = [3 3]'`, is "nice" because it is "short" and "unique". None of the other possible solutions have such nice characteristics.

This problem, "Given the average of two numbers, find the numbers," captures the essence of many ill-posed and underdetermined problems. Computer tomography, which is the life-saving business of generating images from X-ray, magnetic resonance, and other scanners, is really a grown-up version of this question. Additional constraints, like minimum norm or fewest nonzero components, or "good looking picture," have to be specified to make it a reasonable mathematical and computational task.

By the way, I first learned about this "World's Simplest Impossible Problem" from Don Morrison, who also started the University of New Mexico's Computer Science Department, invented the Fast Fourier Transform before Cooley and Tukey, and, years ago, got me to move to New Mexico. Thanks for all those things, Don.

I have to confess that I wrote this anecdote about posing a problem at a MathWorks lunch before it really happened. The results of the actual experiment were even more interesting. As I had expected, everybody did grumble and complain about my problem. Then one person said "9 and -3." She had obviously picked up on the idea. Another person said "3 and 1." Luckily, he is not responsible for any numeric portions of MATLAB. But three people said "2 and 4." That is certainly another "nice" answer, but the constraints it satisfies are more subtle. They have something to do with requiring the solution to have integer components that are distinct, but near the average. It's harder to state and compute in MATLAB without just giving the answer.

Okay, now I'm thinking off three numbers whose average is $\pi$. What are those numbers?

Get
the MATLAB code

Published with MATLAB® R2018b

We will have a two-part minisymposium on "Bohemian Matrices" at ICIAM2019, the International Congress on Industrial and Applied Mathematics in Valencia, Spain, July 15-19. This is an outline of my talk.... read more >>

]]>We will have a two-part minisymposium on "Bohemian Matrices" at ICIAM2019, the International Congress on Industrial and Applied Mathematics in Valencia, Spain, July 15-19. This is an outline of my talk.

The colorful name "Bohemian Matrices" is the responsibility of Rob Corless and Steven Thornton of the University of Western Ontario. Quoting the web site.

When this project was in its early stages, our focus was on random integer matrices of bounded height. We originally used the phrase "bounded height integer matrix eigenvalues" to refer to our work. This led to the acronym BHIME which isn't too far off from "Bohemian".

The `gallery` is a collection of interesting matrices, maintained in MATLAB by Nick Higham. Many of them are Bohemian, or Bohemian wannabees. A catalog for the collection is available with

```
help gallery_
```

GALLERY Higham test matrices. [out1,out2,...] = GALLERY(matname, param1, param2, ...) takes matname, a string that is the name of a matrix family, and the family's input parameters. See the listing below for available matrix families. Most of the functions take an input argument that specifies the order of the matrix, and unless otherwise stated, return a single output. For additional information, type "help private/matname", where matname is the name of the matrix family. binomial Binomial matrix -- multiple of involutory matrix. cauchy Cauchy matrix. chebspec Chebyshev spectral differentiation matrix. chebvand Vandermonde-like matrix for the Chebyshev polynomials. chow Chow matrix -- a singular Toeplitz lower Hessenberg matrix. circul Circulant matrix. clement Clement matrix -- tridiagonal with zero diagonal entries. compar Comparison matrices. condex Counter-examples to matrix condition number estimators. cycol Matrix whose columns repeat cyclically. dorr Dorr matrix -- diagonally dominant, ill-conditioned, tridiagonal. (One or three output arguments, sparse) dramadah Matrix of ones and zeroes whose inverse has large integer entries. fiedler Fiedler matrix -- symmetric. forsythe Forsythe matrix -- a perturbed Jordan block. frank Frank matrix -- ill-conditioned eigenvalues. gcdmat GCD matrix. gearmat Gear matrix. grcar Grcar matrix -- a Toeplitz matrix with sensitive eigenvalues. hanowa Matrix whose eigenvalues lie on a vertical line in the complex plane. house Householder matrix. (Three output arguments) integerdata Array of arbitrary data from uniform distribution on specified range of integers invhess Inverse of an upper Hessenberg matrix. invol Involutory matrix. ipjfact Hankel matrix with factorial elements. (Two output arguments) jordbloc Jordan block matrix. kahan Kahan matrix -- upper trapezoidal. kms Kac-Murdock-Szego Toeplitz matrix. krylov Krylov matrix. lauchli Lauchli matrix -- rectangular. lehmer Lehmer matrix -- symmetric positive definite. leslie Leslie matrix. lesp Tridiagonal matrix with real, sensitive eigenvalues. lotkin Lotkin matrix. minij Symmetric positive definite matrix MIN(i,j). moler Moler matrix -- symmetric positive definite. neumann Singular matrix from the discrete Neumann problem (sparse). normaldata Array of arbitrary data from standard normal distribution orthog Orthogonal and nearly orthogonal matrices. parter Parter matrix -- a Toeplitz matrix with singular values near PI. pei Pei matrix. poisson Block tridiagonal matrix from Poisson's equation (sparse). prolate Prolate matrix -- symmetric, ill-conditioned Toeplitz matrix. qmult Pre-multiply matrix by random orthogonal matrix. randcolu Random matrix with normalized cols and specified singular values. randcorr Random correlation matrix with specified eigenvalues. randhess Random, orthogonal upper Hessenberg matrix. randjorth Random J-orthogonal (hyperbolic, pseudo-orthogonal) matrix. rando Random matrix with elements -1, 0 or 1. randsvd Random matrix with pre-assigned singular values and specified bandwidth. redheff Matrix of 0s and 1s of Redheffer. riemann Matrix associated with the Riemann hypothesis. ris Ris matrix -- a symmetric Hankel matrix. sampling Nonsymmetric matrix with integer, ill conditioned eigenvalues. smoke Smoke matrix -- complex, with a "smoke ring" pseudospectrum. toeppd Symmetric positive definite Toeplitz matrix. toeppen Pentadiagonal Toeplitz matrix (sparse). tridiag Tridiagonal matrix (sparse). triw Upper triangular matrix discussed by Wilkinson and others. uniformdata Array of arbitrary data from standard uniform distribution wathen Wathen matrix -- a finite element matrix (sparse, random entries). wilk Various specific matrices devised/discussed by Wilkinson. (Two output arguments) GALLERY(3) is a badly conditioned 3-by-3 matrix. GALLERY(5) is an interesting eigenvalue problem. Try to find its EXACT eigenvalues and eigenvectors. See also MAGIC, HILB, INVHILB, HADAMARD, PASCAL, ROSSER, VANDER, WILKINSON.

The matrices in the `See also` line were available before `gallery` was introduced. I think of them as part of the gallery.

My favorite matrix in the gallery is

G = gallery(5)

G = -9 11 -21 63 -252 70 -69 141 -421 1684 -575 575 -1149 3451 -13801 3891 -3891 7782 -23345 93365 1024 -1024 2048 -6144 24572

Let's compute its eigenvalues.

format long e eig(G)

ans = -3.472940132398842e-02 + 2.579009841174434e-02i -3.472940132398842e-02 - 2.579009841174434e-02i 1.377760760018629e-02 + 4.011025813393478e-02i 1.377760760018629e-02 - 4.011025813393478e-02i 4.190358744843689e-02 + 0.000000000000000e+00i

How accurate are these? What are the exact eigenvalues?

I will give readers of this post, and attendees at my talk, some time to contemplate these questions by delaying the answers until the end.

The `triw` matrix demonstrates that Gaussian elimination cannot reliably detect near singularity.

dbtype 1:2 private/triw

1 function T = triw(n, alpha, k, classname) 2 %TRIW Upper triangular matrix discussed by Wilkinson and others.

```
n = 12;
T = gallery('triw',n)
```

T = 1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 0 1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 0 0 1 -1 -1 -1 -1 -1 -1 -1 -1 -1 0 0 0 1 -1 -1 -1 -1 -1 -1 -1 -1 0 0 0 0 1 -1 -1 -1 -1 -1 -1 -1 0 0 0 0 0 1 -1 -1 -1 -1 -1 -1 0 0 0 0 0 0 1 -1 -1 -1 -1 -1 0 0 0 0 0 0 0 1 -1 -1 -1 -1 0 0 0 0 0 0 0 0 1 -1 -1 -1 0 0 0 0 0 0 0 0 0 1 -1 -1 0 0 0 0 0 0 0 0 0 0 1 -1 0 0 0 0 0 0 0 0 0 0 0 1

`T` is already upper triangular. It is the `U` of its own LU-decomposition. Since there are no small pivots on the diagonal, we might be tempted to conclude that `T` is well conditioned. However, let

```
x = 2.^(-(1:n-1))';
format rat
x(n) = x(n-1)
```

x = 1/2 1/4 1/8 1/16 1/32 1/64 1/128 1/256 1/512 1/1024 1/2048 1/2048

This vector is nearly a null vector for `T`.

Tx = T*x

Tx = 0 0 0 0 0 0 0 0 0 0 0 1/2048

The 1-norm condition number of `T` is at least as large as

norm(x,1)/norm(Tx,1)

ans = 2048

That is growing like `2^n`.

The `moler` matrix demonstrates that Cholesky decomposition of a symmetric, positive definite matrix cannot reliably detect near singularity either. The matrix can either be obtained from the gallery,

```
M = gallery('moler',n);
```

or generated from `T`,

```
format short
M = T'*T
```

M = 1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 2 0 0 0 0 0 0 0 0 0 0 -1 0 3 1 1 1 1 1 1 1 1 1 -1 0 1 4 2 2 2 2 2 2 2 2 -1 0 1 2 5 3 3 3 3 3 3 3 -1 0 1 2 3 6 4 4 4 4 4 4 -1 0 1 2 3 4 7 5 5 5 5 5 -1 0 1 2 3 4 5 8 6 6 6 6 -1 0 1 2 3 4 5 6 9 7 7 7 -1 0 1 2 3 4 5 6 7 10 8 8 -1 0 1 2 3 4 5 6 7 8 11 9 -1 0 1 2 3 4 5 6 7 8 9 12

I was surprised when Nick named the matrix after me. The comments in the code indicate he found the matrix in a 30-year old book by a friend of mine, John Nash.

dbtype 9:14 private/moler

9 % Nash [1] attributes the ALPHA = -1 matrix to Moler. 10 % 11 % Reference: 12 % [1] J. C. Nash, Compact Numerical Methods for Computers: Linear 13 % Algebra and Function Minimisation, second edition, Adam Hilger, 14 % Bristol, 1990 (Appendix 1).

John doesn't remember when he heard about the example from me and I don't remember where or where I first came across it.

We know that `T` should be the Cholesky decomposition of `M`, but let's make sure.

assert(isequal(chol(M),T))

We have already seen that, despite the fact that it has no small elements, `T` is badly conditioned. `M` must be at least as badly conditioned as `T`. In fact, it is far worse. A good null vector is provided by

[~,v] = condest(M);

Let's look at `v` alongside our `x`. They are nearly the same.

```
format long
[v, x]
```

ans = 0.500000178813913 0.500000000000000 0.250000089406957 0.250000000000000 0.125000223517391 0.125000000000000 0.062500469386522 0.062500000000000 0.031250949948913 0.031250000000000 0.015626905485761 0.015625000000000 0.007816313765488 0.007812500000000 0.003913878927961 0.003906250000000 0.001968383554413 0.001953125000000 0.001007079958072 0.000976562500000 0.000549316340766 0.000488281250000 0.000366210893844 0.000488281250000

The 1-norm condition of `M` is at least as large as

norm(v,1)/norm(M*v,1)

ans = 2.796202999840670e+06

It turns out that is growing faster than `4^n`.

But wait, you say, use column pivoting. Indeed, QR with column pivoting will detect the ill-condition of `T` or `M`. But the `kahan` matrix in the Gallery is a version of `T`, rescaled so that the columns all have nearly the same norm. Even column pivoting is fooled.

format short K = gallery('kahan',7)

K = 1.0000 -0.3624 -0.3624 -0.3624 -0.3624 -0.3624 -0.3624 0 0.9320 -0.3377 -0.3377 -0.3377 -0.3377 -0.3377 0 0 0.8687 -0.3148 -0.3148 -0.3148 -0.3148 0 0 0 0.8097 -0.2934 -0.2934 -0.2934 0 0 0 0 0.7546 -0.2734 -0.2734 0 0 0 0 0 0.7033 -0.2549 0 0 0 0 0 0 0.6555

I am very fond of the Rosser matrix.

R = rosser

R = 611 196 -192 407 -8 -52 -49 29 196 899 113 -192 -71 -43 -8 -44 -192 113 899 196 61 49 8 52 407 -192 196 611 8 44 59 -23 -8 -71 61 8 411 -599 208 208 -52 -43 49 44 -599 411 208 208 -49 -8 8 59 208 208 99 -911 29 -44 52 -23 208 208 -911 99

Before the success of the QR algorithm, the Rosser matrix was a challenge for matrix eigenvalue routines. Today we can use it to show off the Symbolic Math toolbox.

```
R = sym(R);
p = charpoly(R,'x')
f = factor(p)
e = solve(p)
```

p = x^8 - 4040*x^7 + 5080000*x^6 + 82518000*x^5 - 5327676250000*x^4 + 4287904631000000*x^3 - 1082852512000000000*x^2 + 106131000000000000*x f = [ x, x - 1020, x^2 - 1040500, x^2 - 1020*x + 100, x - 1000, x - 1000] e = 0 1000 1000 1020 510 - 100*26^(1/2) 100*26^(1/2) + 510 -10*10405^(1/2) 10*10405^(1/2)

Several years ago, I was totally surprised when I happened to compute the SVD of a nonsymmetric version of the Hilbert matrix.

```
format long
n = 20;
[I,J] = ndgrid(1:n);
P = 1./(I - J + 1/2);
svd(P)
```

ans = 3.141592653589794 3.141592653589794 3.141592653589794 3.141592653589793 3.141592653589793 3.141592653589792 3.141592653589791 3.141592653589776 3.141592653589108 3.141592653567518 3.141592652975738 3.141592639179606 3.141592364762737 3.141587705293030 3.141520331851230 3.140696048746875 3.132281782811693 3.063054039559126 2.646259806143500 1.170504602453951

Where did all those `pi` 's come from? Seymour Parter was able to explain the connection between this matrix and a classic theorem of Szego about the coefficients in the Fourier series of a square wave. The gallery has this matrix as `gallery('parter',n)`.

I must include the Pascal triangle.

```
help pascal_
```

PASCAL Pascal matrix. P = PASCAL(N) is the Pascal matrix of order N. P a symmetric positive definite matrix with integer entries, made up from Pascal's triangle. Its inverse has integer entries. PASCAL(N).^r is symmetric positive semidefinite for all nonnegative r. P = PASCAL(N,1) is the lower triangular Cholesky factor (up to signs of columns) of the Pascal matrix. P is involutory, that is, it is its own inverse. P = PASCAL(N,2) is a rotated version of PASCAL(N,1), with sign flipped if N is even. P is a cube root of the identity matrix. Reference: N. J. Higham, Accuracy and Stability of Numerical Algorithms, Second edition, Society for Industrial and Applied Mathematics, Philadelphia, 2002, Sec. 28.4.

P = pascal(12,2)

P = -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 11 10 9 8 7 6 5 4 3 2 1 0 -55 -45 -36 -28 -21 -15 -10 -6 -3 -1 0 0 165 120 84 56 35 20 10 4 1 0 0 0 -330 -210 -126 -70 -35 -15 -5 -1 0 0 0 0 462 252 126 56 21 6 1 0 0 0 0 0 -462 -210 -84 -28 -7 -1 0 0 0 0 0 0 330 120 36 8 1 0 0 0 0 0 0 0 -165 -45 -9 -1 0 0 0 0 0 0 0 0 55 10 1 0 0 0 0 0 0 0 0 0 -11 -1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0

With this sign pattern and arrangement of binomial coefficients, `P` is a cube root of the identity.

I = P*P*P

I = 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1

Joe Grcar will give us an interesting graphic.

dbtype 7:10 private/grcar

7 % References: 8 % [1] J. F. Grcar, Operator coefficient methods for linear equations, 9 % Report SAND89-8691, Sandia National Laboratories, Albuquerque, 10 % New Mexico, 1989 (Appendix 2).

```
n = 12;
Z = gallery('grcar',n)
```

Z = 1 1 1 1 0 0 0 0 0 0 0 0 -1 1 1 1 1 0 0 0 0 0 0 0 0 -1 1 1 1 1 0 0 0 0 0 0 0 0 -1 1 1 1 1 0 0 0 0 0 0 0 0 -1 1 1 1 1 0 0 0 0 0 0 0 0 -1 1 1 1 1 0 0 0 0 0 0 0 0 -1 1 1 1 1 0 0 0 0 0 0 0 0 -1 1 1 1 1 0 0 0 0 0 0 0 0 -1 1 1 1 1 0 0 0 0 0 0 0 0 -1 1 1 1 0 0 0 0 0 0 0 0 0 -1 1 1 0 0 0 0 0 0 0 0 0 0 -1 1

n = 100; Z = gallery('grcar',n); e = eig(Z); plot(e,'o')

Did you compute the exact eigenvalues of `gallery(5)` yet? Remember, the matrix is

G

G = -9 11 -21 63 -252 70 -69 141 -421 1684 -575 575 -1149 3451 -13801 3891 -3891 7782 -23345 93365 1024 -1024 2048 -6144 24572

MATLAB says the eigenvalues are

e = eig(G)

e = -0.034729401323988 + 0.025790098411744i -0.034729401323988 - 0.025790098411744i 0.013777607600186 + 0.040110258133935i 0.013777607600186 - 0.040110258133935i 0.041903587448437 + 0.000000000000000i

But, look at the powers of `G`. Since `G` has integer entries (it is Bohemian), its entries can be computed without roundoff. The fourth power is

G^4

ans = 0 0 0 0 0 -84 168 -420 1344 -5355 568 -1136 2840 -9088 36210 -3892 7784 -19460 62272 -248115 -1024 2048 -5120 16384 -65280

And the fifth power is

G^5

ans = 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

All zero.

The Cayley-Hamilton theorem tells us the characteristic polynomial must be

`x^5`

The eigenvalue, the only root of this polynomial, is zero with algebraic multiplicity 5.

MATLAB is effectively solving

`x^5` = roundoff

Hence the computed eigenvalues come from

`x` = (roundoff)^(1/5)

It's hard to recognize them as perturbations of zero.

Get
the MATLAB code

Published with MATLAB® R2018b

This is a summary of my talk at the conference Celebrating the Centenary of James H. Wilkinson's Birth at the University of Manchester, May 29.... read more >>

]]>This is a summary of my talk at the conference Celebrating the Centenary of James H. Wilkinson's Birth at the University of Manchester, May 29.

Jim Wilkinson knew that a matrix with badly conditioned simple eigenvalues must be close to a matrix with multiple eigenvalues. How can you find the nearest matrix where some of the eigenvalues have found their mates?

Much of this talk is my interpretation of the 2011 paper by Rafikul Alam, Shreemayee Bora, Ralph Byers and Michael Overton, "Characterization and construction of the nearest defective matrix via coalescence of pseudospectral components," in *Linear Algebra and its Applications*, vol. 435, pp. 494-513. Overton's accompanying code `neardefmat` is over 1000 lines of MATLAB® and is available from his web page.

In last week's blog post, "An Eigenvalue Sensitivity Example", I resurrected this example from the 1980 User's Guide. For a moment, let's forget how it was constructed and start with

A = [ -64 82 21; 144 -178 -46; -771 962 248]

A = -64 82 21 144 -178 -46 -771 962 248

The computed eigenvalues are

```
format long
lambda = eig(A)
```

lambda = 3.000000000003868 0.999999999998212 1.999999999997978

The exact eigenvalues are 1, 2 and 3. We've lost about four figures. This is predicted by the eigenvalue condition numbers,

```
format short
kappa = condeig(A)
```

kappa = 833.1092 450.7228 383.7564

These condition numbers come from the right eigenvector matrix

[X,~] = eig(A)

X = 0.0891 0.0735 -0.1089 -0.1782 -0.1923 0.1634 0.9800 0.9786 -0.9805

And the left eigenvector matrix

Y = inv(X)

Y = -521.9612 628.5984 162.7621 265.1820 -353.5760 -88.3940 -257.0058 275.3634 73.4302

The rows of `Y` are left eigenvectors. They are normalized so that their dot product with the right eigenvectors is one.

Y*X

ans = 1.0000 -0.0000 -0.0000 -0.0000 1.0000 -0.0000 -0.0000 0.0000 1.0000

A function introduced in 2017b that computes the 2-norm of the columns of a matrix comes in handy at this point.

kappa = (vecnorm(Y').*vecnorm(X))'

kappa = 833.1092 450.7228 383.7564

EigTool is open MATLAB software for analyzing eigenvalues, pseudospectra, and related spectral properties of matrices. It was developed at Oxford from 1999 - 2002 by Thomas Wright under the direction of Nick Trefethen. It is now maintained by Mark Embree.

Here is the pseudospectrum of our example matrix.

eigtool(A)

You can see the three eigenvalues at 1, 2 and 3. The contours outline the regions where the eigenvalues can move when the matrix is perturbed by a specified amount.

This animated gif adjusts the region around the eigenvalues at 2 and 3 until the regions coalesce in a saddle point near 2.4.

To see more precisely where they coalesce, we use the function `neardefmat`, nearest defective matrix, written by Michael Overton, a professor at NYU. The outputs from `neardefmat` are `B`, the nearest defective matrix, and `z`, the resulting double eigenvalue. By default `neardefmat` prints some useful details about the computation.

```
format long
[B,z] = neardefmat(A)
```

neardefmat: lower bound on distance is 0.000299834 and upper bound is 0.00618718 neardefmat: finished, lowest saddle found is 2.41449 + 0 i with dist = 0.000375171, mu = 0, |u'*v| = 3.67335e-06 singular vector residual norms are 1.8334e-13, 5.17298e-14 this saddle point was found in search # 3 distance to nearest defective matrix B found is 0.000375171 neardefmat: two closest computed eigenvalues of B differ by 0.00269716 and condition numbers of these eigenvalues are 272230, 272669 B = 1.0e+02 * -0.639999774916363 0.819999582124016 0.210002343527353 1.439999736641696 -1.779999511065701 -0.460002742035793 -7.710000068576393 9.620000127314574 2.479999285995845 z = 2.414492300748820

`B` is constructed from `z`, an SVD and a rank one perturbation of `A`.

I = eye(3); [U,S,V] = svd(A - z*I); format short sigma = S(3,3) u = U(:,3) v = V(:,3) format long B = A - sigma*u*v'

sigma = 3.7517e-04 u = -0.6373 0.7457 0.1942 v = 0.0941 -0.1748 0.9801 B = 1.0e+02 * -0.639999774916363 0.819999582124016 0.210002343527353 1.439999736641696 -1.779999511065701 -0.460002742035793 -7.710000068576393 9.620000127314574 2.479999285995845

This is essentially the same perturbation that I used to construct `A` in last week's blog post. The size of the perturbation is

sigma = S(3,3) sigma = norm(B - A)

sigma = 3.751705175569800e-04 sigma = 3.751705175556183e-04

By construction, `B` has a double eigenvalue, but the eigenvalue's condition is infinite, so it is not computed accurately.

eig(B)

ans = 2.417189453582737 2.414492295643361 1.168318252152021

Wilkinson was interested in the Frank family of matrices. Of course, we can find them in the `gallery`.

```
help private/frank
```

FRANK Frank matrix. F = GALLERY('FRANK',N, K) is the Frank matrix of order N. It is upper Hessenberg with determinant 1. If K = 1, the elements are reflected about the anti-diagonal (1,N)--(N,1). The eigenvalues of F may be obtained in terms of the zeros of the Hermite polynomials. They are positive and occur in reciprocal pairs; thus if N is odd, 1 is an eigenvalue. F has FLOOR(N/2) ill-conditioned eigenvalues---the smaller ones.

The 9x9 Frank matrix fits nicely on the page.

```
F = gallery('frank',9)
```

F = 9 8 7 6 5 4 3 2 1 8 8 7 6 5 4 3 2 1 0 7 7 6 5 4 3 2 1 0 0 6 6 5 4 3 2 1 0 0 0 5 5 4 3 2 1 0 0 0 0 4 4 3 2 1 0 0 0 0 0 3 3 2 1 0 0 0 0 0 0 2 2 1 0 0 0 0 0 0 0 1 1

The matrix is clearly far from symmetric. Since it is almost upper triangular you might think the eigenvalues are near the diagonal elements, but they are not. Here they are, sorted.

```
format long
[lambda,p] = sort(eig(F));
lambda
```

lambda = 0.044802721845067 0.082015890201792 0.162582729486328 0.374121140912600 0.999999999999966 2.672931012564413 6.150714797051293 12.192759202150013 22.320072505788492

And here are the corresponding condition numbers.

format short e kappa = condeig(F); kappa = kappa(p)

kappa = 1.4762e+04 2.5134e+04 1.1907e+04 1.5996e+03 6.8615e+01 3.6210e+00 1.5268e+00 2.3903e+00 1.9916e+00

The three smallest eigenvalues have the largest condition numbers. We expect the nearest defective matrix to bring two of them together.

Eigtool of the Frank matrix.

`eigtool(F)`

Focus on the three smallest eigenvalues. The red contour shows the regions containing the two that are coalescing.

Find the nearest defective matrix.

```
format long
[B,z] = neardefmat(F);
```

neardefmat: lower bound on distance is 3.4791e-07 and upper bound is 8.76938e-05 neardefmat: finished, lowest saddle found is 0.06104 + 0 i with dist = 5.01214e-07, mu = 0, |u'*v| = 6.09133e-10 singular vector residual norms are 3.6435e-15, 6.12511e-16 this saddle point was found in search # 1 distance to nearest defective matrix B found is 5.01214e-07 neardefmat: two closest computed eigenvalues of B differ by 1.7808e-06 and condition numbers of these eigenvalues are 3.62946e+08, 3.62937e+08

The two smallest eigenvalues have merged, but they are hyper-sensitive.

sort(eig(B))

ans = 0.061039353359070 0.061041134158948 0.168054858422710 0.373370423277969 1.000016557229929 2.672931171283966 6.150714794349185 12.192759202129578 22.320072505788680

If we're not too strict on the tolerance, the eigenvector matrix `X` is rank deficient, implying the matrix `B` is defective.

[X,~] = eig(B); rank(X,1.e-8)

ans = 8

We can get a picture of the perturbation with the `imagesc` function. First, the matrix itself. You can see its upper Hessenberg (nearly upper triangular) structure.

`imagesc(F)`

Now, the picture of the perturbation shows that it is concentrated in the upper right-hand corner.

`imagesc(F-B)`

What perturbation of a symmetric matrix makes its eigenvalues coalesce? A great example was provided by Wilkinson himself.

```
help wilkinson
W = wilkinson(9)
```

WILKINSON Wilkinson's eigenvalue test matrix. W = WILKINSON(n) is J. H. Wilkinson's eigenvalue test matrix, Wn+. It is a symmetric, tridiagonal matrix with pairs of nearly, but not exactly, equal eigenvalues. The most frequently used case is WILKINSON(21). W = WILKINSON(n,CLASSNAME) returns a matrix of class CLASSNAME, which can be either 'single' or 'double' (the default). Example: WILKINSON(7) is 3 1 0 0 0 0 0 1 2 1 0 0 0 0 0 1 1 1 0 0 0 0 0 1 0 1 0 0 0 0 0 1 1 1 0 0 0 0 0 1 2 1 0 0 0 0 0 1 3 Reference page in Doc Center doc wilkinson W = 4 1 0 0 0 0 0 0 0 1 3 1 0 0 0 0 0 0 0 1 2 1 0 0 0 0 0 0 0 1 1 1 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 1 1 1 0 0 0 0 0 0 0 1 2 1 0 0 0 0 0 0 0 1 3 1 0 0 0 0 0 0 0 1 4

`W` is symmetric and tridiagonal with ones on the off-diagonals. In order to have a multiple eigenvalue such a matrix has to have a pair of zeroes somewhere on the sub- and super-diagonal. But this matrix does not have any small off-diagonal elements. Nevertheless, the two largest eigenvalues are very close to each other.

lambda = eig(W)

lambda = -1.125422415673319 0.254718759825861 0.952584219075215 1.822717080887109 2.178284739549981 3.177282919112892 3.247396472578982 4.745281240174140 4.747156984469142

Since the matrix is symmetric, the left and right eigenvectors are orthogonal and transposes of each other. All the eigenvalues are perfectly well conditioned. Where should we look for the coalescing eigenvalues?

kappa = condeig(W)

kappa = 1.000000000000000 1.000000000000000 1.000000000000000 1.000000000000000 1.000000000000000 1.000000000000000 1.000000000000000 1.000000000000000 1.000000000000000

`Eigtool` of Wilkinson(9). There are only seven circles for a 9x9 matrix -- the two circles on the right each contain two eigenvalues.

`eigtool(W)`

Let's zoom way in on the rightmost circle.

The -3.04 contours are just barely kissing. This is a different kind of saddle point -- a nonsmooth one. You will have to study the Alam *et al* paper to see why this is so important.

Let's see what `neardefmat` has to say about Wilkinson's matrix.

[B,z] = neardefmat(W);

neardefmat: A is normal within tolerance distance to nearest defective matrix B is 0.000937872 neardefmat: two closest computed eigenvalues of B differ by 2.01677e-09 and condition numbers of these eigenvalues are 465036, 465036

Sure enough, `W` is symmetric and hence normal. A rather different algorithm is required.

The perturbed eigenvalues have small imaginary parts.

eig(B)

ans = -1.125422415673318 + 0.000000000000000i 0.254718759825861 + 0.000000000000000i 0.952584219075214 + 0.000000000000000i 1.822717080887109 + 0.000000000000000i 2.178284739549983 + 0.000000000000000i 3.177282919112890 + 0.000000000000000i 3.247396472578982 + 0.000000000000000i 4.746219112321641 + 0.000000001008386i 4.746219112321641 - 0.000000001008386i

Finally, here are images of the symmetric, tridiagonal matrix W and the very nonsymmetric perturbation `W-B` that brings two of its eigenvalues together.

Get
the MATLAB code

Published with MATLAB® R2018b

On May 29-30, I plan to attend a conference, organized by Nick Higham, at the University of Manchester. The title of the conference is... read more >>

]]>On May 29-30, I plan to attend a conference, organized by Nick Higham, at the University of Manchester. The title of the conference is

**Celebrating the Centenary of James H. Wilkinson's Birth**

I am giving a talk about one of Wilkinson's favorite topics, how can perturbations of a matrix with sensitive eigenvalues produce a defective matrix with a multiple eigenvalue. I uncovered this example from my original Fortran MATLAB User's Guide, a University of New Mexico technical report, dated 1981, four years before MathWorks. The text of this blog post is copied from that guide. I've inserted comments where today's MATLAB scales the eigenvectors differently.

In this example, we construct a matrix whose eigenvalues are moderately sensitive to perturbations and then analyze that sensitivity. We begin with the statement

B = [3 0 7; 0 2 0; 0 0 1]

B = 3 0 7 0 2 0 0 0 1

Obviously, the eigenvalues of B are 1, 2 and 3 . Moreover, since B is not symmetric, these eigenvalues are slightly sensitive to perturbation. (The value B(1,3) = 7 was chosen so that the elements of the matrix A below are less than 1000.)

We now generate a similarity transformation to disguise the eigenvalues and make them more sensitive.

L = [1 0 0; 2 1 0; -3 4 1] M = L\L'

L = 1 0 0 2 1 0 -3 4 1 M = 1 2 -3 -2 -3 10 11 18 -48

The matrix M has determinant equal to 1 and is moderately badly conditioned. The similarity transformation is

A = M*B/M

A = -64.0000 82.0000 21.0000 144.0000 -178.0000 -46.0000 -771.0000 962.0000 248.0000

Because det(M) = 1 , the elements of A would be exact integers if there were no roundoff. So,

A = round(A)

A = -64 82 21 144 -178 -46 -771 962 248

This, then, is our test matrix. We can now forget how it was generated and analyze its eigenvalues.

[X,D] = eig(A)

X = 0.0891 0.0735 -0.1089 -0.1782 -0.1923 0.1634 0.9800 0.9786 -0.9805 D = 3.0000 0 0 0 1.0000 0 0 0 2.0000

% { Today the eigenvectors are scaled to have unit length. % Classic MATLAB did not scale the eigenvectors. It got % % X = % % -.0891 3.4903 41.8091 % .1782 -9.1284 -62.7136 % -.9800 46.4473 376.2818 % }

Since A is similar to B, its eigenvalues are also 1, 2 and 3. They happen to be computed in another order by the EISPACK subroutines. The fact that the columns of X, which are the eigenvectors, are so far from being orthonormal is our first indication that the eigenvalues are sensitive. To see this sensitivity, we display more figures of the computed eigenvalues.

```
format long
diag(D)
```

ans = 3.000000000003868 0.999999999998212 1.999999999997978

We see that, on this computer, the last five {today: four} significant figures are contaminated by roundoff error. A somewhat superficial explanation of this is provided by

```
format short
cond(X)
```

ans = 1.7690e+03

% { Classic: % % ANS = % % 3.2216e+05 % }

The condition number of X gives an upper bound for the relative error in the computed eigenvalues. However, this condition number is affected by scaling.

X = X/diag(X(3,:)) cond(X)

X = 0.0909 0.0751 0.1111 -0.1818 -0.1965 -0.1667 1.0000 1.0000 1.0000 ans = 1.7692e+03

Rescaling the eigenvectors so that their last components are all equal to one has two consequences. The condition of X is decreased by over two orders of magnitude. {Not today.} (This is about the minimum condition that can be obtained by such diagonal scaling.) Moreover, it is now apparent that the three eigenvectors are nearly parallel.

More detailed information on the sensitivity of the individual eigenvalues involves the left eigenvectors.

Y = inv(X') Y'*A*X

Y = -511.5000 259.5000 252.0000 616.0000 -346.0000 -270.0000 159.5000 -86.5000 -72.0000 ans = 3.0000 -0.0000 -0.0000 0.0000 1.0000 0.0000 0 -0.0000 2.0000

We are now in a position to compute the sensitivities of the individual eigenvalues.

for j = 1:3, c(j) = norm(Y(:,j))*norm(X(:,j)); end c

c = 833.1092 450.7228 383.7564

These three numbers are the reciprocals of the cosines of the angles between the left and right eigenvectors. It can be shown that perturbation of the elements of A can result in a perturbation of the j-th eigenvalue which is c(j) times as large. In this example, the first eigenvalue has the largest sensitivity.

We now proceed to show that A is close to a matrix with a double eigenvalue. The direction of the required perturbation is given by

E = -1.e-6*Y(:,1)*X(:,1)'

E = 1.0e-03 * 0.0465 -0.0930 0.5115 -0.0560 0.1120 -0.6160 -0.0145 0.0290 -0.1595

With some trial and error which we do not show, we bracket the point where two eigenvalues of a perturbed A coalesce and then become complex.

eig(A + .4*E) eig(A + .5*E)

ans = 1.1500 2.5996 2.2504 ans = 2.4067 + 0.1753i 2.4067 - 0.1753i 1.1866 + 0.0000i

Now, a bisecting search, driven by the imaginary part of one of the eigenvalues, finds the point where two eigenvalues are nearly equal.

r = .4; s = .5; while s-r > 1.e-14 t = (r+s)/2; d = eig(A+t*E); if imag(d(1)) == 0, r = t; else, s = t; end end format long t

t = 0.450380734135428

Finally, we display the perturbed matrix, which is obviously close to the original, and its pair of nearly equal eigenvalues.

A+t*E eig(A+t*E)

ans = 1.0e+02 * -0.639999790572959 0.819999581145917 0.210002303697455 1.439999747786789 -1.779999495573578 -0.460002774345322 -7.710000065305207 9.620000130610412 2.479999281642729 ans = 2.415743144226897 2.415738627741217 1.168517777651195

The first two eigenvectors of A + t*E are almost indistinguishable indicating that the perturbed matrix is almost defective.

```
[X,D] = eig(A+t*E)
format short
cond(X)
```

X = 0.094108719644215 -0.094108788388666 -0.070056238584537 -0.174780805492871 0.174780755585658 0.194872753681838 0.980099596427929 -0.980099598727050 -0.978323429806240 D = 2.415743144226897 0 0 0 2.415738627741217 0 0 0 1.168517777651195 ans = 3.9853e+08

Get
the MATLAB code

Published with MATLAB® R2018b

Let me tell you about MATLAB's controversial function `rat`.... read more >>

Let me tell you about MATLAB's controversial function `rat`.

The function `rat` is older than MathWorks. It was one of the 73 functions in my original Fortran MATLAB. Here is the help entry from 40 years ago.

```
<>
help rat
```

RAT An experimental function which attempts to remove the roundoff error from results that should be "simple" rational numbers. RAT(X) approximates each element of X by a continued fraction of the form

a/b = d1 + 1/(d2 + 1/(d3 + ... + 1/dk))

with k <= len, integer di and abs(di) <= max . The default values of the parameters are len = 5 and max = 100. RAT(len,max) changes the default values. Increasing either len or max increases the number of possible fractions. <A,B> = RAT(X) produces integer matrices A and B so that

A ./ B = RAT(X)

Some examples:

long T = hilb(6), X = inv(T) <A,B> = rat(X) H = A ./ B, S = inv(H)

short e d = 1:8, e = ones(d), A = abs(d'*e - e'*d) X = inv(A) rat(X) display(ans)

Let's try that second example with modern MATLAB and `format rat`.

format short d = 1:6 e = ones(size(d)), A = abs(d'*e - e'*d) format short e X = inv(A) format rat X

d = 1 2 3 4 5 6 e = 1 1 1 1 1 1 A = 0 1 2 3 4 5 1 0 1 2 3 4 2 1 0 1 2 3 3 2 1 0 1 2 4 3 2 1 0 1 5 4 3 2 1 0 X = -4.0000e-01 5.0000e-01 8.3267e-17 -3.3307e-17 -9.9920e-17 1.0000e-01 5.0000e-01 -1.0000e+00 5.0000e-01 3.6978e-33 8.8818e-17 -5.5511e-17 8.3267e-17 5.0000e-01 -1.0000e+00 5.0000e-01 1.4063e-16 -9.2519e-17 -3.3307e-17 3.6978e-33 5.0000e-01 -1.0000e+00 5.0000e-01 5.5511e-17 -9.9920e-17 8.8818e-17 1.4063e-16 5.0000e-01 -1.0000e+00 5.0000e-01 1.0000e-01 -5.5511e-17 -9.2519e-17 5.5511e-17 5.0000e-01 -4.0000e-01 X = Columns 1 through 5 -2/5 1/2 * * * 1/2 -1 1/2 * * * 1/2 -1 1/2 * * * 1/2 -1 1/2 * * * 1/2 -1 1/10 * * * 1/2 Column 6 1/10 * * * 1/2 -2/5

Which version of `X` do you prefer? The first, printed with a floating point format so the roundoff errors are readily apparent. Or the second, printed with a rational format that seeks to disguise the roundoff. For this example, almost everybody would vote for the rational format.

But I have learned that it is better to explain floating point arithmetic than it is to try to hide it. Besides, most important MATLAB calculations involve larger matrices, and eigenvalues, differential equations and signal processing. They don't have such tidy results.

So, we don't have much use for `rat` and `format rat` today, but they're still there and, as I said, they're moderately controversial.

The core of `rat` generates continued fractions by repeatedly subtracting off the integer part and taking the reciprocal of what is left. If there is nothing left to reciprocate, the process terminates because the input is a rational number. If the input is irrational, the process does not terminate.

The controversary comes from what is meant by integer part. To generate proper continued fractions, integer part should be `floor`. That always leaves a positive fraction to reciprocate. But 40 years ago, I got clever and used `round` instead of `floor`. That means it may take fewer terms to obtain a specified accuracy, but the continued fractions are, shall I say, *unorthodox*. Let's see more examples.

Let's generate increasingly accurate rational approximations to $\pi$. With two output arguments, `rat` unwinds the continued fraction to produce two integers whose ratio has the same value. The accuracy of `rat` is determined by an optional tolerance, `rat(x,tol)`. The default tolerance is `1.0e-6*norm(x)`.

Requesting low accuracy produces a familiar approximation of $\pi$.

```
form = ' tol = %6.1e \n r = %s \n = %d/%d\n';
tol = 1.0e-2;
r = rat(pi,tol);
[n,d] = rat(pi,tol);
fprintf(form,tol,r,n,d)
```

tol = 1.0e-02 r = 3 + 1/(7) = 22/7

The default tolerance produces a less familiar, but elegant approximation.

tol = 1.0e-6; r = rat(pi,tol); [n,d] = rat(pi,tol); fprintf(form,tol,r,n,d)

tol = 1.0e-06 r = 3 + 1/(7 + 1/(16)) = 355/113

Somewhat higher accuracy produces a negative term.

tol = 1.0e-8; r = rat(pi,tol); [n,d] = rat(pi,tol); fprintf(form,tol,r,n,d)

tol = 1.0e-08 r = 3 + 1/(7 + 1/(16 + 1/(-294))) = 104348/33215

Another negative term.

tol = 1.0e-11; r = rat(pi,tol); [n,d] = rat(pi,tol); fprintf(form,tol,r,n,d)

tol = 1.0e-11 r = 3 + 1/(7 + 1/(16 + 1/(-294 + 1/(3 + 1/(-4))))) = 1146408/364913

If we're careful with the curly braces, Latex and MathJax can typeset that last continued fraction and produce a graphic for this post.

{{1}\over{{7+ {{1}\over{{16+ {{1}\over{{-294+ {{1}\over{{3+ {{1}\over{-4}}}}}}}}}}}}}}

$${{1}\over{{7+ {{1}\over{{16+ {{1}\over{{-294+ {{1}\over{{3+ {{1}\over{-4}}}}}}}}}}}}}}$$

I have another function, `ratp`, that is `rat` with `round` replaced by `floor` so that it produces proper, but longer, continued fractions with positive terms. Let's see how it compares to `rat`. Since the MATLAB quantity `pi` that approximates $\pi$ is rational, we can set `tol` to zero without fear of an infinite loop.

form = ' tol = %6.1e \n r = %s \n p = %s \n\n'; for tol = 10.^[-2 -6 -8 -11 -inf] fprintf(form,tol,rat(pi,tol),ratp(pi,tol)) end

tol = 1.0e-02 r = 3 + 1/(7) p = 3 + 1/(7) tol = 1.0e-06 r = 3 + 1/(7 + 1/(16)) p = 3 + 1/(7 + 1/(15 + 1/(1))) tol = 1.0e-08 r = 3 + 1/(7 + 1/(16 + 1/(-294))) p = 3 + 1/(7 + 1/(15 + 1/(1 + 1/(292)))) tol = 1.0e-11 r = 3 + 1/(7 + 1/(16 + 1/(-294 + 1/(3 + 1/(-4))))) p = 3 + 1/(7 + 1/(15 + 1/(1 + 1/(292 + 1/(1 + 1/(1 + 1/(1 + 1/(2)))))))) tol = 0.0e+00 r = 3 + 1/(7 + 1/(16 + 1/(-294 + 1/(3 + 1/(-4 + 1/(5 + 1/(-15))))))) p = 3 + 1/(7 + 1/(15 + 1/(1 + 1/(292 + 1/(1 + 1/(1 + 1/(1 + 1/(2 + 1/(1 + 1/(3 + 1/(1 + 1/(14))))))))))))

The world's most irrational number, as measured by the length of continued fraction approximations, is $\phi$, the golden ratio. (Proof: solve $x = 1 + 1/x$.) The coefficients in the proper continued fraction are all ones, thereby achieving the slowest possible convergence. For the same accuracy, `ratp(phi,tol)` requires about twice as many terms as `rat(phi,tol)`.

```
format short
phi = (1 + sqrt(5))/2
tol = 1.0e-4
r = rat(phi)
p = ratp(phi)
```

phi = 1.6180 tol = 1.0000e-04 r = '2 + 1/(-3 + 1/(3 + 1/(-3 + 1/(3 + 1/(-3 + 1/(3 + 1/(-3)))))))' p = '1 + 1/(1 + 1/(1 + 1/(1 + 1/(1 + 1/(1 + 1/(1 + 1/(1 + 1/(1 + 1/(1 + 1/(1 + 1/(1 + 1/(1 + 1/(1 + 1/(1))))))))))))))'

tol = 0; lengthr = length(rat(phi,tol)) lengthp = length(ratp(phi,tol))

lengthr = 154 lengthp = 297

(For my international readers, DIY stands for "Do It Yourself".) Create your own `ratp` by making a one-word edit in a copy of `toolbox\matlab\specfun\rat.m`. Compare the two functions on `sqrt(2)`. Warning: all is not what appears at first glance.

Get
the MATLAB code

Published with MATLAB® R2018b

*This post is by my colleague Cosmin Ionita.*... read more >>

*This post is by my colleague Cosmin Ionita.*

The `'makima'` cubic interpolation method was recently introduced in MATLAB® in the R2017b release as a new option in `interp1`, `interp2`, `interp3`, `interpn`, and `griddedInterpolant`. Its implementation is not user visible; thus, we have been receiving inquiries from our users about the specifics of this new cubic method.

In the following, we address our users' inquiries by answering these questions:

- What is the
`'makima'`formula? - How does
`'makima'`compare with MATLAB's other cubic methods?

In a nutshell, `'makima'` is short for *modified Akima piecewise cubic Hermite interpolation*. It represents a MATLAB-specific modification of Akima's derivative formula and has the following key properties:

- It produces undulations which find a nice middle ground between
`'spline'`and`'pchip'`. - It is a local cubic interpolant which generalizes to 2-D grids and higher-dimensional n-D grids.
- It increases the robustness of Akima's formula in the edge case of equal side slopes.
- It eliminates a special type of overshoot arising when the data is constant for more than two consecutive nodes.

For each interval $[x_i~x_{i+1})$ of an input data set of nodes $x$ and values $v$, piecewise cubic Hermite interpolation finds a cubic polynomial which not only interpolates the given data values $v_i$ and $v_{i+1}$ at the interval's nodes $x_i$ and $x_{i+1}$, but also has specific derivatives $d_i$ and $d_{i+1}$ at $x_i$ and $x_{i+1}$. For more details we refer to the interpolation chapter in Cleve's *Numerical Computing with MATLAB* textbook.

The key to cubic Hermite interpolation is the choice of derivatives $d_i$.

In his 1970 paper, Hiroshi Akima introduced a derivative formula which *avoids excessive local undulations:*

*H. Akima, "A New Method of Interpolation and Smooth Curve Fitting Based on Local Procedures", JACM, v. 17-4, p.589-602, 1970.*

Let $\delta_i =\left(v_{i+1} -v_i \right)/\left(x_{i+1} -x_i \right)$ be the slope of the interval $[x_i~x_{i+1})$. Akima's derivative at $x_i$ is defined as:

$$d_i = \frac{|\delta_{i+1} - \delta_i|\delta_{i-1} + |\delta_{i-1} - \delta_{i-2}|\delta_i} {|\delta_{i+1} - \delta_i| + |\delta_{i-1} - \delta_{i-2}|},$$

and represents a weighted average between the slopes $\delta_{i-1}$ and $\delta_i$ of the intervals $[x_{i-1}~x_i)$ and $[x_i~x_{i+1})$:

$$d_i = \frac{w_1}{w_1+w_2}\delta_{i-1} + \frac{w_2}{w_1+w_2}\delta_i \,, \quad\quad w_1 = |\delta_{i+1} - \delta_i|, \quad\quad w_2 = |\delta_{i-1} - \delta_{i-2}|.$$

Notice that Akima's derivative at $x_i$ is computed locally from the five points $x_{i-2}$, $x_{i-1}$, $x_i$, $x_{i+1}$, and $x_{i+2}$. For the end points $x_1$ and $x_n$, it requires the slopes $\delta_{-1}, \delta_0$ and $\delta_n ,\delta_{n+1}$. Since these slopes are not available in the input data, Akima proposed using quadratic extrapolation to compute them as $\delta_0 = 2\delta_1 -\delta_2$, $\delta_{-1} = 2\delta_0 -\delta_1$ and $\delta_n = 2\delta_{n-1} -\delta_{n-2}$, $\delta_{n+1} = 2\delta_n -\delta_{n-1}$.

But wait!

MATLAB already has two cubic Hermite interpolation methods (see Cleve's blog Splines and Pchips):

`'spline'`computes the derivatives by imposing the constraint of continuous second derivatives (this guarantees a very smooth interpolation result),`'pchip'`computes the derivatives by imposing local monotonicity in each interval $[x_i~x_{i+1})$ (this preserves the shape of the data).

How does Akima's formula compare with `'spline'` and `'pchip'`?

Akima's formula finds a nice middle ground between `'spline'` and `'pchip'`:

- Its undulations (or wiggles) have smaller amplitudes than
`'spline'`. - It is not as aggressive as
`'pchip'`in reducing the undulations.

Here's a representative example comparing these three cubic Hermite interpolants:

```
x1 = [1 2 3 4 5 5.5 7 8 9 9.5 10];
v1 = [0 0 0 0.5 0.4 1.2 1.2 0.1 0 0.3 0.6];
xq1 = 0.75:0.05:10.25;
compareCubicPlots(x1,v1,xq1,false,'NE')
```

Akima's formula certainly produces nice results. However, we are not yet ready to fully adopt it.

**Edge case of equal side slopes**

When the lower slopes are equal, $\delta_{i-2} =\delta_{i-1}$, and the upper slopes are equal, $\delta_i =\delta_{i+1}$, both the numerator and denominator become equal to 0 and Akima's formula returns `NaN`. Akima himself recognized this problem and proposed taking the average of the lower and upper slopes for this edge case: $d_i =\left(\delta_{i-1} +\delta_i \right)/2$.

Let's try this fix for the following data set where the intervals $[3~4)$ and $[4~5)$ have equal slopes $\delta_3 =\delta_4 =1$ and the intervals $[5~6)$ and $[6~7)$ also have equal slopes $\delta_5 =\delta_6 =0$:

```
x2 = 1:8;
v2 = [-1 -1 -1 0 1 1 1 1];
xq2 = 0.75:0.05:8.25;
compareCubicPlots(x2,v2,xq2,false,'SE')
ylim([-1.2 1.2])
```

In this case, Akima's `NaN` derivative at $x_5=5$ gets replaced with the average slope $d_5 =\left(\delta_4 +\delta_5 \right)/2=0.5$.

But we are still not ready!

When should we switch from Akima's formula to the averaging formula? What if the slopes are *almost* equal?

Let's change the above example such that $\delta_5 =\varepsilon$ and $\delta_6 =-\varepsilon$ by setting $v_6 =1+\varepsilon$ , where $\varepsilon =2^{-52}$ represents machine epsilon. We can now compare Akima's full formula with the averaging formula at $d_5$ for two data sets that differ only by $\varepsilon$:

% v2eps: new data set which differs only by an eps in v2(6) v2eps = v2; v2eps(6) = v2eps(6) + eps; plot(x2,v2,'ko','LineWidth',2,'MarkerSize',10,'DisplayName','Input data v2') hold on plot(xq2,akima(x2,v2,xq2),'Color',[1 2/3 0],'LineWidth',2,... 'DisplayName','v2(6) = 1 uses averaging formula') plot(xq2,akima(x2,v2eps,xq2),'-.','Color',[1 2/3 0],'LineWidth',2,... 'DisplayName','v2(6) = 1+eps uses Akima''s formula') hold off legend('Location','SE') title('Akima interpolant for (almost) equal side slopes' ) ylim([-1.2 1.2])

The averaging formula (solid line) returns the derivative $d_5=0.5$ while Akima's formula (dashed line) returns $d_5=1$. Therefore, there is an unwanted non-negligible difference at around $x_5=5$ between the two Akima interpolants, even though the two underlying data sets differ only by $\varepsilon$ at $v_6$.

**Requirement (1)**

We do not want our Akima interpolant to switch to a different formula for edge cases.

**Overshoot if data is constant for more than two consecutive nodes**

The previous example was strategically chosen to also reveal another property of the Akima interpolant: if the data is constant for more than two consecutive nodes, like $v_5=v_6=v_7=1$ in the previous example, then the Akima interpolant may produce an overshoot, namely, the interpolant's undulation in the $[5~6)$ interval above.

This special type of overshoot is not desirable in many engineering applications.

The same concerns hold for the undershoot present in the $[2~3)$ interval above.

**Requirement (2)**

We do not want our Akima interpolant to produce overshoot (or undershoot) when the data is constant for more than two consecutive nodes.

In MATLAB, `'makima'` cubic Hermite interpolation addresses requirements (1) and (2) outlined above.

To eliminate overshoot and avoid edge cases of both numerator and denominator being equal to 0, we modify Akima's derivative formula by tweaking the weights $w_1$ and $w_2$ of the slopes $\delta_{i-1}$ and $\delta_i$:

$$ d_i = \frac{w_1}{w_1+w_2}\delta_{i-1} + \frac{w_2}{w_1+w_2}\delta_i \,, \quad\quad w_1 = |\delta_{i+1} - \delta_i| + |\delta_{i+1} + \delta_i|/2, \quad\quad w_2 = |\delta_{i-1} - \delta_{i-2}| + |\delta_{i-1} + \delta_{i-2}| / 2. $$

**This modified formula represents the 'makima' derivative formula used in MATLAB:**

- The addition of the $\left|\delta_{i+1} +\delta_i \right|/2$ and $\left|\delta_{i-1} +\delta_{i-2} \right|/2$ terms forces $d_i =0$ when $\delta_i =\delta_{i+1} =0$, i.e., $d_i =0$ when $v_i =v_{i+1} =v_{i+2}$. Therefore, it eliminates overshoot when the data is constant for more than two consecutive nodes.
- These new terms also address the edge case of equal side slopes discussed above, $\delta_{i-2} =\delta_{i-1}$ and $\delta_i =\delta_{i+1}$. However, there is one case which slips through: if the data is constant $v_{i-2} = v_{i-1} = v_i = v_{i+1} = v_{i+2}$, then the four slopes are all equal to zero and we get $d_i = \mathrm{NaN}$. For this special case of constant data, we set $d_i =0$.

Let's try the `'makima'` formula on the above overshoot example:

```
compareCubicPlots(x2,v2,xq2,true,'SE')
ylim([-1.2 1.2])
```

Indeed, `'makima'` does not produce an overshoot if the data is constant for more than two nodes ($v_5=v_6=v_7=1$ above).

But what does this mean for the undulations we saw in our first example?

```
compareCubicPlots(x1,v1,xq1,true,'NE')
```

Notice that `'makima'` closely follows the result obtained with Akima's formula. In fact, the results are so similar that it is hard to tell them apart on the plot.

Therefore, `'makima'` still preserves Akima's desirable properties of being a nice middle ground between `'spline'` and `'pchip'` in terms of the resulting undulations.

But we are still not done!

Akima's formula and our modified `'makima'` formula have another desirable property: they generalize to higher dimensional n-D gridded data. Akima noticed this property in his 1974 paper.

*H. Akima, "A Method of Bivariate Interpolation and Smooth Surface Fitting Based on Local Procedures", Communications of the ACM, v. 17-1, p.18-20, 1974.*

For example, to interpolate on a 2-D grid $\left(x,y\right)$, we first apply the `'makima'` derivative formula separately along $x$ and $y$ to obtain two directional derivatives for each grid node. Then, we also compute the cross-derivative along $xy$ for each grid node. The cross-derivative formula first computes the 2-D divided differences corresponding to the 2-D grid data and applies the `'makima'` derivative formula along $x$ on these differences; then, it takes the result and applies the derivative formula along $y$. The derivatives and cross-derivatives are then plugged in as coefficients of a two-variable cubic Hermite polynomial representing the 2-D interpolant.

The same procedure applies to higher dimensional n-D grids with $n\ge 2$ and requires computing cross-derivatives along all possible cross-directions. Therefore, `'makima'` is supported not only in `interp1`, but also in `interp2`, `interp3`, `interpn`, and `griddedInterpolant`.

Here's how 2-D `'makima'` interpolation compares with 2-D `'spline'` interpolation on gridded data generated with the `peaks` function:

[X1,Y1,V1] = peaks(5); [Xq1,Yq1] = meshgrid(-3:0.1:3,-3:0.1:3); Vqs1 = interp2(X1,Y1,V1,Xq1,Yq1,'spline'); surf(Xq1,Yq1,Vqs1) axis tight title('2-D ''spline''') snapnow Vqm1 = interp2(X1,Y1,V1,Xq1,Yq1,'makima'); surf(Xq1,Yq1,Vqm1) axis tight title('2-D ''makima''')

Notice the smaller undulations (or wiggles) generated by `'makima'`.

Finally, let's try `'makima'` on an example with a few 2-D peaks where the data has sharp edges and steps:

V2 = zeros(10); V2(2:5,2:5) = 3/7; % First step V2(6:7,6:7) = [1/4 1/5; 1/5 1/4]; % Middle step V2(8:9,8:9) = 1/2; % Last step V2 = flip(V2,2); [Xq2,Yq2] = meshgrid(1:0.2:10,1:0.2:10); Vqs2 = interp2(V2,Xq2,Yq2,'spline'); surf(Xq2,Yq2,Vqs2) axis tight title('2-D ''spline''') snapnow Vqm2 = interp2(V2,Xq2,Yq2,'makima'); surf(Xq2,Yq2,Vqm2) axis tight title('2-D ''makima''')

The `'makima'` result has smaller undulations (or wiggles) than `'spline'.`

We conclude with a summary of properties of interest for the three cubic Hermite interpolation methods supported in MATLAB:

summary = table( ... ["C2"; "Uses all points"; "Not-a-knot condition"; "Yes"], ... ["C1"; "Uses 4 points"; "One-sided slope"; "No"], ... ["C1"; "Uses 5 points"; "Quadratic extrapolation"; "Yes"], ... 'VariableNames', ... ["spline" "pchip" "makima"], ... 'RowNames', ... ["Continuity"; "Locality"; "End points"; "Supports n-D"]); disp(summary)

spline pchip makima ______________________ _________________ _________________________ Continuity "C2" "C1" "C1" Locality "Uses all points" "Uses 4 points" "Uses 5 points" End points "Not-a-knot condition" "One-sided slope" "Quadratic extrapolation" Supports n-D "Yes" "No" "Yes"

function compareCubicPlots(x,v,xq,showMakima,legendLocation) % Plot 'pchip', 'spline', Akima, and 'makima' interpolation results vqp = pchip(x,v,xq); % same as vq = interp1(x,v,xq,'pchip') vqs = spline(x,v,xq); % same as vq = interp1(x,v,xq,'spline') vqa = akima(x,v,xq); if showMakima vqm = makima(x,v,xq); % same as vq = interp1(x,v,xq,'makima') end plot(x,v,'ko','LineWidth',2,'MarkerSize',10,'DisplayName','Input data') hold on plot(xq,vqp,'LineWidth',4,'DisplayName','''pchip''') plot(xq,vqs,'LineWidth',2,'DisplayName','''spline''') plot(xq,vqa,'LineWidth',2,'DisplayName','Akima''s formula') if showMakima plot(xq,vqm,'--','Color',[0 3/4 0],'LineWidth',2, ... 'DisplayName','''makima''') end hold off xticks(x(1)-1:x(end)+1) legend('Location',legendLocation) title('Cubic Hermite interpolation') end

function vq = akima(x,v,xq) %AKIMA Akima piecewise cubic Hermite interpolation. % % vq = AKIMA(x,v,xq) interpolates the data (x,v) at the query points xq % using Akima's cubic Hermite interpolation formula. % % References: % % H. Akima, "A New Method of Interpolation and Smooth Curve Fitting % Based on Local Procedures", JACM, v. 17-4, p.589-602, 1970. % % AKIMA vs. PCHIP vs. SPLINE: % % - Akima's cubic formula is a middle ground between SPLINE and PCHIP: % It has lower-amplitude wiggles than SPLINE, but is not as agressive % at reducing the wiggles as PCHIP. % - Akima's cubic formula and SPLINE generalize to n-D grids. % % Example: Compare AKIMA with MAKIMA, SPLINE, and PCHIP % % x = [1 2 3 4 5 5.5 7 8 9 9.5 10]; % v = [0 0 0 0.5 0.4 1.2 1.2 0.1 0 0.3 0.6]; % xq = 0.75:0.05:10.25; % vqs = spline(x,v,xq); % vqp = pchip(x,v,xq); % vqa = akima(x,v,xq); % vqm = makima(x,v,xq); % % figure % plot(x,v,'ko','LineWidth',2,'MarkerSize',10), hold on % plot(xq,vqp,'LineWidth',4) % plot(xq,vqs,xq,vqa,xq,vqm,'--','LineWidth',2) % legend('Data','''pchip''','''spline''','Akima''s formula','''makima''') % title('Cubic Hermite interpolation in MATLAB') % % See also MAKIMA, SPLINE, PCHIP. % mailto: cosmin.ionita@mathworks.com assert(isvector(x) && isvector(v) && (numel(x) == numel(v))) assert(numel(x) > 2) x = x(:).'; v = v(:).'; % Same shapes as in pchip.m and spline.m % Compute Akima slopes h = diff(x); delta = diff(v)./h; slopes = akimaSlopes(delta); % Evaluate the piecewise cubic Hermite interpolant vq = ppval(pwch(x,v,slopes,h,delta),xq); end

function s = akimaSlopes(delta) % Derivative values for Akima cubic Hermite interpolation % Akima's derivative estimate at grid node x(i) requires the four finite % differences corresponding to the five grid nodes x(i-2:i+2). % % For boundary grid nodes x(1:2) and x(n-1:n), append finite differences % which would correspond to x(-1:0) and x(n+1:n+2) by using the following % uncentered difference formula correspondin to quadratic extrapolation % using the quadratic polynomial defined by data at x(1:3) % (section 2.3 in Akima's paper): n = numel(delta) + 1; % number of grid nodes x delta_0 = 2*delta(1) - delta(2); delta_m1 = 2*delta_0 - delta(1); delta_n = 2*delta(n-1) - delta(n-2); delta_n1 = 2*delta_n - delta(n-1); delta = [delta_m1 delta_0 delta delta_n delta_n1]; % Akima's derivative estimate formula (equation (1) in the paper): % % H. Akima, "A New Method of Interpolation and Smooth Curve Fitting % Based on Local Procedures", JACM, v. 17-4, p.589-602, 1970. % % s(i) = (|d(i+1)-d(i)| * d(i-1) + |d(i-1)-d(i-2)| * d(i)) % / (|d(i+1)-d(i)| + |d(i-1)-d(i-2)|) weights = abs(diff(delta)); weights1 = weights(1:n); % |d(i-1)-d(i-2)| weights2 = weights(3:end); % |d(i+1)-d(i)| delta1 = delta(2:n+1); % d(i-1) delta2 = delta(3:n+2); % d(i) weights12 = weights1 + weights2; s = (weights2./weights12) .* delta1 + (weights1./weights12) .* delta2; % To avoid 0/0, Akima proposed to average the divided differences d(i-1) % and d(i) for the edge case of d(i-2) = d(i-1) and d(i) = d(i+1): ind = weights1 == 0 & weights2 == 0; s(ind) = (delta1(ind) + delta2(ind)) / 2; end

function vq = makima(x,v,xq) %MAKIMA Modified Akima piecewise cubic Hermite interpolation. % % We modify Akima's formula to eliminate overshoot and undershoot % when the data is constant for more than two consecutive nodes. % % vq = MAKIMA(x,v,xq) interpolates the data (x,v) at the query points xq % using a modified Akima cubic Hermite interpolation formula. % % References: % % H. Akima, "A New Method of Interpolation and Smooth Curve Fitting % Based on Local Procedures", JACM, v. 17-4, p.589-602, 1970. % % MAKIMA vs. PCHIP vs. SPLINE: % % - MAKIMA is a middle ground between SPLINE and PCHIP: % It has lower-amplitude wiggles than SPLINE, but is not as agressive % at reducing the wiggles as PCHIP. % - MAKIMA and SPLINE generalize to n-D grids. % % Example: No overshoot and undershoot with MAKIMA when the data is % constant for more than two consecutive nodes % % x = 1:7; % v = [-1 -1 -1 0 1 1 1]; % xq = 0.75:0.05:7.25; % vqs = spline(x,v,xq); % vqp = pchip(x,v,xq); % vqa = akima(x,v,xq); % vqm = makima(x,v,xq); % % figure % plot(x,v,'ko','LineWidth',2,'MarkerSize',10), hold on % plot(xq,vqp,'LineWidth',4) % plot(xq,vqs,xq,vqa,xq,vqm,'LineWidth',2) % legend('Data','''pchip''','''spline''','Akima''s formula',... % '''makima''','Location','SE') % title('''makima'' has no overshoot in x(1:3) and x(5:7)') % % Example: Compare MAKIMA with AKIMA, SPLINE, and PCHIP % % x = [1 2 3 4 5 5.5 7 8 9 9.5 10]; % v = [0 0 0 0.5 0.4 1.2 1.2 0.1 0 0.3 0.6]; % xq = 0.75:0.05:10.25; % vqs = spline(x,v,xq); % vqp = pchip(x,v,xq); % vqa = akima(x,v,xq); % vqm = makima(x,v,xq); % % figure % plot(x,v,'ko','LineWidth',2,'MarkerSize',10), hold on % plot(xq,vqp,'LineWidth',4) % plot(xq,vqs,xq,vqa,xq,vqm,'--','LineWidth',2) % legend('Data','''pchip''','''spline''','Akima''s formula','''makima''') % title('Cubic Hermite interpolation in MATLAB') % % See also AKIMA, SPLINE, PCHIP. % mailto: cosmin.ionita@mathworks.com assert(isvector(x) && isvector(v) && (numel(x) == numel(v))) assert(numel(x) > 2) x = x(:).'; v = v(:).'; % Same shapes as in pchip.m and spline.m % Compute modified Akima slopes h = diff(x); delta = diff(v)./h; slopes = makimaSlopes(delta); % Evaluate the piecewise cubic Hermite interpolant vq = ppval(pwch(x,v,slopes,h,delta),xq); end

function s = makimaSlopes(delta) % Derivative values for modified Akima cubic Hermite interpolation % Akima's derivative estimate at grid node x(i) requires the four finite % differences corresponding to the five grid nodes x(i-2:i+2). % % For boundary grid nodes x(1:2) and x(n-1:n), append finite differences % which would correspond to x(-1:0) and x(n+1:n+2) by using the following % uncentered difference formula correspondin to quadratic extrapolation % using the quadratic polynomial defined by data at x(1:3) % (section 2.3 in Akima's paper): n = numel(delta) + 1; % number of grid nodes x delta_0 = 2*delta(1) - delta(2); delta_m1 = 2*delta_0 - delta(1); delta_n = 2*delta(n-1) - delta(n-2); delta_n1 = 2*delta_n - delta(n-1); delta = [delta_m1 delta_0 delta delta_n delta_n1]; % Akima's derivative estimate formula (equation (1) in the paper): % % H. Akima, "A New Method of Interpolation and Smooth Curve Fitting % Based on Local Procedures", JACM, v. 17-4, p.589-602, 1970. % % s(i) = (|d(i+1)-d(i)| * d(i-1) + |d(i-1)-d(i-2)| * d(i)) % / (|d(i+1)-d(i)| + |d(i-1)-d(i-2)|) % % To eliminate overshoot and undershoot when the data is constant for more % than two consecutive nodes, in MATLAB's 'makima' we modify Akima's % formula by adding an additional averaging term in the weights. % s(i) = ( (|d(i+1)-d(i)| + |d(i+1)+d(i)|/2 ) * d(i-1) + % (|d(i-1)-d(i-2)| + |d(i-1)+d(i-2)|/2) * d(i) ) % / ( (|d(i+1)-d(i)| + |d(i+1)+d(i)|/2 ) + % (|d(i-1)-d(i-2)| + |d(i-1)+d(i-2)|/2) weights = abs(diff(delta)) + abs((delta(1:end-1)+delta(2:end))/2); weights1 = weights(1:n); % |d(i-1)-d(i-2)| weights2 = weights(3:end); % |d(i+1)-d(i)| delta1 = delta(2:n+1); % d(i-1) delta2 = delta(3:n+2); % d(i) weights12 = weights1 + weights2; s = (weights2./weights12) .* delta1 + (weights1./weights12) .* delta2; % If the data is constant for more than four consecutive nodes, then the % denominator is zero and the formula produces an unwanted NaN result. % Replace this NaN with 0. s(weights12 == 0) = 0; end

Copyright 2019 The MathWorks, Inc.

Get
the MATLAB code

Published with MATLAB® R2018b

Why are manhole covers round? It is so they won't fall through the hole they are intended to cover. They have the same diameter regardless of where it is measured. If the hole has a slightly smaller diameter, it is not possible to orient the cover so that it will fall through. A square or rectangular cover can be turned slightly and it will easily fit through the hole.... read more >>

]]>Why are manhole covers round? It is so they won't fall through the hole they are intended to cover. They have the same diameter regardless of where it is measured. If the hole has a slightly smaller diameter, it is not possible to orient the cover so that it will fall through. A square or rectangular cover can be turned slightly and it will easily fit through the hole.

The circle is not the only planar curve that has this constant width property. The Reuleaux Triangle is the next best-known curve of constant width. I want to describe an interactive MATLAB® program that generates generalizations of the Reuleaux Triangle.

*Image Credit: Wikipedia*

Franz Reuleaux (1829-1905) was a German mechanical engineer, a professor at the *Eidgenössische Technische Hochschule Zürich*, then a lecturer at the Berlin Royal Technical Academy and ultimately President of the Academy. In addition to his triangle and many other mechanical creations, he is known for the theory of *kinematic pairs*.

Here is the Reuleaux triangle in action in an animated gif showing one full rotation. The sides of the "triangle" are circular arcs; each is one-sixth of a circle centered at the opposite vertex.

The triangle is rotating, but not about the origin at its center. These plots were made with the MATLAB command

`axis tight square`

so the bounding box is a square with the length of a side equal to the diameter of the curve. The `x` and `y` labels show this diameter. The fact that the diameter does not change as the figure rotates is the defining property of constant width.

These animations vary two parameters, `n` and `delta`. In the next two figures there are `n` blue dots at the vertices of a regular polygonal. There are `n` large arcs; each is one `n`-th of a circle centered at the opposite blue dot. They join the `n` small arcs, each is one n-th of a circle of radius `delta` centered at the nearby blue dot. The resulting curve is smooth, except when `delta` is zero.

The diameter is always the distance from one large circular arc to the opposite small circular arc, no matter how the curve is rotated.

Here `n` is kept constant at three and `delta` varies.

Here `delta` is kept constant at 0.05 and `n` varies over odd integers.

The blue dots do not have to form a regular polygon. Here they form a 3-4-5 triangle.

I will add an interactive app, `constant_width`, to Cleve's Laboratory. You can choose `regular` or `irregular` spacing, vary `n` and `delta`, and rotate the resulting curves at any speed.

I was intrigued by a paper by Stanley Rabinowitz. He describes how curves of constant width can be generated by so-called "support functions." An example is

$$ p(\theta) = a \cos^2{(\frac{k\theta}{2})} + b $$

Parametric equations for the curve are

$$ x = p(\theta) \cos{\theta} - p'(\theta) \sin{\theta} $$

$$ y = p(\theta) \sin{\theta} + p'(\theta) \cos{\theta} $$

Not all choices of the parameters produce convex curves, but $a = 2$, $b = 40$, and $k = 7$ produce a satisfactory heptagonal figure. A constant step size $h = \pi/(12 k)$ yields

Rabinowitz goes on to create a curve of constant width that is actually a *polynomial* of degree eight in the two variables $x$ and $y$.

Several countries mint coins shaped by curves of constant width. They are difficult to counterfeit and work perfectly well in vending machines. Here are a few.

*Image credit: cp4space*

Wankel engines, which incorporate a Reuleaux triangle, were invented in the early 1950s by Felix Wankel, a German engineer.

*Image credit: gizmodo*

Get
the MATLAB code

Published with MATLAB® R2018b

Sedona Arizona is a perfect site for the first MathWorks excursion into lifestyle products.... read more >>

]]>Sedona Arizona is a perfect site for the first MathWorks excursion into lifestyle products.

According to Wikipedia

Sedona is a city ... in the northern Verde Valley region of the U.S. state of Arizona. As of the 2010 census, its population was 10,031. Sedona's main attraction is its array of red sandstone formations. The formations appear to glow in brilliant orange and red when illuminated by the rising or setting sun. The red rocks form a popular backdrop for many activities, ranging from spiritual pursuits to the hundreds of hiking and mountain biking trails.

Sedona is the site of the first MathWorks venture into retail marketing. Similar to the successful Apple outlets, this store features products made with MATLAB that can enhance a visitor's Sedona spiritual experience.

One popular item, intended to attract visitors to the store, is biorhythms. For only $1.99, a customer who reveals their birthday can have a printout of a personal biorhythm. As a special attraction, on a customer's birthday the biorhythm is free.

Here is an example, the biorhythm of someone whose 30th birthday is April 1, 2019. People with this birthday will have an exceptionally fortuitous 10-day period. Two and three days prior to the birthday their physical and emotional capabilities will be at their peak. Their intellectual capability will be at its 33-day maximum a week after the birthday.

biorhythm

For many visitors to the Sedona MATLAB store, the highlight will be their personal energy vortices. For $29.99, a customer who supplies their complete horological background, including the dominant zoological sign at the moment of birth, can obtain a topographic map of the Sedona area with regions of maximum personal spiritual energy highlighted. The service includes the exact longitude and latitude of these personal vortices for use in GPS devices.

Here are examples of energy vortices for certain people.

Another example is a vortex easily reached from the Schuerman Mountain trail.

The view from the Schuerman vortex.

photo credit: Mark Chisholm

Get
the MATLAB code

Published with MATLAB® R2018b

s = s + x(i)*y(i)or "daxpys", double a x plus y,

y(i) = a*x(i) + y(i)In either case, that's one multiplication, one addition, and a handful of indexing, load and store operations. Taken altogether, that's two floating point operations, or two

t = zeros(1,30); m = zeros(1,30); n = 0:500:15000; tspan = 1.0;

while ~get(stop,'value') k = randi(30); nk = n(k); A = randn(nk,nk); b = randn(nk,1); cnt = 0; tok = 0; tic while tok < tspan x = A\b; cnt = cnt + 1; tok = toc; end t(k) = t(k) + tok/cnt; m(k) = m(k) + 1; end

t = t./m; gigaflops = 2/3*n.^3./t/1.e9;This randomly varies the size of the matrix between 500-by-500 and 15,000-by-15,000. If it takes less than one second to compute

plot(n,t,'o')I get The times are increasing like

giga = 2/3*n.^3./t/1.e9I conclude that for

Published with MATLAB® R2018b
]]>

My MATLAB® program, `amaze`, generates mazes by combining old friends, `numgrid` and `delsq`, with a new friend, the `graph` object. Let's see how we make this example maze.... read more >>

My MATLAB® program, `amaze`, generates mazes by combining old friends, `numgrid` and `delsq`, with a new friend, the `graph` object. Let's see how we make this example maze.

Select a maze size. Our example is 15-by-15.

n = 15;

We are going to make two graphs, the barriers and the cells. The barriers, `B`, is the graph of the classic five-point discrete Laplacian on an n-by-n square grid. Each node is the interior is connected to its four closest neighbors and each node on the boundary has two or three neighbors. The functions `delsq` and `numgrid` have been in the `sparfun` directory since sparse matrices were introduced in 1992. The `graph` object was introduced in R2015b. We are going to carve a path through this grid.

Db = delsq(numgrid('S',n+2)); B = graph(logical(Db),'omitselfloops');

The cells, `C`, are also based on the five-point discrete Laplacian, but on an (n-1)-by-(n-1) grid. Initially `C` has just the nodes, no edges. When we plot `B` and `C` together, the nodes of `C` are centered in the squares of `B`.

m = n-1; Dc = delsq(numgrid('S',m+2)); C = graph(logical(Dc),'omitselfloops');

available = 1:m^2; % Nodes we haven't visited yet. branches = []; branching = 'middle'; tree = zeros(0,2); % Depth first search. p = 1;

Now we make a series of random walks on `C`.

while true % Break when available is empty

available(p) = 0; if ~any(available) break end

Each step is from a node to a random choice of one of its neighbors that has not already been visited.

```
[~,~,ngh] = find(available(neighbors(C,p)));
if ~isempty(ngh)
```

idx = randi(length(ngh)); % Random choice. q = ngh(idx); % Next cell. tree(end+1,:) = [p q];

At the same time, the edge in the barrier `B` that is blocking the step is removed.

[i,j] = wall(p,q); B = rmedge(B,i,j);

We keep a list of the nodes where a random choice between two or three neighbors is made.

if length(ngh) > 1 branches(end+1) = p; end p = q;

Eventually the walk reaches a dead end, a node whose neighbors have all been visited. Here is the first such walk in our example.

When we reach a dead end, we backtrack to one of the nodes on our list of nodes where we had a choice of available neighbors. Here another kind of choice is made. We can go back to the *first* node on the list, to the node in the *middle* of the list, or to the *last* node on the list.

else switch branching case 'first' idx = 1; case 'last' idx = length(branches); otherwise idx = round(length(branches)/2); end p = branches(idx); branches(idx) = []; end

Here is the result of backtracking to the middle and then walking to a second dead end.

Here is the result after several more backtracks.

```
end
C = graph(tree(:,1),tree(:,2));
```

Near the end of this process, the length of a path before a dead end and a backtrack becomes shorter and shorter. A path may involve a single point. Eventually all the nodes in `C` are covered. Here is the final result.

When the entire path is plotted with `linestyle = 'none'` we have the maze.

We can easily "solve" the maze by asking for the shortest path from one corner to another.

path = shortestpath(C,1,m^2);

(*Exercise*: How does the `branching` parameter affect the average length of this `path`?)

I have a program called `amazing` that allows you to investigate `amaze`. I've submitted it to the MATLAB Central File Exchange and will put the link here. I will also add `amazing` to version 4.40 of Cleve's Laboratory.

Thanks to Pat Quillen who provided invaluable help.

Get
the MATLAB code

Published with MATLAB® R2018b