The MATLAB functions for the numerical evaluation of integrals has evolved from `quad`, through `quadl` and `quadgk`, to today's `integral`.... read more >>

The MATLAB functions for the numerical evaluation of integrals has evolved from `quad`, through `quadl` and `quadgk`, to today's `integral`.

Several years ago I regarded myself as an expert on numerical methods for approximating integrals. But then I stopped paying attention for a while. So this is an opportunity to get reacquainted. The task involves a real-valued function $f(x)$ of a real variable $x$, defined on an interval $a \le x \le b$. We seek to compute the value of the definite integral,

$$ \int_a^b \ f(x) \ dx $$

I wrote a chapter on the subject in my textbook *Numerical Computing with MATLAB*. Here is a link to that chapter. The chapter is titled "Quadrature". One of the first things I have to do in that chapter is explain that *quadrature* is an old-fashioned term derived from the notion that the value of a definite integral can be approximated by plotting the function on graph paper and then counting the number of little squares that fall underneath the graph. Here is a coarse picture.

For most of the time since the beginning of MATLAB, up until when I wrote the chapter on quadrature twelve years ago, the MATLAB function for evaluating definite integrals was one named `quad`. It was based on adaptive Simpson's method, pioneered by my grad school buddy Bill McKeeman to show off recursion in Algol.

These two figures illustrate the basis for adaptive Simpson's method.

The algorithm works recursively on subintervals of the original interval [a,b]. If the values obtained from the 3-point Simpson's rule and the 5-point composite Simpson's rule agree to within a specified tolerance (which they clearly do not for these figures), then an extrapolation taken from the two values is accepted as the value for the integral over the subinterval. If not, the algorithm is recursively applied to the two halves of the subinterval.

I have used a function named `humps` to demonstrate numerical integration for a long time. It has a strong peak at x = 0.3 and a milder peak at x = 0.9. It has found its way into the MATLAB `demos` directory, so you have it on your own machine. A graph appears below.

$$ h(x) = \frac{1}{(x-0.3)^2+.01} + \frac{1}{(x-0.9)^2+.04} - 6 $$

An interactive graphical app named `quadgui` is included with the programs for *Numerical Computing with MATLAB*, available from MATLAB Central File Exchange. It shows `quad` in action. Here is a snapshot of `quad` part way through computing the integral of `humps` over the interval [0,1]. Integration over the interval to the left of 1/4 is finished. The integration over the subinterval [1/4,3/8] has just been successfully completed. A few samples to the right of 3/8 were computed during the initialization. The integrand has been evaluated 19 times.

Here is the final result. The default tolerance for `quad` is $10^{-6}$ and the value of the integral obtained to six significant digits is 29.8583. The integrand has been evaluated 93 times. Most of those evaluations took place near the peaks where the `humps` varies rapidly. This indicates that the adaptivity works as expected.

`Humps` is a rational function. It is possible to use the Symbolic Toolbox to obtain the indefinite and then the definite integral analytically.

syms x h = 1/((x-3/10)^2+1/100) + 1/((x-9/10)^2+4/100) - 6 I = int(h,x) D = subs(I,1)-subs(I,0) format long Q = double(D)

h = 1/((x - 9/10)^2 + 1/25) + 1/((x - 3/10)^2 + 1/100) - 6 I = 10*atan(10*x - 3) - 6*x + 5*atan(5*x - 9/2) D = 5*atan(1/2) + 10*atan(3) + 10*atan(7) + 5*atan(9/2) - 6 Q = 29.858325395498674

This confirms the six figure value shown by `quadqui`

About the time in 2004 that I was writing the *Quadrature* chapter for *NCM*, the `quadl` function showed up. This uses the *Lobatto* rule.

$$ \int_{-1}^1 \ f(x) \ dx \ \approx \ w_1 f(-1) + w_2 f(-x_1) + w_2 f(x_1) + w_1 f(1) $$

with $x_1 = 1/\sqrt{5}, \ w_1 = 1/6, \ w_2 = 5/6$.

`Quadl` also employs the associated *Kronrod* rule, which provides 13 additional irrational abscissa and weights. The Lobatto-Kronrod rule is much higher order and consequently more accurate than Simpson's rule. As a result `quadl` is usually faster and more accurate than `quad`. But I just devoted a couple of exercises in the book to `quadl`. And I think it never became very popular.

In the mathematical software business when we were developing libraries using Fortran or C, we evaluated the cost of numerical quadrature by counting the number of function evaluations. In 2006 Larry Shampine, the SMU Professor who is a consultant to MathWorks and principal author of our ODE Suite, published the paper referenced below where he changed the ground rules for quadrature functions in the MATLAB environment. Since MATLAB function evaluations are vectorized, we should count the number of vectorized evaluations, not the number of individual evaluations. This is a better predictor of execution time.

Shampine's paper includes a MATLAB function that he calls `quadva`, where `va` stands for "vectorized adaptive". His code uses a 15-point Gauss-Kronrod formula. The Gauss formula is not only higher order than the Lobatto formula with the same number of points, it has the additional advantage that it does not require evaluating the integrand at the end points of an interval. This makes it possible to treat infinite intervals and singularites at the end points.

When we added Shampine's `quadva` to MATLAB, we called it `quadgk`, for "Gauss-Kronrod".

Let's instrument `humps` to investigate the behavior of our quadrature routines on just one simple test case. The functions `quad` and `quadl` call the integrand once with a vector argument to initialize the integration and then call it more times with a shorter vector argument during the adaptive process. This is such an easy problem that `quadgk` gets all the information that it needs during the initialization.

In the following table, "init" is the length of the vector in the initialization call, "kount" is the number of times that `humps` is called after the initialization, "length" is the length of the vector argument in those calls, "err" is the error in the final result, and "time" is the execution time in milliseconds.

$$\begin{array}{lrrrrr} & \mbox{init} & \mbox{kount} & \mbox{length} & \mbox{err\ \ \ \ } & \mbox{time}\\ \mbox{quad} & 7 & 69 & 2 & 7.3 \cdot 10^{-7} & 0.87 \\ \mbox{quadl} & 13 & 37 & 5 & 1.8 \cdot 10^{-10} & 0.78 \\ \mbox{quadgk} & 150 & 0 & - & 2.5 \cdot 10^{-14} & 0.34 \end{array}$$

We see that the three functions perform very differently in this one experiment. The oldest function, `quad`, actually requires the fewest total function evaluations, if we count them in the traditional scalar sense, but it is the least accurate and the slowest. The newest function, `quadgk`, requires only one vectorized function evaluation, and is the fastest and the most accurate.

In 2012, with release 2012a, we came to the realization that most new MATLAB users who wanted to evaluate an integral had a hard time finding the function named "quadgk" that does the job. For all these years, we numerical analysts have been explaining to the rest of the world our quaint term *quadrature*. If that weren't enough, we have to tack on initials of mathematicians. So, we introduced `integral`. Actually `quadgk` is still available. `Integral` is just an easier to find and easier to use version of `quadgk`.

L.F. Shampine, "Vectorized Adaptive Quadrature in MATLAB, *Journal of Computational and Applied Mathematics* 211 (2006), pp.131–140, available at <http://faculty.smu.edu/shampine/rev2vadapt.pdf>

Get
the MATLAB code

Published with MATLAB® R2016a

Gil Strang has produced a MOOC-style video course on Differential Equations and Linear Algebra. I have added some videos about the MATLAB ODE suite. The series is available from the MathWorks Web site, MIT OpenCourseWare and several other popular sources.... read more >>

]]>Gil Strang has produced a MOOC-style video course on Differential Equations and Linear Algebra. I have added some videos about the MATLAB ODE suite. The series is available from the MathWorks Web site, MIT OpenCourseWare and several other popular sources.

MIT's Professor Gilbert Strang has taught 18.06, Linear Algebra, regularly for many years. Over fifteen years ago MIT videotaped the lectures. MIT was also taping 8.01, the first-year physics course, taught by the flamboyant Walter Lewin, in the same lecture room the hour before. They were careful with the videography, with professional camera work, good lighting and good audio. There was no computer graphics, no PowerPoint. Just excellent lectures, using big chalk on actual black boards. You can still see Gil's original lectures by Googling "18.06", or by following this link.

MIT decided to make the lectures freely available on the web, with no strings attached. That was an unusual move at the time. The two courses proved to be very popular. They led to the creation of what MIT now calls OpenCourseWare and were in part responsible for the MOOC, Massive Open Online Courses, revolution in higher education.

Gil told me the experience changed his life. He realized he could reach millions of students instead of just one classroom.

Strang recently published a new book, Differential Equations and Linear Algebra.

He invited MathWorks to join him in producing a series of videos on the subject. I joined the project to describe how MATLAB provides tools for solving ODEs.

Here is a link to the home page for the project at the MathWorks. Learn Differential Equations.

And here is a link to a video where Gil and I have a short chat. Gil and Cleve chat.

Gil produced 55 videos, averaging about 15 minutes each. There are eight topic areas, which map to and complement the textbook, DELA.

- Introduction
- First Order Equations
- Second Order Equations
- Graphical and Numerical Methods
- Vector Spaces and Subspaces
- Eigenvalues and Eigenvectors
- Applied Mathematics and ATA
- Fourier and Laplace Transforms

Here is a screen shot from one of the videos to give you a feel for the style of the presentation.

My part of the series is much smaller, only a dozen videos, also about 15 minutes each. The first video is about Euler's method. My Euler's method code is `ode1.m`. It is 16 lines long; half of the lines are comments. The last three videos are examples. One of them is based on my blog post about the tumbling box ode. Gil also discusses the tumbling box in one of his videos. I have included several hands on exercises that you can try out by downloading the MATLAB code available on the MATLAB Central File Exchange, solving ODEs in MATLAB.

Here is a screen shot from my discussion of stiffness.

The home page for the project at the MathWorks is <http://www.mathworks.com/academia/courseware/learn-differential-equations.html>.

The series is available from MIT OpenCourseWare under the title *Learn Differential Equations: Up Close with Gilbert Strang and Cleve Moler*. The OCW course number is 18-009. The URL is <http://ocw.mit.edu/resources/res-18-009-learn-differential-equations-up-close-with-gilbert-strang-and-cleve-moler-fall-2015>.

The videos are also available through the following distribution sites.

YouTube: <http://www.youtube.com/playlist?list=PLUl4u3cNGP63oTpyxCMLKt_JmB0WtSZfG>.

iTunes: <http://itunes.apple.com/us/itunes-u/id1102923623>.

Internet Archive: <http://archive.org/details/MITRES18-009F15>.

Get
the MATLAB code

Published with MATLAB® R2016a

The equations generating a `surf` plot of the Moebius strip can be parameterized and the parameters allowed to take on expanded values. The results are a family of surfaces that I have been displaying for as long as I have had computer graphics available.... read more >>

The equations generating a `surf` plot of the Moebius strip can be parameterized and the parameters allowed to take on expanded values. The results are a family of surfaces that I have been displaying for as long as I have had computer graphics available.

The animation shows the formation of a classic Moebius strip. Start with a long, thin rectangular strip of material. Bend it into a cylinder. Give the two ends a half twist. Then join them together. The result is a surface with only one side and one boundary. One could traverse the entire length of the surface and return to the starting point without ever crossing the edge.

I am going to investigate the effect of three free parameters:

- $c$: curvature of the cylinder.
- $w$: width of the strip.
- $k$: number of half twists.

These parameters appear in a mapping of variables $s$ and $t$ from the square $-1 <= s,t <= 1$ to a generalized Moebius strip in 3D.

$$ r = (1-w) + w \ s \ \sin{(k \pi t/2)} $$

$$ x = r \ \sin{(c \pi t)}/c $$

$$ y = r \ (1 - (1-\cos{(c\pi t))}/c) $$

$$ z = w \ s \ \cos{(k \pi t/2)} $$

Here is the code. The singularity at $c = 0$ can be avoided by using `c = eps`.

```
type moebius
```

function [x,y,z,t] = moebius(c,w,k) % [x,y,z,t] = moebius(c,w,k) % [x,y,z] = surface of moebius strip, use t for color % c = curvature, 0 = flat, 1 = cylinder. % w = width of strip % k = number of half twists if c == 0 c = eps; end m = 8; n = 128; [s,t] = meshgrid(-1:2/m:1, -1:2/n:1); r = (1-w) + w*s.*sin(k/2*pi*t); x = r.*sin(c*pi*t)/c; y = r.*(1 - (1-cos(c*pi*t))/c); z = w*s.*cos(k/2*pi*t); end

The animation begins with a flat, untwisted strip of width 1/4, that is

$$ c = 0, \ w = 1/4, \ k = 0 $$

The classic Moebius strip is reached with

$$ c = 1, \ w = 1/4, \ k = 1 $$

The final portion of the animation simply changes the viewpoint, not the parameters.

Let's play with colormaps. I'll look at the classic Moebius strip and make the default MATLAB colormap periodic by appending a reversed copy. I think this looks pretty nice.

map = [parula(256); flipud(parula(256))];

The Hue-Saturation-Value color map has received a lot of criticism in MathWorks blogs recently, but it works nicely in this situation. This figure has $k = 4$, so there are two full twists. That's very hard to see from this view. The colors help, but it's important to be able to rotate the view, which we can't do easily in this blog.

Here is the one with three half twists. Like the classic Moebius strip, this has only one side. I've grown quite fond of this guy, so I'll show two views. The HSV color map, with its six colors, works well.

Source: IEEE

I have written in this blog about the glorious marriage between MATLAB and the Dore graphics system on the ill-fated Ardent Titan computer almost 30 years ago. I somehow generated this image using MATLAB, Dore, these parametrized Moebius equations and the Titan graphics hardware in 1988. This was fancy stuff back then. The editors of IEEE Computer were really excited about this cover.

I've never been able to reproduce it. I don't know what the parameters were.

Kermit Sigmon was a professor at the University of Florida. He was an excellent teacher and writer. He wrote a short *MATLAB Primer* that was enormously popular and widely distributed in the early days of the web. It went through several editions and was translated into a couple of other languages.

For the most part, people had free copies of the *MATLAB Primer* without fancy covers. But CRC Press obtained the rights to produce a reasonably priced bound paperback copy which featured a cover with $k = 5$.

Kermit passed away in 1997 and Tim Davis took over the job of editing the *Primer*. He generated more elaborate graphics for his covers.

Get
the MATLAB code

Published with MATLAB® R2016a

A model of the human gait, developed by Nikolaus Troje, is a five-term Fourier series with vector-valued coefficients that are the principal components for data obtained in motion capture experiments involving subjects walking on a treadmill.... read more >>

]]>A model of the human gait, developed by Nikolaus Troje, is a five-term Fourier series with vector-valued coefficients that are the principal components for data obtained in motion capture experiments involving subjects walking on a treadmill.

Today's post is about work done several years ago by Nikolaus Troje when he was at Ruhr-Universität in Bochum, Germany. He has since moved to Queens University in Ontario, Canada, where he is head of the Biomotion Lab. Here is their Web page.

Troje studies the human gait. How exactly do we move when we walk? He also studies how we perceive others when they walk. Can we tell if people are happy or sad just by the way they walk? Can we determine if they are male or female?

To get a feeling for his work, take a look at this demo: BMLgender. When I first saw this demo it was being done with MATLAB, so I asked Troje about it and I have written and talked about it many times ever since. (Today, the demo is computed in the web browser with Java Script.)

Motion capture employs high speed cameras, lights, and reflectors attached to human subjects -- German grad students -- as they walk on a closed track or a treadmill. Once a steady, comfortable pace is established for a single subject, 20 or so steps are recorded. At 120 frames per second, this results in typically 1440 frames. The data from the reflectors are collated into 15 standard locations on the head, arms, torso and legs. As these points move in three dimensions the raw data captured is an array of time series with 45 rows and roughly 1440 columns.

Here are four of the subjects, with the data scaled so that they are all walking at the same speed.

The walking motion is periodic so Troje's model is a five-term Fourier series describing the position of the 15 spots in three dimensions.

$$ p(t) = v_1 + v_2 \cos{\omega t} + v_3 \sin{\omega t} + v_4 \cos{2 \omega t} + v_5 \sin{2 \omega t} $$

Here $\omega$, the scalar frequency, is the only nonlinear parameter. It is determined by Fourier analysis of the data. The coefficients $v_1$ through $v_5$ are *vectors* with 45 components. Once $\omega$ is computed the $v$'s are obtained by *Principal Component Analysis* (also known as *Singular Value Decomposition*) of the data. PCA effectively reduces the 45-by-1440 data matrix to a 45-by-5 coefficient matrix.

The frequency $\omega$ and the coefficient matrix $V$ with columns $v_1, \ldots, v_5$ are computed just once. Then the motion can be computed for any time $t$ with a matrix-vector product.

$$ p(t) = V \ [1 \ \cos{\omega t} \ \sin{\omega t} \ \cos{2 \omega t} \ \sin{2 \omega t}]^T $$

This model usually reproduces the raw data with better than 95% accuracy.

There were several dozen subjects, each with his or her characteristic coefficient matrix $V$. Collect all these $V$'s together into one giant set of coefficients. Then use PCA again to reduce this to one average set of coefficients. You get *eigenwalker*, the *everyman*. Or at least the *typisch deutsche junge Akademiker*.

The *eigenwalker* has five vector coefficients $v_1, \ldots, v_5$. These components don't have any inherent names -- they're just vectors that come out of the principal component analysis. But when we look at how they affect the motion, we can give them names. The first component, $v_1$, is simply the static position. It determines how *eigenwalker* looks at rest. Here are the other four components, one at a time.

The second coefficient, $v_2$, dominants the forward motion so we've called it "stride". The third coefficient, $v_3$, controls a sideways motion so it gets named "sway". The coefficients $v_4$ and $v_5$ drive vertical motions that are in-phase and out-of-phase with the stride; they get fairly arbitrary names of "hop" and "bounce".

When these four motions are added to the at-rest position, they give *eigenwalker*, which is a fairly realistic picture of the human gait.

The PCA does not pick out gender as one of dominant components. But half of the subjects were female and half were male. So we can compute a matrix of coefficients $F$ for the females and a matrix of coefficients $M$ for the males. The top two figures below are the female *eigenwalker* and the male *eigenwalker*. Without looking at the bottom two figures, can you tell which is which?

The bottom two figures overemphasize the gender bias. They run the model with the difference between $F$ and $M$ tripled. Now you can see which is which?

Disney and Pixar know this -- they have to exaggerate the motion to get the point across.

This is Troje's research. How do humans determine the gender of other humans just by watching them walk? Even better -- just by watching a stick figure, or 15 dots, walk?

The *Numerical Computing with MATLAB* collection of programs includes an app that lets you see the *eigenwalker* in action. Sliders vary the strength of the five components and the gender bias. You can also rotate the viewpoint. You will need two files, `walker.m` and `walkers.mat`, available from the MATLAB Central File Exchange.

This has been my own description of an approach that Troje took over a decade ago. He has since replaced Principle Components Analysis with Fourier analysis to represent the movements of an individual walker, and has applied PCA to all sorts of other walker properties on the population level. But I like his original approach because I am such a big fan of the SVD. The following references are mostly about Troje's more recent work.

<http://jov.arvojournals.org/article.aspx?articleid=2192503>

<http://www.biomotionlab.ca/walking.php>

<http://www.biomotionlab.ca/Text/WDP2002_Troje.pdf>

<http://www.biomotionlab.ca/Text/trojeChapter07b.pdf>

Get
the MATLAB code

Published with MATLAB® R2016a

Recent theoretical, observational and computational results establish the possibility that gravitational waves produced by the dark energy created at the dawn of the universe affect the clock rate of silicon digital processors operating at very low temperatures.... read more >>

]]>Recent theoretical, observational and computational results establish the possibility that gravitational waves produced by the dark energy created at the dawn of the universe affect the clock rate of silicon digital processors operating at very low temperatures.

Ed Plum is a professor in the Institute for Theoretical Physics at the U. C. Santa Barbara. He is also a long-time friend of mine and an avid MATLAB fan.

Ed and his colleagues are publishing an important paper today. Here is a link to the journal, *Physical Review Communications*. Ed has shared an advanced copy with me. They have found convincing evidence that the gravitational waves predicted by Einstein's General Theory of Relativity and emanating from "Dark Energy" affect the operating frequency of silicon computer chips. The effect can be observed in experiments operating at temperatures below 2 degrees Kelvin.

The announcement on February 11 of the detection of gravitational waves by the LIGO team thrilled scientists and the public at large worldwide. LIGO stands for Laser Interferometer Gravitational-Wave Observatory. Here is the link to the LIGO home page. Here is their scientific paper in *Physical Review Letters*. Excellent survey articles were published in *The New York Times Sunday Review* and *The New Yorker*.

Here is a Wikipedia animated graphic depicting gravitational waves in the article on General Relativity.

The LIGO experiment is designed to detect powerful one-time astrophysical events such as colliding black holes. The gravitational waves involved have wave lengths that require L-shaped detectors with legs four kilometers long. A wave passes through the detector just once, in a fraction of a second. There are two detectors, located in Washington State and Louisiana.

Photo: <http://www.ligo.caltech.edu >

This beautiful graphic from the *Phys. Rev. Letters* paper shows the event detected by the LIGO team. The top two figures show the actual signals recorded at the two labs. The next two figures show the results from supercomputer runs for the numerical solution of Einstein's equations. The small pair of figures show that the difference between the observations and the simulations is noise. The two color figures show the "chirps", spectrograms plotting frequency versus time. These chirps are readily audible to the human ear.

Source: <http://www.ligo.caltech.edu>

The most commonly accepted model of cosmology hypothesizes *dark energy*, formed at the origin of the universe and permeating all of space, but invisible to optical observation. According to Einstein's General Theory of Relativity, this energy constantly produces short wave length, very low energy gravitational waves. These waves continuously interact with all mass in the universe, but on a scale that is almost impossible to detect.

Plum and his colleagues have evidence from both numerical simulations and actual observations that the dark energy gravitational waves affect the operating frequency of commercial microprocessors. Very slight changes in processor speed can be detected, provided the silicon is exceptionally pure and the observations are made at very low temperatures. It's like immersing a high-end laptop in liquid nitrogen. Here is one of the detectors.

Photo: Plum *et al*

Some time ago, Ed gave me access to one of the first signals his group had detected and we have been quietly including it as one of sample sounds on recent MATLAB distributions. The file is `chirp.mat`. We have reversed the signal to disguise it until now. You should be able to load a copy on your own machine.

clear which chirp.mat load chirp.mat y = flipud(y)'; t = (1:length(y))/Fs; whos

C:\Program Files\MATLAB\R2016a\toolbox\matlab\audiovideo\chirp.mat Name Size Bytes Class Attributes Fs 1x1 8 double t 1x13129 105032 double y 1x13129 105032 double

The first plot is the entire signal. Since the waves are arriving continuously, we see several peaks. The second plot zooms in on one of the peaks.

```
subplot(2,1,1)
plot(t,y)
axis([0 t(end) -1 1])
subplot(2,1,2)
plot(t,y)
axis([.75 .80 -1 1])
xlabel('Time (secs)')
```

I hope you can run this `sound` yourself and hear the chirps. If you are just reading this on the Web, you will have to imagine what this sounds like.

sound(y,Fs)

Here is the picture of one chirp, the spectrogram of about one-eighth of the signal.

clf reset spectrogram(y(1:2048),kaiser(128,18),120,128,Fs,'yaxis')

This is just be beginning. We expect to learn much more about our universe as we are able to refine the silicon in the microprocessors and lower the temperature of the operating environment.

Get
the MATLAB code

Published with MATLAB® R2016a

An extraordinarily creative Danish mathematician, inventor, and poet who often wrote under the Old Norse pseudonym "Kumbel" meaning "tombstone." A direct descendant of the Dutch naval hero of the 16th century who had the same name, Piet Hein was born in Copenhagen and studied at the Institute for Theoretical Physics of the University of Copenhagen (later the Niels Bohr Institute) and the Technical University of Denmark. ... read more >>

]]>

Let's quote the entry about Piet Hein in David Darling's *Enclopedia of Science*, http://www.daviddarling.info/encyclopedia/H/Hein.html

An extraordinarily creative Danish mathematician, inventor, and poet who often wrote under the Old Norse pseudonym "Kumbel" meaning "tombstone." A direct descendant of the Dutch naval hero of the 16th century who had the same name, Piet Hein was born in Copenhagen and studied at the Institute for Theoretical Physics of the University of Copenhagen (later the Niels Bohr Institute) and the Technical University of Denmark.

[Piet Hein] is famed for his many mathematical games, including Hex, Tangloids, Polytaire, TacTix, and the Soma Cube, his advocacy of the superellipse curve in applications as diverse as city planning and furniture making, and his thousands of short, aphoristic poems called Grooks.

This is one line from one of the Grooks.

"Problems worthy of attack prove their worth by hitting back!" -- Piet Hein

My wife and I purchased one of these Piet Hein pocelain baking dishes years ago during a visit to Denmark. It's been a mainstay in our kitchen ever since. This photo is not from our kitchen, it's from a page at http://www.piethein.com featuring http://www.piethein.com/category/dishes-white-22.

Hein's design of the baking dish, as well as many other items, is based on what he called the *superellipse*, a curve that lies between old-fashioned ellipses and rectangles. The important parameter $p$ controls the roundness. Additional parameters $a$ and $b$ determine the aspect ratio. A superellipse with these parameters is the set of points $(x,y)$ satisfing

$$ |\frac{x}{a}|^p + |\frac{y}{b}|^p = 1 $$

If $a$ and $b$ are both equal to one, the equation is

$$ |x|^p + |y|^p = 1 $$

And if $p$ is equal to two it becomes

$$ x^2 + y^2= 1 $$

which defines a circle.

Here is the MATLAB definition.

F = @(x,y) abs(x/a).^p + abs(y/b).^p - 1;

Let's vary `p` between 0.1 and 10 and make an animated gif, pausing briefly as we pass through `p = 1` and `p = 2`.

```
type make_superellipse_movie
```

gif_frame('superellipse_movie.gif') for p = [0.1 0.1 0.1:0.1:1 1 1 1:0.2:2 2 2 2:.5:10] F = @(x,y) abs(x).^p + abs(y).^p - 1; ezplot(F,[-1.2 1.2]) axis square title(['|x|^p + |y|^p = 1, p = ',sprintf('%.1f',p)]) gif_frame end

All of Piet's designs have `p = 2.5`.

p = 2.5; a = 1.5; F = @(x,y) abs(x/a).^p + abs(y).^p - 1; ezplot(F,[-1.8 1.8]) axis equal title(sprintf('Superellipse, p = %.1f',p))

Different scale factors on `x` and `y` make the shape rectangular and `p` greater than two makes it not exactly an ellipse.

I could have said that the superellipse is the unit ball in two dimensions of the vector p-norm,

$$ ||v||_p = \big( \sum_i |v_i]^p \big)^{1/p} $$

In MATLAB this is `norm(v,p)`. In other words the function `F` starring in the first animated gif could have been defined by

F = @(x,y) norm([x,y],p) - 1;

When `p = 1` the unit "ball" becomes the unit diamond. And for `p < 1` the quantity `norm(v,p)` is no longer actually a norm. It fails to satisfy the triangle inequality.

Homework: what is `norm(v,p)` for *negative* `p`, and when `p = -Inf`?

Use the Images option in Google to search for all of the different images of "Soma Cube" available on the Internet. It's awesome.

Here is one of those images, available from http://mathcats.org/ideabank/geometrykits.html. This is all possible shapes with a reentrant corner that can be made out of three or four "cubelets". There are a total of 27 of those cubelets. The goal is to make a 3-by-3-by-3 larger cube.

My father's hobby was woodworking. He had a very nice shop. One year when I was a kid, he and I made 12 Soma Cube puzzle sets like this one to give to our family for Christmas. My set is now over 65 years old, but I still have it in my office at home.

The following is the shortest code segment I've ever presented in over 2-1/2 years of Cleve's Corner blog. In the command window of your own copy of MATLAB you can enter.

soma

This launches an old demo written by Bill McKeeman. Click away, to your heart's content. In this animated gif, watch the solution number change . It's looping over only the first 12 out of the 240 possible puzzle solutions that McKeeman's algorithm finds.

McKeeman and I were buddies in grad school. We spent too many hours perfecting our Some Cube skills. In fact, I can still put the puzzle together behind my back.

Get
the MATLAB code

Published with MATLAB® R2015b

Today's blog post is a complete working MATLAB program investigating the crossed ladders problem. Download a copy of the program via the link at the end. Publish it again with the `publish` command or the `publish` editor tab.... read more >>

Today’s blog post is a complete working MATLAB program investigating the crossed ladders problem. Download a copy of the program via the link at the end. Publish it again with the `publish` command or the `publish` editor tab.

Explore the interaction of eight variables satisfying five nonlinear equations.

The eight variables are:

- $a$, length of one ladder.
- $b$, length of the other ladder.
- $c$, height of the point where they cross.
- $u$, width of one base triangle.
- $v$, width of the other base triangle.
- $w$, width of the alley.
- $x$, height of the point where one ladder meets the wall.
- $y$, height of the point where the other ladder meets the other wall.

The five equations are:

$$ a^2 = x^2 + w^2 $$

$$ b^2 = y^2 + w^2 $$

$$ \frac{x}{w} = \frac{c}{v} $$

$$ \frac{y}{w} = \frac{c}{u} $$

$$ w = u + v $$

The initial values provide the smallest solution that is all integers. Are there any other integer solutions within the range of this app?

```
function ladders_app
```

% The Classic Crossed Ladders Puzzle. % See Cleve's Corner, Feb. 29, 2016. % <http://blogs.mathworks.com/cleve/2016/02/29/the-classic-crossed-ladders-puzzle> % Starting values are a minimal sum integer solution. w0 = 56; x0 = 105; y0 = 42; w = w0; x = x0; y = y0; c = x*y/(x+y); u = c*w/y; v = c*w/x; % Starting toggles have both labels off. letters = false; numbers = false; initialize_figure;

Snapshot with `letters = true` and a wider valley.

Snapshot of the mirror image of the previous values with `numbers = true`. Four of the eight values are still integers.

function initialize_figure % Size of the square view h = 160; clf shg set(gcf,'menubar','none', ... 'numbertitle','off', ... 'name','Crossed Ladders', ... 'windowbuttondownfcn',@down, ... 'windowbuttonupfcn',@up) axis([-h/8 7/8*h -h/8 7/8*h]) axis square lines buttons labels end %initialize_figure

Called repeatedly as the mouse moves.

function motion(varargin) % WindowsButtonMotionFunction [p,q] = read_mouse; u = c*w/y; v = c*w/x; % Which control point? if q < c/2 % Drag the alley width horizonally. w = p; elseif p < u/2 % Drag the ladder against the left wall vertically. x = q; elseif p > w-v/2 % Drag both the width and the other ladder. w = p; y = q; else % Drag the crossing point, creating lots of action. c = q; u = p; v = w-u; x = c*w/v; y = c*w/u; end c = x*y/(x+y); lines end % motion

Play dot-to-dot.

function lines cla ms = 3*get(gcf,'defaultlinemarkersize'); % Large dots u = c*w/y; % Five markers at the vertices. line([0 0],[0 0],'Marker','.','MarkerSize',ms) line([w w],[0 0],'Marker','.','MarkerSize',ms,'Color','k') line([0 0],[x x],'Marker','.','MarkerSize',ms,'Color','k') line([w w],[y y],'Marker','.','MarkerSize',ms,'Color','k') line([u u],[c c],'Marker','.','MarkerSize',ms,'Color','k') % Connect the markers with six lines. line([0 w],[0 0]) line([0 0],[0 x]) line([w w],[0 y]) line([0 w],[x 0]) line([w 0],[y 0]) line([u u],[c 0],'LineStyle','-.') box on end % lines

function down(varargin) % Called at the start of mouse movement. % Activate the motion function. cla set(gcf,'windowbuttonmotionfcn',@motion) end % down

function up(varargin) % Called at the end of mouse movement. % Deactivate motion function. set(gcf,'windowbuttonmotionfcn',[]) set(gcf,'windowbuttondownfcn',@down) labels end % up

function buttons % Two toggles and three pushbuttons. posit = [.88 .64 .10 .06]; delta = [ 0 .08 0 0]; % letters uicontrol('style','toggle', ... 'units','normalized', ... 'position',posit, ... 'string','letters', ... 'value',letters, ... 'callback',@letters_cb); % numbers uicontrol('style','toggle', ... 'units','normalized', ... 'position',posit-delta, ... 'string','numbers', ... 'value',numbers, ... 'callback',@numbers_cb); % mirror uicontrol('style','pushbutton', ... 'units','normalized', ... 'position',posit-2*delta, ... 'string','mirrror', ... 'callback',@mirror_cb); % reset uicontrol('style','pushbutton', ... 'units','normalized', ... 'position',posit-3*delta, ... 'string','reset', ... 'callback',@reset_cb); % close uicontrol('style','pushbutton', ... 'units','normalized', ... 'position',posit-4*delta, ... 'string','close', ... 'callback','close(gcf)'); end % buttons

function labels % Label lines with either letters or numbers. u = c*w/y; v = c*w/x; lines if letters f = '%s'; xs = 'x '; ys = ' y'; ws = ' w '; as = ' a'; bs = 'b '; us = ' u '; vs = ' v '; cs = 'c '; top = 'Drag a black dot.'; elseif numbers % This is the only place values of a and b are needed. a = sqrt(x^2 + w^2); b = sqrt(y^2 + w^2); z = [a b c u v w x y]; if all(abs(z-round(z)) < .001) f = '%.0f'; else f = '%.1f'; end if sum(z) < 50 || sum(z) > 5000 scream end xs = x; ys = y; ws = w; as = a; bs = b; us = u; vs = v; cs = c; top = ['sum = ' sprintf(f,sum(z))]; else top = 'Drag a black dot.'; end if letters || numbers text(0,x/2,sprintf(f,xs),'horiz','right') text(w,y/2,sprintf(f,ys),'horiz','left') text(w/2,-3,sprintf(f,ws),'horiz','center') text(w/2,x/2,sprintf(f,as),'horiz','left') text(w/2,y/2,sprintf(f,bs),'horiz','right') text(u/2,3,sprintf(f,us),'horiz','center') text(w-v/2,3,sprintf(f,vs),'horiz','center') text(u,c/2,sprintf(f,cs),'horiz','right') end title(top) end % labels

function letters_cb(varargin) % Called when letters button is toggled. letters = get(findobj(gcf,'string','letters'),'value'); % Make sure numbers is false. numbers = false; set(findobj(gcf,'string','numbers'),'value',numbers) labels end % letter_cb

function numbers_cb(varargin) % Called when numbers button is toggled. numbers = get(findobj(gcf,'string','numbers'),'value'); % Make sure letters is false. letters = false; set(findobj(gcf,'string','letters'),'value',letters) labels end % numbers_cb

function mirror_cb(varargin) % Called when mirror button is pushed. % Interchange x and y. t = x; x = y; y = t; % c = x*y/(x+y) is unchanged lines labels end % mirror_cb

function reset_cb(varargin) % Called when reset button is pushed. % Restore initial values. w = w0; x = x0; y = y0; c = x*y/(x+y); u = c*w/y; v = c*w/x; lines labels end % reset_cb

function [p,q] = read_mouse % Current horizontal and vertical coordinates of the mouse. pq = get(gca,'currentpoint'); p = pq(1,1); q = pq(1,2); end % read_mouse

end % ladders_app

Published with MATLAB® R2016a

]]>Here is a classic puzzle. A pair of ladders leaning against the sides of an alley form a lopsided cross. Each ladder is propped against the base of one wall and leans against the opposite wall. If one ladder is 30 feet long, the other 20 feet long, and the point where they cross 10 feet above the ground, how wide is the alley?... read more >>

]]>Here is a classic puzzle. A pair of ladders leaning against the sides of an alley form a lopsided cross. Each ladder is propped against the base of one wall and leans against the opposite wall. If one ladder is 30 feet long, the other 20 feet long, and the point where they cross 10 feet above the ground, how wide is the alley?

The lengths of the ladders are denoted by $a$ and $b$. The height of the crossing is $c$. The heights of the points where the ladders reach the walls are $x$ and $y$. The bases of the right triangles whose common height is the crossing are $u$ and $v$. Finally, the width of the alley is $w$. Here is the picture.

```
ladders_diagram('vars')
```

These eight variables are connected by five equations. The Pythagorean Theorem provides two equations.

$$ a^2 = x^2 + w^2 $$

$$ b^2 = y^2 + w^2 $$

The ratios between sides of similar triangles provide two more.

$$ \frac{x}{w} = \frac{c}{v} $$

$$ \frac{y}{w} = \frac{c}{u} $$

The fifth equation is a simple linear relation that ties the two triangles together.

$$ w = u + v $$

Five equations involving eight variables leaves three degrees of freedom. In principle we could specify any three of the variables and solve for the other five. But that may or may not be easy. We will have two different situations in what follows. The puzzle posed in the opening paragraph specifies $a$, $b$ and $c$, and asks to find $w$. This is nontrivial. On the other hand, our graphic app is driven by $x$, $y$ and $w$. When these three are pinned down by mouse clicks, the other five follow immediately.

Combining the last three equations generates an elegant relation. This isn't a new equation; it's a consequence of the ones we already have.

$$ \frac{1}{c} = \frac{1}{x} + \frac{1}{y} $$

This says that the height of the crossing is one-half of the *harmonic mean* of the heights of the two ladders. The equation is a familiar one in *optics*, where it is known as the *thin lens equation*. It relates the location of an image to the object distance and focal length.

The following picture shows a solution of these equations where all eight variables have integer values. In fact, if we collect all eight variables into a vector $s$.

$$ s = [a, b, c, u, v, w, x, y] $$

And use the 1-norm to measure the "size" of a solution.

$$ ||s||_1 = \sum_{i=1}^8 s_i $$

Then ITOT ("It Turns Out That") this is smallest solution with all integer elements.

```
ladders_diagram('vals')
```

Combining Pythagoras and optics provides a single nonlinear equation for $w$ in terms of $a$, $b$ and $c$.

$$ \frac{1}{\sqrt{a^2 - w^2}}+\frac{1}{\sqrt{b^2 - w^2}}=\frac{1}{c} $$

It is not difficult to use the one-dimensional MATLAB zero finder `fzero` with fixed values for `a`, `b` and `c` to compute a solution of this equation.

a = 30; b = 20; c = 10; F = @(w) 1./sqrt(a^2 - w.^2) + 1./sqrt(b^2 - w.^2) - 1/c w = fzero(F,c)

F = @(w)1./sqrt(a^2-w.^2)+1./sqrt(b^2-w.^2)-1/c w = 12.3119

So this is the answer to the puzzle. For the prescribed values of the lengths of the ladders and the height of the crossing point, the width of the alley has to be about `12.3` feet.

A graph of `F(w)` in the vicinity of its zero shows why `fzero` has no trouble.

ezplot(F,[10,15]) set(gca,'xaxislocation','origin') line(w,0,'marker','.','markersize',24,'color','k')

It is easy to solve our single nonlinear equation numerically to compute $w$, but to find an analytic solution is a formidable challenge, even for computer algebra systems. Historically, the approach has focused on an equivalent quartic polynomial.

Take $a = 40$, $b = 30$ and $c = 20$. Let

$$ z = w^2 $$

Multiply through by the expressions in the denominators to put everything on one level.

$$ 10\, \sqrt{400 - z} + 10\, \sqrt{900 - z} = \sqrt{400 - z}\, \sqrt{900 - z} $$

Now square both sides to get rid of the sqrt's on the right. Then rearrange terms and square everything again to eliminate the remaining sqrt's. If you're careful, you will eventually reach a polynomial of degree 4 in $z$.

$$ z^4 - 2200\, z^3 + 1630000\, z^2 - 454000000\, z + 38500000000 = 0 $$

But we're only halfway there. We still have to solve this quartic. In principle it is possible to do this analytically, but let's again abandon algebraic and resort to numeric techniques.

poly = [ 1, -2200, 1630000, -454000000, 38500000000]; z = roots(poly)

z = 1.0e+02 * 8.4877 + 0.5881i 8.4877 - 0.5881i 3.5087 + 0.0000i 1.5158 + 0.0000i

The repeated squaring has produced extraneous roots. ITOT that we can recover $w$ from the fourth one.

w = sqrt(z(4))

w = 12.3119

It is reassuring to find the same width.

I am having a lot of fun with an interactive graphical experience where you vary any one of four parameters and see how it affects the others. You can change the height of the points where the ladders hit the wall, or change the width of the alley. You will find that dragging any one of these three control points affects only a couple of the other quantities.

You will see more action when you vary the crossing point. This changes all of the values except the width. Of course, it is not physically realistic to expect to alter the lengths of the ladders by changing where they cross, but that would actually have to happen if they were constrained to meet the walls.

The complete program for this app is the subject for my blog in two weeks. The title is "Investigating the Classic Crossed Ladders Puzzle". This link should be good after March 14.

Thanks to Ned Gulley of MathWorks for suggesting that this classic puzzle warrants another look.

Get
the MATLAB code

Published with MATLAB® R2016a

We say that a deck of playing cards is *completely shuffled* if it is impossible to predict which card is coming next when they are dealt one at a time. So a completely shuffled deck is like a good random number generator. We saw in my previous post that a perfect faro shuffle fails to completely shuffle a deck. But a *riffle shuffle*, with some randomness in the process, can produce complete shuffling. How many repeated riffle shuffles does that take?... read more >>

We say that a deck of playing cards is *completely shuffled* if it is impossible to predict which card is coming next when they are dealt one at a time. So a completely shuffled deck is like a good random number generator. We saw in my previous post that a perfect faro shuffle fails to completely shuffle a deck. But a *riffle shuffle*, with some randomness in the process, can produce complete shuffling. How many repeated riffle shuffles does that take?

*"Stop Shuffling and Deal."*

I learned about this topic from my friend Nick Trefethen, a Professor at Oxford and a MATLAB enthusiast. In 2000 he and his father wrote a paper, "How Many Shuffles to Randomize a Deck of Cards?". They refer to work by Persi Diaconis and colleagues. Diaconis is a Professor at Stanford and a world-class magician who has written extensively about the probability involved in activities like card shuffling and coin flipping.

I want to describe some simulations that reproduce the cutoff phenomenon investigated by the Trefethens and Diaconis.

The analysis is based on the number of *rising sequences* in the shuffled deck. A rising sequence is a maximal subset of cards in increasing numerical order. A deck is the union of its rising seguences. Here is an example, a small deck with just eight cards.

v = [3 1 4 5 7 2 8 6];

There are three rising sequences is this example.

[1 2]; [3 4 5 6]; [7 8];

We define the *signature* of a deck to be the number of rising sequences. Here is a tricky, but elegant, way to compute the signature of a deck `v`. Sort the deck to find the permutation `vbar` that reverses `v`, then count the number of places where `vbar` has a negative first difference.

```
type signature
```

function sig = signature(v); [~,vbar] = sort(v); sig = sum(diff(vbar)<0)+1; end

sigv = signature(v)

sigv = 3

Let's see how this works on the example.

[~,vbar] = sort(v) d = diff(vbar)<0; disp('d = ') disp([' ', sprintf('%6.0f', d)]), sig = sum(d)+1

vbar = 2 6 1 3 4 8 5 7 d = 0 1 0 0 0 1 0 sig = 3

An unshuffled deck of `n` cards, fresh out of the package, is represented by the vector `v = 1:n`. This has one long rising sequence, so the signature is 1. If you reverse the order, each card by itself is a singleton rising sequence, so the signature of the deck is `n`. These are the minimum and maximum possible values for the signature.

n = 52; v = 1:n; sig_min = signature(v) sig_max = signature(fliplr(v))

sig_min = 1 sig_max = 52

As we saw in my previous post, the "out" perfect shuffle keeps the first and last cards in first and last place, while the "in" perfect shuffle inserts these cards in the second and next-to-last places. Here are their permutation vectors.

m = n/2; pout = reshape([1:m; m+(1:m)],[1,n]); pin = reshape([m+(1:m); 1:m],[1,n]);

Both shuffles produce decks with two rising sequences obtained from the two halves of the deck. So their signatures are both equal to 2.

sig_pout = signature(v(pout)) sig_pin = signature(v(pin))

sig_pout = 2 sig_pin = 2

These signatures for `pout` and `pin` are way too small; the presence of just two rising sequences indicates perfect shuffles do not produce completely shuffled results.

A completely shuffled deck could be obtained by spreading the cards out on the table and then picking them up one at a time in a random order. This would be a random permutation of the integers `1:52`. Let's examine the signature of random permutations. We'll generate a million samples, then normalize the histogram of the resulting signatures so that it becomes a probability distribution, $\sigma_\infty$. This is the signature distribution of completely shuffled decks. Our goal is to get close to this distribution with riffle shuffling. (In order to reduce publication time, distributions are computed separately.)

```
type compute_sigma_inf
```

ksamples = 1000000; s = zeros(ksamples,1); for k = 1:ksamples v = randperm(n); s(k) = signature(v); end sigma_inf = hist(s,0:n)/length(s); save sigma_inf sigma_inf

```
load sigmas_52
plot(0:n,sigma_inf)
plot_details_inf
```

The distribution $\sigma_\infty$ is centered at 26 and 27 and falls off quickly from these two points. Here are the indices and values of all probabilities greater than .005. This tells us that completely shuffled decks have 26 or 27 increasing sequences 18.4% of the time each and that fewer than 21 or more than 32 increasing sequences are very unlikely.

support(sigma_inf,.005)

21 22 23 24 25 26 27 28 29 30 31 32 0.0061 0.0193 0.0480 0.0942 0.1471 0.1840 0.1841 0.1467 0.0936 0.0475 0.0196 0.0062

There are two sources of randomness in this model. The deck is cut at a depth determined by a random number chosen from a binomial distribution centered at its midpoint. Cards are then riffled back together. At each step the probability that the next card drops from either cut is proportional to the size of that cut.

```
type riffle
```

function vout = riffle(v) % Riffle Shuffle. % vout = riffle(v) % Cut deck roughly in half. n = length(v); m = binornd(n,0.5); % Binomial distribution % Interleave cuts. vout = zeros(1,n); i = 1; j = m+1; for k = 1:n % Choose which cut releases next card. % Probability is proportional to size of cut. if rand*(n-k+1) < m-i+1 vout(k) = v(i); i = i+1; else vout(k) = v(j); j = j+1; end end end

Here is one riffle shuffle of a fresh deck. The result almost certainly has a signature equal to 2, although the signature could be only 1 in the extremely unlikely event that the entire first cut falls in one clump before any of the second cut.

```
v = riffle(1:n);
deck_view(v)
title(['one riffle, signature = ' int2str(signature(v))])
```

Here are three riffles of a fresh deck. The result is likely to have a signature equal to 8, although any smaller value is possible.

```
v = riffle(riffle(riffle(1:52)));
deck_view(v)
title(['3 riffles, signature = ' int2str(signature(v))])
snapnow
close
```

Find the probability distribution generated by repeating a riffle shuffle *r* times. Denote this by $\sigma_r$.

```
type repeated_riffles
```

function sigma_r = repeated_riffles(n,r,ksamples) % sigma_r = repeated_riffles(n,r,ksamples) % Signature distribution of r repeated riffle shuffles. s = zeros(ksamples,1); for k = 1:ksamples v = 1:n; for j = 1:r v = riffle(v); end s(k) = signature(v); end sigma_r = hist(s,0:n)/length(s); end

For each value `r` between `0` and `rmax`, sample the signature $\sigma_r$ observed when a deck is subjected to `r` repeated riffle shuffles.

```
rmax = 14;
type compute_sigmas
```

sigma = zeros(rmax+1,n+1); ksamples = 1000000; for r = 0:rmax sigma(r+1,:) = repeated_riffles(n,r,ksamples) end save sigmas

```
load sigmas_52
```

For example, here is $\sigma_3$. Almost 95% of the distribution is concentrated at the maximum signature, 8. Slightly over 5% of the support has a signature in the range of 5 through 7.

```
r = 3;
plot(0:52,sigma(r+1,:),'.-')
plot_details_r
```

Here is the support of $\sigma_3$.

support(sigma(r+1,:),0)

5 6 7 8 0.000001 0.000424 0.050308 0.949267

Here are all of the $\sigma_r$. As expected, they appear to be heading towards $\sigma_{\infty}$.

```
plot(sigma','.-')
plot_details_sigma
```

The *total variation norm* is a very general concept that measures the distance between two probability distributions. In our situation it is simply the 1-norm of their difference, scaled by `1/2` so that it lies between 0 and 1.

norm_tv = @(s1,s2)(1/2)*norm(s1-s2,1)

norm_tv = @(s1,s2)(1/2)*norm(s1-s2,1)

For each value `r` between `0` and `rmax`, sample the signature $\sigma_r$ observed when a deck is subjected to `r` repeated riffle shuffles. Compare this with the completely shuffled distribution $\sigma_{\infty}$ and compute the total variation.

tv = zeros(rmax+1,1); for r = 0:rmax tv(r+1) = norm_tv(sigma(r+1,:), sigma_inf); end disp('tv =') fprintf('%12.6f\n',tv)

tv = 1.000000 1.000000 1.000000 1.000000 1.000000 0.923751 0.613194 0.333949 0.166832 0.084924 0.042375 0.021786 0.010285 0.004999 0.001935

Make an animated gif.

```
type make_tv_movie
```

gif_frame('html/tv_movie.gif'); for r = 0:rmax plot(0:n,[sigma(r+1,:); sigma_inf]','.-') plot_details_movie gif_frame(2) end clf gif_frame(2)

How many riffle shuffles are required to produce a completely shuffled deck? Plot the total variation as a function of the number of repetitions. The value of `tv` hovers very close to 1.0 for the first four shuffles, but then between five and nine shuffles it drops quickly to below 0.1. This is the *cutoff* phenomenon.

```
plot(0:rmax,tv,'o-')
plot_details_cutoff
```

It takes seven shuffles for the total variation to fall below 0.5. The conclusion is that about seven riffle shuffles are usually necessary and sufficient to produce a completely shuffled deck. Diaconis reaches this conclusion with a rigorous probabilistic analysis. Trefethen sees the cutoff with an analysis involving matrix powers and pseudospectra.

Here are the animation and the cutoff plot for four decks, so n = 208. The distribution of the random permutations, $\sigma_{\infty}$, is more spread out. But the shape of the cutoff is the same, it just occurs three shuffles later.

The conclusion here is that about ten riffle shuffles are usually necessary and sufficient to completely shuffle four decks.

n = 208; load sigmas_208 disp('tv =') fprintf('%12.6f\n',tv)

tv = 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 0.999555 0.913595 0.601832 0.328962 0.168368 0.083952 0.042310 0.021535

```
plot(0:rmax,tv,'o-')
plot_details_cutoff
```

Thanks very much to Nick Trefethen for help with this post.

L. N. Trefethen and L. M. Trefethen, "How Many Shuffles to Randomize a Deck of Cards?", *Proceedings: Mathematical, Physical and Engineering Sciences 456*, 2000. <http://www.jstor.org/stable/2665604>. <http://people.maths.ox.ac.uk/trefethen/publication/PDF/2000_87.pdf>.

Gudbjorn F. Jonsson and Lloyd N. Trefethen, "A Numerical Analyst Looks at the Cutoff Phenomenon in Card Shuffling and Other Markov Chains", *Numerical Analysis* (ed. D. F. Griffiths, D. J.Higmam and G. A. Watson), Addison-Wesley, 1998. <http://people.maths.ox.ac.uk/trefethen/publication/PDF/1997_74.pdf>.

Dave Bayer and Persi Diaconis, "Trailing the Dovetail Shuffle to its Lair", *Annals of Probability 2*, 1992. <http://projecteuclid.org/download/pdf_1/euclid.aoap/1177005705>.

A terrific series of videos from an interview with Persi Diaconis on Numberphile.

Get
the MATLAB code

Published with MATLAB® R2016a

When a deck of playing cards is shuffled perfectly, the result is not random. A perfect shuffle places the cards in a mathematically precise order. As a result, when the most common version of a perfect shuffle is repeated eight times, the deck returns to its original state.... read more >>

]]>When a deck of playing cards is shuffled perfectly, the result is not random. A perfect shuffle places the cards in a mathematically precise order. As a result, when the most common version of a perfect shuffle is repeated eight times, the deck returns to its original state.

Fresh out of its package, a deck of playing cards always has the Ace of Spades first, followed by the rest of the spades. The hearts, in order, are next, followed by the clubs (shown here in green), and finally the diamonds (blue). The King of Diamonds is last. Here is a view of a fresh deck.

v = 1:52; deck_view(v) title(0)

A perfect shuffle is also known as a *faro shuffle*, because it is used frequently in the classic card game, Faro. The deck is cut into two packs, each containing exactly half of the cards. Then the two packs are merged by precisely interleaving the cards.

There are two variants. In an *out-faro* shuffle, the merger is started with the top card from the first half and eventually completed with the last card from the second half. So with an out-faro shuffle of a fresh deck, the Ace of Spades remains on top and the King of Diamonds remains on the bottom. An *in-faro* shuffle is started with the second half. So the Ace of Spades and the King of Diamonds are inserted in the second and next to last positions respectively.

A shuffle is a permutation of the elements of a vector representing the deck. Here is the index vector that produces an out-faro shuffle. A rectangular matrix is created with `1:26` in the first row and `27:52` in the second. So its first column is `[1;27]`, its second column is `[2;28]`, and so on. The `reshape` operation, which is carried out by columns, then generates the desired index vector.

pout = reshape([1:26; 27:52],[1,52])

pout = Columns 1 through 13 1 27 2 28 3 29 4 30 5 31 6 32 7 Columns 14 through 26 33 8 34 9 35 10 36 11 37 12 38 13 39 Columns 27 through 39 14 40 15 41 16 42 17 43 18 44 19 45 20 Columns 40 through 52 46 21 47 22 48 23 49 24 50 25 51 26 52

Here is the result of two successive out-faro shuffles of a fresh deck.

v = 1:52; v = v(pout); deck_view(v) title(1)

v = v(pout); deck_view(v) title(2)

The four Aces have moved to the front and the four Kings are bringing up the rear.

Let's start over with a fresh deck, apply the `pout` permutation eight times, and capture the results in an animated gif. The deck is returned to its original state in just 8 steps.

v = 1:52; deck_view(v) title(0) for t = 1:8 v = v(pout) deck_view(v) title(t) end

Here is the indexing vector for the in-faro suffle.

pin = reshape([27:52; 1:26],[1,52])

pin = Columns 1 through 13 27 1 28 2 29 3 30 4 31 5 32 6 33 Columns 14 through 26 7 34 8 35 9 36 10 37 11 38 12 39 13 Columns 27 through 39 40 14 41 15 42 16 43 17 44 18 45 19 46 Columns 40 through 52 20 47 21 48 22 49 23 50 24 51 25 52 26

It takes 52 in-faro shuffles to get back to the original state.

Both of the faro shuffles can be represented by permutation matrices. Can you see how their `spy` plots differ?

```
I = eye(52);
close
Pout = I(pout,:);
spy(Pout)
title('Pout')
```

```
Pin = I(pin,:);
spy(Pin)
title('Pin')
```

It turns out that the smallest value of `t` for which the matrix power `P^t` is equal to the identity matrix is `t = 8` for `P = Pout` and `t = 52` for `P = Pin`.

All of this is explained by *eigenvalues*. The matrix `Pout` has order 52, but only 8 distinct eigenvalues, namely the 8-th roots of unity.

plot(eig(Pout),'o') title('eig(Pout)') axis(1.25*[-1 1 -1 1]) axis('square')

On the other hand, `Pin` has 52 distinct eigenvalues, the 52-nd roots of unity.

plot(eig(Pin),'o') title('eig(Pin)') axis(1.25*[-1 1 -1 1]) axis('square')

The eigenvalues of `Pout` occur with different multiplicites. Can you explain where these counts come from?

e = eig(Pout); z = exp(pi/4*1i) count = zeros(1,8); for k = 0:7 count(k+1) = sum(abs(e-z^k)<52*eps); end disp(' ') disp('multiplicites = ') disp(' 1 z i z^3 -1 z^5 -i z^7') disp(count)

z = 0.7071 + 0.7071i multiplicites = 1 z i z^3 -1 z^5 -i z^7 9 6 6 6 7 6 6 6

Here is the listing of `deck_view`.

```
type deck_view
```

function deck_view(v) % deck_view(v) displays the card deck v. % Prepare the figure window. f = clf; f.Position = [200 300 600 150]; f.NumberTitle = 'off'; f.ToolBar = 'none'; f.MenuBar = 'none'; f.PaperPositionMode = 'auto'; % Prepare the axes. ax = axes; ax.Position = [0 0 1 .88]; ax.XLim = [-1 54]; ax.YLim = [-1 15]; ax.XTick = ''; ax.YTick = ''; ax.Box = 'on'; % Colors for spades, hearts, clubs, diamonds. shcd = [0 0 0; 1 0 0; 0 .5 0; 0 0 .5]; % Plot one text character per card. pips = 'A23456789TJQK'; for i = 1:52 j = mod(v(i)-1,13)+1; k = floor((v(i)-1)/13)+1; t = text(i,j,pips(j)); t.FontSize = 10; t.Color = shcd(k,:); end end

Kevin Houston is a mathematician at the University of Leeds in the UK who can also do perfect shuffles. Here is a link to his YouTube video. Perfect Shuffle.

Get
the MATLAB code

Published with MATLAB® R2015b