Find the optimum traveling salesman tour through the capitals of the 48 contiguous states in the USA.... read more >>

]]>Find the optimum traveling salesman tour through the capitals of the 48 contiguous states in the USA.

I am continuing with the workspace from my previous blog post.

John Burkardt provides a second data set, state_capitals_ll.txt, giving the longitudes and latitudes of capitals of the 48 contiguous states. (Ignore AK, DC, HI, PR and US.)

[lon,lat] = get_capital_coordinates;

When I use these coordinates for the graph I created in my previous post, we get a more accurate picture of the USA.

G = graph(A,cellstr(names)); plot(G,'xdata',lon,'ydata',lat);

The travelling salesman problem, the TSP, was mathematically formulated in the 19th century. The problem is to find the closed circuit of a list of cities that travels the shortest total distance.

For 25 years MATLAB releases have included a simple demo program named `travel` that finds an approximate solution to the TSP through a few dozen randomly chosen points. The background of the plot, an outline the USA, is just for artistic effect.

I've made the `travel` demo into a function named `traveler` whose input is the distance matrix `D`, the pairwise distances between the cities. The output is a path length `miles` and permutation vector `p` so that `miles` is the length of the route produced by `lon(p)` and `lat(p)`.

Each time it is called, `travel` starts with a random permutation of the cities. Then, for several thousand iterations, it tries to reduce the length of the trip by revising the permutation in each of two simple ways. One scheme is to move one city to another position in the ordering. This example starts with the order `[1:7]`. It then moves city `4` to follow city `5`, so the order becomes `[1 2 3 5 4 6 7]`. This reduces the length from `9.71` to `7.24`.

half_fig traveler_scheme1

The second scheme is based on the observation that any time a route crosses itself, the length can be reduced by reversing the crossed segment. We don't actually look for crossings, we just try reversing random chunks of the ordering. This example reverses the `[3 4 5]` in `[1:7]` to get `[1 2 5 4 3 6 7]`, thereby reducing the length from `9.47` to `7.65`.

traveler_scheme2

These two heuristics are not powerful to guarantee that `traveler` will find the global shortest path. It's easy to get stuck in local minima. And the random path approach means that `traveler` is inefficient for more than a few dozen cities.

But nobody knows how to find a guaranteed shortest path for the TSP and the code for `traveler` is only 58 lines. I will include the code at the end of this post and in Cleve's Laboratory on the MATLAB Central File Exchange.

close

The input to `traveler` is the distance matrix.

D = distances(lon,lat);

The inputs to the function `distances` are the vectors `lon` and `lat`, the longitude and latitude of the cities. The output is the 48-by-48 matrix `D`, the great circle distance between pairs of cities. With some spherical trigonometry, `lon` and `lat`, which are essentially angles measured in degrees, are converted to miles, the lengths of "as the crow flies" geodesics on the surface of the earth. The radius of the earth is a scale factor. Here is the core of the function.

dbtype 11:20 distances.m

11 R = 3959; % radius of the earth (mi.) 12 n = length(lon); 13 D = zeros(n,n); 14 for k = 1:n 15 for j = 1:k-1 16 D(k,j) = R*acos(sind(lat(k))*sind(lat(j)) + ... 17 cosd(lat(k))*cosd(lat(j))*cosd(lon(k)-lon(j))); 18 D(j,k) = D(k,j); 19 end 20 end

Which two capitals are nearest each other?

Dmin = min(min(D(D>0))) [k,j] = find(D == Dmin) cities = names(k)

Dmin = 34.9187 k = 37 19 j = 19 37 cities = 2×1 string array "RI" "MA"

Providence, Rhode Island and Boston, Massachusetts ("our fair city") are less than 35 miles apart.

What about the other extreme?

Dmax = max(max(D)) [k,j] = find(D == Dmax) cities = names(k)

Dmax = 2.6629e+03 k = 17 4 j = 4 17 cities = 2×1 string array "ME" "CA"

As we might have expected, the capitals of Maine and California are the farthest apart, 2663 miles. That would require one prodigious crow.

Based on some experience that I will describe shortly, I am going to set the random number seed to `347` before we take our road trip with `traveler`.

rng(347) [miles,p] = traveler(D); miles = round(miles)

miles = 10818

That's an average of

avg = miles/48

avg = 225.3750

miles per leg.

Let's highlight our graph of neighboring states with our traveling salesman path linking capital cities.

Gp = plot(G,'xdata',lon,'ydata',lat,'edgecolor',turquoise); highlight(Gp,p,'edgecolor',darkred,'linewidth',2) title(miles)

There are four missing links. The path goes from UT to MT, but those states are not neighbors. We must go through ID. And the leg from FL to SC goes through GA. There are two more in the northeast. Fill in those missing links with dashed lines.

missing_links(A,lon,lat,p,darkred)

Zoom in on the northeast.

axis([-84 -68 33 46])

Here the state abbreviations, ordered by this path.

fprintf(fmt,names(p))

LA TX OK NM AZ NV CA OR WA ID MT UT CO WY SD ND MN WI MI OH WV PA NY VT ME NH MA RI CT NJ DE MD VA NC SC FL AL GA TN KY IN IL IA NE KS MO AR MS LA

It is interesting to look at the last step. After 3298 iterations, `travel` arrives at this situation with a path length of 10952.

The link connecting WV to VA crosses the link connecting NC to PA. It takes `travel` 222 more iterations, but eventually a permutation that reverses the path from VA and PA is tried. This connects WV to PA and NC to VA to produce the path shown in the previous plot. The total path length is decreased by 134 miles to 10818.

I am confident that we have found the shortest path, but I don't have a proof.

I ran 1000 different initializations of `rng`. Eight runs produced the path with length 10818 that we found with `rng(347)`. None of the runs found a shorter path. Sixteen more runs found a quite different path with length 10821, only three miles longer. Here is this next best path.

I must include an animated gif here. This initially shows every sixth step. It switches to every step when the path length is less than 12000. Few legs in random paths match edges in the neighbor graph so initially most lines are dashed. As we near the minimum, more solid lines appear.

Chances that, without help, `traveler` will find the shortest route through the 48 capital cities are less than 1 in 100. Here is your opportunity to help guide the search. The `traveler_game` is a modification of `travel` that plots every successful step and that provides controls so that you to back up a few steps and try the random strategy again.

I will include the `traveler_game` in the next update, version 3.70, of Cleve's Laboratory on the MATLAB Central File Excange. That may take a few days.

There are five push buttons.

- > Take one minimization step.
- >> Take repeated steps, until a local minimum is reached.
- < Reverse one step.
- << Reverse many steps.
- ^ Start over in a new figure window.

Here are four typical runs.

Here are the codes for `traveler`, `distances`, and `path_length`.

type traveler type distances type path_length

function [len,p] = traveler(D) %TRAVELER Functional form of old MATLAB demo, "travel". % A pretty good, but certainly not the best, solution of the % Traveling Salesman problem. Form a closed circuit of a % number of cities that travels the shortest total distance. % % [len,p] = traveler(D). % Input: D = distances between cities with coordinates x and y. % Output: p a permutation, x(p) and y(p) is a path with length len. % Copyright 1984-2018 The MathWorks, Inc. n = size(D,1); p = randperm(n); len = path_length(p,D); for iter = 1:10000 % Try to reverse a portion. pt1 = floor(n*rand)+1; pt2 = floor(n*rand)+1; lo = min(pt1,pt2); hi = max(pt1,pt2); q = 1:n; q(lo:hi) = q(hi:-1:lo); pnew = p(q); lennew = path_length(pnew,D); if lennew < len p = pnew; len = lennew; iterp = iter; end % Try a single point insertion pt1 = floor(n*rand)+1; pt2 = floor((n-1)*rand)+1; q = 1:n; q(pt1) = []; q = [q(1:pt2) pt1 q((pt2+1):(n-1))]; pnew = p(q); lennew = path_length(pnew,D); if lennew < len p = pnew; len = lennew; iterp = iter; end end % Close the permutation. p(end+1) = p(1); end function D = distances(lon,lat) % D = distances(lon,lay) % Input: vectors lon and lat, the longitude and latitude of the cities. % Output: D(k,j) is the distance (in miles) between cities k and j. % Copyright 1984-2018 The MathWorks, Inc. % Great circle distance matrix between pairs of cities. % https://en.wikipedia.org/wiki/Great-circle_distance. R = 3959; % radius of the earth (mi.) n = length(lon); D = zeros(n,n); for k = 1:n for j = 1:k-1 D(k,j) = R*acos(sind(lat(k))*sind(lat(j)) + ... cosd(lat(k))*cosd(lat(j))*cosd(lon(k)-lon(j))); D(j,k) = D(k,j); end end end function len = path_length(p,D) % len = path_length(p,D) % Calculate current path length for traveling salesman problem. % This function calculates the total length of the current path % p in the traveling salesman problem. % % This is a vectorized distance calculation. % % Create two vectors: p and q = p([n 1:(n-1)]). % The first is the list of first cities in any leg, and the second % is the list of destination cities for that same leg. % (the second vector is an element-shift-right version of the first) % % Use column indexing into D to create a vector of % lengths of each leg which can then be summed. n = length(p); q = p([n 1:(n-1)]); len = sum(D(q + (p-1)*n)); end

Get
the MATLAB code

Published with MATLAB® R2018a

Create a MATLAB® graph object from information about neighbors among the 48 contiguous states in the USA.... read more >>

]]>Create a MATLAB® graph object from information about neighbors among the 48 contiguous states in the USA.

A few years ago, MathWorks asked Tom and Ray Magliozzi, MIT alumni and hosts of the popular show Car Talk on National Public Radio, to come to our offices and meet some visiting students. They said yes on one condition. They figured MathWorkers would be a good source of new material for the puzzler segment of their show. It turned out that we had only one puzzle to offer that they hadn't heard before. But one was enough, and they came for the visit.

Here is that puzzle.

*
My friend tells me that she is about to take a road trip.
She plans to visit each of the 48 contiguous states in the
USA exactly once. She tells me she is going to start someplace
in the West and doesn't tell me where she is going to finish.
However, I am able to agree to meet her in the last state she
will visit. Which state is that?
*

I will reveal the state later in this post, but try to answer the puzzle before I do. (Readers who are not familiar with USA geography don't have to worry -- few Americans can provide the answer.)

The web site of John Burkardt at Florida State University is a tremendous repository of software, data sets and other goodies. I am going to use state neighbors. The first column in this file provides the two-letter abbreviations of the names of the 48 contiguous states.

S = get_state_neighbors; names = S{:,1}; fmt = [repmat('%s ',1,16) '\n']; fprintf(fmt,names)

AL AZ AR CA CO CT DE FL GA ID IL IN IA KS KY LA ME MD MA MI MN MS MO MT NE NV NH NJ NM NY NC ND OH OK OR PA RI SC SD TN TX UT VT VA WA WV WI WY

The list is ordered alphabetically by the full state name, so Wyoming is last.

last = names(end)

last = "WY"

Burkardt's data also gives us the neighbors of each state. We can form an adjacency matrix `A` where `A(i,j)` is one if states `i` and `j` share a common border and zero otherwise.

A = adjacency(S);

For example, Wyoming has six neighbors.

last_neighbors = names(A(end,:))

last_neighbors = 6×1 string array "CO" "ID" "MT" "NE" "SD" "UT"

Every good adjacency deserves a `spy` plot.

spy(A)

The last row (and column) has six nonzeros for Wyoming's six neighbors.

Which states have the most neighbors?

num_neigh = sum(A); max_num_neigh = max(num_neigh) max_neigh = names(num_neigh == max_num_neigh)

max_num_neigh = 8 max_neigh = 2×1 string array "MO" "TN"

Missouri and Tennessee each have eight neighbors.

Which state has the fewest neighbors?

min_num_neigh = min(num_neigh) min_neigh = names(num_neigh == min_num_neigh)

min_num_neigh = 1 min_neigh = "ME"

Maine has only one neighbor, namely New Hampshire. This provides the answer to the MathWorks Car Talk puzzle. My friend's trip must either start or finish in in Maine. Once she's in Maine, she can't leave without going through New Hampshire. Since she is starting in the West, she must finish in Maine.

Every adjacency matrix also has an associated graph.

G = graph(A,cellstr(names))

G = graph with properties: Edges: [105×1 table] Nodes: [48×1 table]

The default layout for plotting this graph is `force`, which uses attractive and repulsive forces on the nodes.

plot(G)

Interchange the axes by switching the viewpoint.

view(-90,-90)

The plot now resembles the United States. California, Oregon and Washington are on the west coast. Florida is a peninsula in the southeast. The New England states dominate the northeast. MO and TN show off their many neighbors.

Utah, Arizona, New Mexico and Colorado are the corners of one of two quadrilaterals in the graph, mimicking their actual meeting at Four Corners. The other quadrilateral involves the four states that surround Lake Michigan.

All this without any coordinates, just knowing the connectivity of neighboring states.

My next post will add the coordinates of the state capitals and tackle the classic Traveling Salesman Problem.

Thanks to John Burkhardt for his data sets and for pointing out that the wording of the puzzle needs to preclude starting in Maine.

Get
the MATLAB code

Published with MATLAB® R2018a

I probably first saw this matrix in 1960 in John Todd's class at Caltech. But I forgot about it until Tahar Loulou jogged my memory with a comment following my blog post in late May. It deserves a place in our gallery of interesting matrices.... read more >>

]]>I probably first saw this matrix in 1960 in John Todd's class at Caltech. But I forgot about it until Tahar Loulou jogged my memory with a comment following my blog post in late May. It deserves a place in our gallery of interesting matrices.

I'm afraid that I don't know very much about T. S. Wilson himself. The only paper we could find on the Internet is the joint paper about pendulum dampers listed in the references. This paper, which does not mention the matrix, was published in 1940. The authors were from a British company called D. Napier & Son Ltd. The earliest appearance of the matrix that we know of is the 1946 paper by J. Morris referenced below. It is possible that by this time Wilson and Morris were colleagues at the Royal Aircraft Establishment in Farnborough in the UK.

I think that Wilson and Todd must have known each other. Todd acknowledges Wilson in several papers and books.

One of the references we used in Todd's class is the booklet written by Marvin Marcus also referenced below. The booklet has a short section collecting some "particular" matrices. A_13 is our old friend the Hilbert matrix and A_14 is Wilson's matrix.

Here is the matrix. Let's use the Live Editor® and the Symbolic Math Toolbox®.

W = sym([5 7 6 5; 7 10 8 7; 6 8 10 9; 5 7 9 10]);

The first noteworthy property is the determinant. The determinant has to be an integer because all the matrix elements are integers.

d = det(W);

This is the smallest possible determinant for a nonsingular matrix with integer entries.

Since the determinant is one, the inverse must also have integer entries.

Winv = round(inv(W));

The matrix is symmetric. Let's use `chol` to verify that it is positive definite.

try R = chol(W); catch disp("Not SPD.") end

Since `chol` was successful we know that `W` is positive definite. Factor out the diagonal to get a nice rational $LDL^T$ factorization.

D = diag(diag(R)); L = R'/D; D = D^2; LDLT = L*D*L';

The characteristic polynomial is a quartic with integer coefficients.

```
p = charpoly(W,'lambda');
```

Since the degree of *p* is less than 5, the exact eigenvalues can be expressed in closed form.

e = eig(W);

We asked for it. Like it or not, that's the exact answer. Since `p` does not factor nicely over the rationals, the closed form is intimidating. It's hard to tell which eigenvalue is the largest, for example. Notice that even though the eigenvalues are real, the formulae involve complex numbers whose imaginary parts necessarily cancel.

For the remainder of this article, let's revert to floating point arithmetic.

W = double(W)

W = 5 7 6 5 7 10 8 7 6 8 10 9 5 7 9 10

Now the computed eigenvalues are approximate, but more useful.

e = eig(W)

e = 0.0102 0.8431 3.8581 30.2887

Is this matrix poorly conditioned?

condW = cond(W)

condW = 2.9841e+03

So the condition is "only" about 3000. That may not seem to be very badly conditioned. But let's look more closely. Set up a linear system where we know the exact solution.

xact = ones(4,1); b = W*xact

b = 23 32 33 31

And solve it.

x = W\b

x = 1.0000 1.0000 1.0000 1.0000

We expect to get all ones, and we do. But now perturb the right-hand side a little.

delta_b = [.01 -.01 .01 -.01]' b_prime = b + delta_b

delta_b = 0.0100 -0.0100 0.0100 -0.0100 b_prime = 23.0100 31.9900 33.0100 30.9900

This perturbation changes the solution drastically.

x_prime = W\b_prime

x_prime = 1.8200 0.5000 0.8100 1.1100

I want to argue that Wilson's matrix *is* badly conditioned, when compared to matrices with these properties.

- 4-by-4
- Symmetric
- Nonsingular
- Integer entries between 1 and 10.

My experiment involves generating one million random matrices with these properties. Random symmetric matrices are created by

U = triu(randi(10,4,4)); A = U + triu(U,1)';

One typical run produced these statistics.

- Number of matrices = 1000000
- Singular = 1820
- Positive definite = 11217
- Condition > condW = 2115
- Maximum condition = 4.8087e4

The distribution of the condition numbers looks like this. We can see that Wilson's matrix is unusual. Only about 0.21 percent of comparable matrices have a larger condition number.

By the way, the matrix with the largest condition number in this run was

condMax

condMax = 1 3 10 10 3 4 8 9 10 8 3 9 10 9 9 3

Thanks to Tahar Loulou whose comment got this all started, to Steve Lord who contributed another comment, and to Jan Peter Schäfermeyer, Nick Higham and Sven Hammarling for their email pointing to references.

R. W. Zdanowich and T. S. Wilson, "The elements of pendulum dampers," *Proceedings of the Institution of Mechanical Engineers*, 143_028_02, pp. 182-210, 1940, link.

J. Morris, "XVI. An escalator process for the solution of linear simultaneous equations", *The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science*, 37:265, pp. 106-120, 1946, link.

Jacques Delforge, "Application of the concept of the conditioning of a matrix to an identification problem," *Comp. and Maths with Appls.*, Vol. 4, pp. 139-148, 1978, link.

E. Bodewig, Matrix calculus, North Holland, Amsterdam, 2nd edition, 1959.

John Todd, "The condition of a certain matrix," *Proceedings of the Cambridge Philosophical Society*, vol. 46, pp. 116-118, 1950, link.

John Todd, Basic Numerical Algorithms: Vol. 2: Numerical Algorithms, Birkhäuser Verlag, Basel, 1977.

Marvin Marcus, Basic Theorems in Matrix Theory, National Bureau of Standards, Applied Mathematics Series, number 57, page 22, 1960, link.

Get
the MATLAB code

Published with MATLAB® R2018a

MIT's Professor Gil Strang gave two talks in one morning recently at the SIAM annual meeting. Both talks derived from his experience teaching a new course at MIT on linear algebra and neural nets. His first talk, "The Structure of a Deep Neural Net", was in a minisymposium titled "Deep Learning and Deep Teaching", which he organized. Another talk in that minisymposium was by Drexel's Professor Pavel Grinfeld on "An Informal Approach to Teaching Calculus." An hour later, Gil's gave his second talk, "Teaching About Learning." It was an invited talk at the SIAM Conference on Applied Mathematics Education.... read more >>

]]>MIT's Professor Gil Strang gave two talks in one morning recently at the SIAM annual meeting. Both talks derived from his experience teaching a new course at MIT on linear algebra and neural nets. His first talk, "The Structure of a Deep Neural Net", was in a minisymposium titled "Deep Learning and Deep Teaching", which he organized. Another talk in that minisymposium was by Drexel's Professor Pavel Grinfeld on "An Informal Approach to Teaching Calculus." An hour later, Gil's gave his second talk, "Teaching About Learning." It was an invited talk at the SIAM Conference on Applied Mathematics Education.

Inspired by Pavel's talk about teaching calculus, Gil began his second talk with some spontaneous remarks. "Can we teach a deep learner calculus?" he wondered. The system might be trained with samples of functions and their derivatives and then be asked to find derivatives of other functions that were not the training set.

Immediately after Gil's second talk, I asked the other MathWorkers attending the meeting if we could take up Gil's challenge. Within a few hours, Mary Fenelon, Christine Tobler and Razvan Carbunescu had the essential portions of the following demonstration of the Neural Networks Toolbox® working at the MathWorks booth.

Our first task is less ambitious than differentiation. It is simply to recognize the shapes of functions. By a function we mean the MATLAB vector obtained by sampling a familiar elementary function at a finite set of ordered random points drawn uniformly from the interval $[-1, 1]$. Derivatives, which we have not done yet, would be divided differences. We use six functions, $x$, $x^2$, $x^3$, $x^4$, $\sin{\pi x}$, and $\cos{\pi x}$. We attach a random sign and add white noise to the samples.

F = {@(x) x, @(x) x.^2, @(x) x.^3, @(x) x.^4, ... @(x) sin(pi*x), @(x) cos(pi*x)}; labels = ["x", "x^2", "x^3", "x^4", "sin(pi*x)", "cos(pi*x)"];

Here are a few definitions and parameters.

Set random number generator state.

rng(2)

Generate uniform random variable on [-1, 1].

randu = @(m,n) (2*rand(m,n)-1);

Generate random +1 or -1.

randsign = @() sign(randu(1,1));

Number of functions.

m = length(F);

Number of repetitions.

n = 1000;

Number of samples in the interval.

nx = 100;

Noise level.

noise = .0001;

C = cell(m,n); for j = 1:n x = sort(randu(nx,1)); for i = 1:m C{i,j} = randsign()*F{i}(x) + noise*randn(nx,1); end end

Let's plot instance one of each function. (With this initialization of `rng`, the $x^2$ and $sin(\pi x)$ curves have negative signs.)

set(gcf,'position',[300 300 300 300]) for i = 1:m plot(x,C{i,1},'.') axis([-1 1 -1 1]) title(labels(i)) snapnow end close

Our deep learning network is one that has proved to be successful in signal processing, text, and other applications with sequential data. There are six layers. The nonlinear activation layer, `relu`, for REctified Linear Unit, is essential. `ReLU(x)` is simply `max(0,x)`. LSTM stands for Long Short-Term Memory. Softmax is a generalization of the logistic function used to compute probabilities.

inputSize = nx; numClasses = m; numHiddenUnits = 100; layers = [ ... sequenceInputLayer(inputSize) reluLayer lstmLayer(numHiddenUnits,'OutputMode','last') fullyConnectedLayer(numClasses) softmaxLayer classificationLayer];

The first option, `'adam'`, is the stochastic optimization algorithm, adaptive moment estimation. An epoch is one forward pass and one backward pass over all of the training vectors, updating the weights. In our experience with this network, six passes is enough.

maxEpochs = 6; miniBatchSize = 27; options = trainingOptions('adam', ... 'ExecutionEnvironment','cpu', ... 'MaxEpochs',maxEpochs, ... 'MiniBatchSize',miniBatchSize, ... 'GradientThreshold',1, ... 'Verbose',0, ... 'Plots','training-progress');

With our setting of the `Plots` option, `trainNetwork` opens a custom figure window that dynamically shows the progress of the optimation.

It takes 15 or 20 seconds to train this network on my laptop. Big time neural nets with more layers and more epochs can make use of GPUs and pools of parallel workers

C = reshape(C',1,[]); Y = repelem(categorical(labels'), n); net = trainNetwork(C,Y,layers,options);

Generate more functions to form a test set.

nt = 100; Ctest = cell(m,nt); for j = 1:nt x = sort(randu(nx,1)); for i = 1:m Ctest{i,j} = randsign()*F{i}(x) + noise*randn(nx,1); end end

Classify the functions in the test set.

```
miniBatchSize = 27;
Ctest = reshape(Ctest',1,[]);
Ytest = repelem(categorical(labels'), nt);
Ypred = classify(net,Ctest,'MiniBatchSize',miniBatchSize);
```

Here are the results. We see scores above 95 percent, except for learning to distinguish between plots of $x^2$ and $x^4$. That's understandable.

T = table(Ypred, Ytest); heatmap(T, 'Ypred', 'Ytest');

Here are typical failed tests. It's not hard to see why the network is having trouble with these.

set(gcf,'position',[300 300 300 300]) for j = 1:m for i = [1:j-1 j+1:m] mismatch = find(Ytest == labels(i) & Ypred == labels(j)); if ~isempty(mismatch) % Plot one of each type of mismatch for k = 1 plot(linspace(-1, 1), Ctest{mismatch(k)},'.') title(labels(i)+" or "+labels(j)) snapnow end end end end close

Thanks to my colleagues attending the SIAM meeting -- Mary, Christine and Razvan. And thanks to Gil for his spontaneous idea.

Get
the MATLAB code

Published with MATLAB® R2018a

David Wilson, from the Auckland University of Technology in New Zealand, has alerted me to this remarkable coincidence. Here is the "Burning Ship" from my post about fractals four weeks ago.... read more >>

]]>David Wilson, from the Auckland University of Technology in New Zealand, has alerted me to this remarkable coincidence. Here is the "Burning Ship" from my post about fractals four weeks ago.

And here is a real photo, entitled *The Long, Long Night*, taken in August 1915 by Frank Hurley, of Antarctic explorer Ernest Shakleton's ship *Endurance*, frozen in the ice in the Weddell Sea.

Photo credit: National Library of Australia, http://nla.gov.au/nla.obj-158931586/view.

Get
the MATLAB code

Published with MATLAB® R2018a

Here are the first 12 integers from an infinite sequence defined by a deceptively simple rule. Can you see the pattern? Try to predict the next number in the sequence.... read more >>

]]>Here are the first 12 integers from an infinite sequence defined by a deceptively simple rule. Can you see the pattern? Try to predict the next number in the sequence.

0, 1, 3, 6, 2, 7, 13, 20, 12, 21, 11, 22

The OEIS, the On-Line Encyclopedia of Integer Sequences, is a mathematical treasure chest. Established in 1964 by Neal J. A. Sloan, the OEIS has over 300,000 entries. Dozens are added each week. Anyone who registers can propose a new entry. A team of 130 volunteer editors then review the proposals and create new entries.

A powerful search engine accepts queries made from short subsequences and finds all the entries in the encyclopedia that contain that query. It's similar to a search engine that identifies a piece of music from a few notes.

Here is the link, https://oeis.org. One way to introduce yourself to the OEIS is to take a look at https://oeis.org/OEIS_pics.html. Or, watch the OEIS movie on YouTube, https://www.youtube.com/watch?v=LCWglXljevY.

Did you try to guess the next number is our sequence yet? You can ask the OEIS. If you enter the 12 numbers in the search window, you get something like this.

So, you can see that the next number is 10.

Each entry in the encyclopedia has an identifier; for this sequence the identifier is A005132. We could also reach A005132 by searching for "A005132" or "Recaman".

Neal Sloan has said this is one of his favorite sequences in the entire collection. He learned about it in a 1991 letter from Columbian mathematician Bernardo Recamán Santos.

In the past few weeks I've become infatuated with Recamán sequence.

Starting with $a_0$ = 0, the sequence is generated by the recurrence

$$a_n = a_{n-1} \pm n$$

The minus sign is preferred. It is used if the resulting $a_n$ is positive and not already in the sequence; otherwise the plus sign is used.

So, for n = 1, 2, 3 the plus sign is required and the rule generates the familiar triangle numbers 1, 3, 6. But when n = 4 it is possible to use the minus sign because

$a_4$ = $a_3$ - 4 = 6 - 4 = 2

is not already in the sequence.

Here is my MATLAB function.

```
type recaman
```

function R = recaman(m) % R = recaman(n) is the first n terms of Recamán's sequence. % OEIS A005132, https://oeis.org/A005132. a = 0; R = a; for n = 1:m-1 if a > n && ~ismember(a-n,R) a = a-n; else a = a+n; end R(n+1) = a; end end

It turns out to be faster for the sequence array `R` to grow during the recurrence than it is to preallocate it with `R = zeros(1,m)`.

Let's visualize the Recamán sequence. A plot made with semicircles connecting the terms is especially illuminating. The first small semicircle, on the left, has diameter one and end-points at the first two elements, 0 and 1. Subsequently, the diameter of the n-th semicircle is n and its end-points are $a_{n-1}$ and $a_n$. This plot ends on the right at $a_{12}$ = 22

I learned about this way of plotting the Recamán sequence from the YouTube channel Numberphile in a video featuring Alex Bellos and Edmund Harriss.

```
set(gcf,'position',[200 200 480 360])
n = 12;
recaman_circles(n)
```

A conventional MATLAB plot of the sequence versus index is also worth a look. The title is the last term, $a_n$.

recaman_plot(n)

The value n = 66 is a good stopping point for the plots. The next value would jump from $a_{66}$ = 91 to $a_{67}$ = 157.

The line plot is showing its characteristic behavior -- stretches of sawtooth oscillation with increasing amplitude separated by chaotic jumps to new "energy levels".

This behavior with 300 terms continues indefinitely.

n = 300; recaman_plot(n)

The computational complexity of `recaman(n)` is `O(n^2)`. It takes about six minutes to generate a million terms on my laptop.

The Recamán recurrence is designed to hit every positive integer. Does it succeed? For each k > 0, is there a value of n for which $a_n$ = k. This is an open question. Sometimes n has to be surprisingly large.

Take k = 4.

k = 4 R = recaman(1000); n = find(k == R)

k = 4 n = 132

Take k = 19.

k = 19 R = recaman(100000); n = find(k == R)

k = 19 n = 99735

Take k = 852655. Wait, don't try that. Nobody knows if the sequence ever hits this value. Benjamin Chaffin has computed 10^230 terms without seeing 852655.

I can't easily include audio in my blog. You'll have to do it yourself. For a real treat, go to https://oeis.org/A005132 and click on "listen".

Let's have a little more fun with the OEIS. There are only a handful of people in the world that know the next term in this sequence.

9, 6, 3, 9, 7, 2, 3, 8

To find out, enter these numbers in the search window at https://oeis.org. Or, bring up your MATLAB and

```
edit membrane
```

Get
the MATLAB code

Published with MATLAB® R2018a

If you follow this blog regularly, you know that I love fractals. I recently spent a pleasant afternoon in Nashua, New Hampshire, where my daughter Teresa introduced me to Gregory Searle, a fractal artist and computer geek. Here is his logo.... read more >>

]]>If you follow this blog regularly, you know that I love fractals. I recently spent a pleasant afternoon in Nashua, New Hampshire, where my daughter Teresa introduced me to Gregory Searle, a fractal artist and computer geek. Here is his logo.

Greg's web page is <http://www.fractalartdesign.com>. There is much to explore. To get started, click on "Fractal Generator" in the menu on the left of the page, and then on the "Fractal Generator" icon in the center of the next page. Or go directly to <http://fractal.fractalartdesign.com>.

You will find yourself in an app, written in JavaScript and running in your browser, where you can endlessly explore the Mandelbrot set and an unlimited collection of variations. You can zoom, pan, change iterators, vary parameters, and select from almost two dozen themes or color maps. You can edit the colormaps and investigate smoothers, renderers, and other effects.

This is the tool that Greg uses to produce his art. He writes, "A painter uses paints; I use a CPU. I am giving you access to my paints! Please be aware that this is always a work in progress. It will change at a whim."

There is a lot of fractal art on the Web. This tool allows you to create your own.

Zooming in on a Mandelbrot or similar fractal by a factor approaching $2^{53}$ reaches the limit of IEEE double precision floating point. So, Greg's generator switches to *double double* precision. This doubles the fraction length and slows the rendering significantly. It does not increase the exponent, but these calculations need more precision, not more range. (See https://blogs.mathworks.com/cleve/2017/05/22/quadruple-precision-128-bit-floating-point-arithmetic.)

Greg uses his generator, and his skill, to create unique single-edition prints. To preserve their value, he does not reveal the parameters and never reprints. He displays his work at galleries and exhibitions around New England. Some of it is shown, and offered for sale, on his web site.

Here is a sample.

*Aquifer*

*Amphibian*

*Clematis*

Here is one of his wallpapers, digital images that are available as a free download for use as a desktop background or device lock screen.

*Wavelets*

Greg did not make the next two images -- I did, using his generator and nonstandard iterators. The color themes and shading are novel. This one is known on the Internet as the "Mandelbar" set because the iterator uses the complex conjugate, with a bar over the $z$,

$$z = \bar{z}^2 + c$$

instead of the traditional

$$z = z^2 + c$$

*Mandelbar*

And this amazing fractal, known as the "Burning Ship", is produced by zooming in on an iterator involving the complex absolute value.

$$z = |z|^2 + c$$

*Burning Ship*

To see much more of the Burning Ship fractal, check out this video, https://www.youtube.com/watch?v=CD9yNFmb2FE.

Get
the MATLAB code

Published with MATLAB® R2018a

I first saw this game years ago in the Exploratorium, a science museum in San Francisco. You are given a computer screen with two color patches and three sliders. One of the color patches is fixed. The sliders control the red, green and blue components of the other patch. Your task is to adjust the sliders until the two patches are the same color. It ain't easy.... read more >>

]]>I first saw this game years ago in the Exploratorium, a science museum in San Francisco. You are given a computer screen with two color patches and three sliders. One of the color patches is fixed. The sliders control the red, green and blue components of the other patch. Your task is to adjust the sliders until the two patches are the same color. It ain't easy.

My code for this game is available from the MATLAB® Central File Exchange.

The Exploratorium used random colors, but I'm going to pick the colors from a list of 143 named colors recognized by modern web browsers. These are wonderful names -- tomato, papaya whip, dark khaki, light sea green, burly wood, chocolate, moccasin, and many more. You can use the names in HTML and CSS code.

I found the names, and corresponding RGB triples, on a web site, color-names, provided by a San Francisco architecture and design firm, Dixon & Moe. Several other web sites, including Wikipedia, have similar lists.

The names are divided into ten groups -- reds, pinks, oranges, yellows, purples, greens, blues, browns, whites, and grays. There are two pairs of synonyms -- "fuchsia" and "magenta" are the same color and "aqua" and "cyan" are the same color. And there are two duplicates -- "light salmon" is listed in both the orange and red groups and "medium slate blue" is in both blue and purple. See if you can find them.

Here is the first screen of the game. I've already clicked on one of the tiles. The name, "blue violet", appears in the title and the color is used in a large patch that invites you to click it.

You can accept the invitation to click, or you can choose other tiles and see their names.

If you choose to click on the large patch, you get the second screen. Now the challenge begins. You must adjust the red, green and blue sliders until both patches have the same color.

If you succeed, you get rewarded with the gold star.

This example shows that "blue violet" has the RGB triple [138 43 226].

Some of the colors are easy to match. For "red", you obviously have to move R to its maximum, 255, and G and B to zero. But most are not so easy. When you need some help, click on the "score" toggle.

The score is the 1-norm of the difference between the RGB vectors of the two patches. That's the sum of the absolute values of the differences between the three components. In this example

left = [138 43 226]; right = [128 128 128]; score = norm(left - right,1)

score = 193

You want to drive this score to zero. The score is a monotone function of the sliders on either side of the optimum setting, so you don't have to worry about local minima.

The YIQ color space is now nearly obsolete, but I revealed my personal fondness for it in a blog post in 2016. YIQ was used by the NTSC color standard for analog TV in North America from 1954 until 2009. Y stands for *luminance*, I for *in phase*, and Q for *quadrature*. The Y component alone drove black-and-white TVs for half a century and is used today to convert RGB to gray scale.

If you have mastered the sliders in RGB space, click on the "rgb/yiq" toggle and try your hand there. Y varies between 0.0 and 1.0, while I and Q vary between obscure positive and negative limits. Again, the score is the 1-norm of the difference, now shown to three decimal places. I find it very difficult to steer to a match, even when I can see the score.

This screen shot was taken when the colors do not match and the score is far from 0.000. The YIQ coordinates for "blue velvet" are *not* [0.633 -0.447 0.302]. But is the variable patch on the right too bright or too dark?

Green is a special case. If you color something with the G component set to its maximum and the R and B components turned off, most people will say it's too bright. Here are two shades of pure green.

In the HTML naming scheme the dark one on the left, with G set to only 128, is called "green", while the light one on the right, with G set to 255, is "lime".

Here are two shades of pink. One is "deep pink", the other is "hot pink".

Both say, "click me". But the text in one is white and the other black. Why?

We want the text in dark colored patches to be white and the text in light colored patches to be black. The decision can be based on the first component of YIQ, *luminance*, which is computed by taking the dot product of a scaled RGB triple with

luma = [0.2989 0.5870 0.1140];

The RGB triple and the luminance for deep pink are

deep = [255 20 147]; deep_luminance = deep/255*luma'

deep_luminance = 0.4107

The corresponding quantities for hot pink are

hot = [255 105 180]; hot_luminance = hot/255*luma'

hot_luminance = 0.6211

The two luminance values fall on either side of 0.5, the luminance of gray. So deep pink gets white text while hot pink gets black text.

Today's homework: how many of our colors are dark and how many are light? Can you guess which are the most frequent? The answer surprised me. The ratio is nearly 2 to 1.

My code for this game is available from the MATLAB Central File Exchange.

```
type match_color_tail
```

%{ reds 1 indianred 205 92 92 2 lightcoral 240 128 128 3 salmon 250 128 114 4 darksalmon 233 150 122 5 lightsalmon 255 160 122 6 crimson 220 20 60 7 red 255 0 0 8 firebrick 178 34 34 9 darkred 139 0 0 pinks 10 pink 255 192 203 11 lightpink 255 182 193 12 hotpink 255 105 180 13 deeppink 255 20 147 14 mediumvioletred 199 21 133 15 palevioletred 219 112 147 oranges 16 lightsalmon 255 160 122 17 coral 255 127 80 18 tomato 255 99 71 19 orangered 255 69 0 20 darkorange 255 140 0 21 orange 255 165 0 yellows 22 gold 255 215 0 23 yellow 255 255 0 24 lightyellow 255 255 224 25 lemonchiffon 255 250 205 26 lightgoldenrodyellow 250 250 210 27 papayawhip 255 239 213 28 moccasin 255 228 181 29 peachpuff 255 218 185 30 palegoldenrod 238 232 170 31 khaki 240 230 140 32 darkkhaki 189 183 107 purples 33 lavender 230 230 250 34 thistle 216 191 216 35 plum 221 160 221 36 violet 238 130 238 37 orchid 218 112 214 38 fuchsia 255 0 255 39 magenta 255 0 255 40 mediumorchid 186 85 211 41 mediumpurple 147 112 219 42 rebeccapurple 102 51 153 43 blueviolet 138 43 226 44 darkviolet 148 0 211 45 darkorchid 153 50 204 46 darkmagenta 139 0 139 47 purple 128 0 128 48 indigo 75 0 130 49 slateblue 106 90 205 50 darkslateblue 72 61 139 51 mediumslateblue 123 104 238 greens 52 greenyellow 173 255 47 53 chartreuse 127 255 0 54 lawngreen 124 252 0 55 lime 0 255 0 56 limegreen 50 205 50 57 palegreen 152 251 152 58 lightgreen 144 238 144 59 mediumspringgreen 0 250 154 60 springgreen 0 255 127 61 mediumseagreen 60 179 113 62 seagreen 46 139 87 63 forestgreen 34 139 34 64 green 0 128 0 65 darkgreen 0 100 0 66 yellowgreen 154 205 50 67 olivedrab 107 142 35 68 olive 128 128 0 69 darkolivegreen 85 107 47 70 mediumaquamarine 102 205 170 71 darkseagreen 143 188 139 72 lightseagreen 32 178 170 73 darkcyan 0 139 139 74 teal 0 128 128 blues 75 aqua 0 255 255 76 cyan 0 255 255 77 lightcyan 224 255 255 78 paleturquoise 175 238 238 79 aquamarine 127 255 212 80 turquoise 64 224 208 81 mediumturquoise 72 209 204 82 darkturquoise 0 206 209 83 cadetblue 95 158 160 84 steelblue 70 130 180 85 lightsteelblue 176 196 222 86 powderblue 176 224 230 87 lightblue 173 216 230 88 skyblue 135 206 235 89 lightskyblue 135 206 250 90 deepskyblue 0 191 255 91 dodgerblue 30 144 255 92 cornflowerblue 100 149 237 93 mediumslateblue 123 104 238 94 royalblue 65 105 225 95 blue 0 0 255 96 mediumblue 0 0 205 97 darkblue 0 0 139 98 navy 0 0 128 99 midnightblue 25 25 112 browns 100 cornsilk 255 248 220 101 blanchedalmond 255 235 205 102 bisque 255 228 196 103 navajowhite 255 222 173 104 wheat 245 222 179 105 burlywood 222 184 135 106 tan 210 180 140 107 rosybrown 188 143 143 108 sandybrown 244 164 96 109 goldenrod 218 165 32 110 darkgoldenrod 184 134 11 111 peru 205 133 63 112 chocolate 210 105 30 113 saddlebrown 139 69 19 114 sienna 160 82 45 115 brown 165 42 42 116 maroon 128 0 0 whites 117 white 255 255 255 118 snow 255 250 250 119 honeydew 240 255 240 120 mintcream 245 255 250 121 azure 240 255 255 122 aliceblue 240 248 255 123 ghostwhite 248 248 255 124 whitesmoke 245 245 245 125 seashell 255 245 238 126 beige 245 245 220 127 oldlace 253 245 230 128 floralwhite 255 250 240 129 ivory 255 255 240 130 antiquewhite 250 235 215 131 linen 250 240 230 132 lavenderblush 255 240 245 133 mistyrose 255 228 225 grays 134 gainsboro 220 220 220 135 lightgray 211 211 211 136 silver 192 192 192 137 darkgray 169 169 169 138 gray 128 128 128 139 dimgray 105 105 105 140 lightslategray 119 136 153 141 slategray 112 128 144 142 darkslategray 47 79 79 143 black 0 0 0 %}

Get
the MATLAB code

Published with MATLAB® R2018a

Camille Jordan (1838-1922)... read more >>

]]>

Camille Jordan (1838-1922)

Photo credit: http://serge.mehl.free.fr/chrono/Jordan.html

(This column first appeared in the May 1994 issue of the MathWorks Newsletter.)

Matrices and differential equations are the key mathematical tools in MATLAB® and Simulink®. The Jordan Canonical Form is the key relationship between matrices and differential equations. So, why doesn't MATLAB use the JCF in any of its computations? In fact, until the Symbolic Math Toolbox came along, we didn't even have a function to compute the JCF.

The difficulty with the Jordan Canonical Form is that it is wildly sensitive to perturbation. Any kind of error -- uncertainty in experimental data, arithmetic roundoff error, linearization of nonlinear functions -- completely changes the JCF and, more importantly, the transformations that generate it.

Let's start with a system of $n$ linear, constant coefficient ordinary differential equations.

$$\dot{x} = Ax$$

Here $x(t)$ is a vector-valued function of $t$, the dot denotes differentiation with respect to $t$, and $A$ is a given $n$ -by- $n$ matrix which is independent of $t$. If you've never heard of -- or want to forget -- the Jordan Canonical Form, then you have to hope there is a matrix $V$ which diagonalizes $A$, that is $V^{-1} A V$ is a diagonal matrix $\Lambda$. Then the change of variables,

$$x = Vy$$

reduces the $n$ -by- $n$ problem to $n$ instances of the 1-by-1 case,

$$\dot{y}_k = \lambda_k y_k$$

The elements of the diagonal matrix $\Lambda$ are the *eigenvalues* of $A$ and the columns of $V$ are the *eigenvectors* of $A$. The $k$-th component of the solution $y(t)$ is the exponential function of $t$ determined by the eigenvalue $\lambda_k$, and the initial value $y_k(0)$.

$$ y_k(t) = y_k(0) e^{\lambda_k t} $$

The components of the solution $x(t)$ of the original problem are then linear combinations of the exponential functions determined by the eigenvalues and eigenvectors and the initial values.

But this fantasy about diagonalizing matrices leaves out some very important differential equations, like the single, second order

$$\ddot{u} = 0$$

The solution, of course, is a straight line -- a linear polynomial in $t$. The solution does *not* involve any exponentials. This equation can be converted to a 2-by-2, first-order system by taking $x$ to be the vector with components

$$x_1 = u$$

$$x_2 = \dot{u}$$

Then the differential equation is

$$\dot{x_1} = x_2$$

$$\dot{x_2} = 0$$

This can be written in matrix form

$$\dot{x} = J x$$

where $J$ is the 2-by-2 matrix generated by the MATLAB statement

J = [ 0 1; 0 0 ]

J = 0 1 0 0

This matrix *cannot* be diagonalized; no change of variable which preserves the differential equation will produce a zero in the upper right hand corner. That had better be the case, because the differential equation has polynomial, not exponential, solutions. The matrix `J` is the simplest example of a nondiagonal Jordan Canonical Form. The eigenvalues -- in this case, a pair of zeros - are on the diagonal, but that crucial `1` off the diagonal couples the two differential equations.

A different change of variable leads to a matrix whose structure is not quite so obvious. Let

$$x_1 = u + \dot{u}$$

$$x_2 = u - \dot{u}$$

Then the equation becomes

$$\dot{x} = Ax$$

where $A$ is generated by the MATLAB statement

A = [ 1/2 -1/2; 1/2 -1/2 ]

A = 0.5000 -0.5000 0.5000 -0.5000

This new change of variables really hasn't altered things very much. The Jordan Canonical Form of `A` is the matrix `J` from our first system. The differential equation still has solutions which are linear polynomials in $t$; it can't be transformed into a pair of decoupled equations with exponential solutions.

But now, let's perturb the matrix A very slightly to

A = [ 0.5 -0.49999999; 0.49999999 -0.5 ]

A = 0.500000000000000 -0.499999990000000 0.499999990000000 -0.500000000000000

The eigenvalues are changed from a pair of zeros to

lambda = eig(A)

lambda = 1.0e-04 * 1.0000 -1.0000

A change of $10^{-8}$ in the matrix elements changes the eigenvalues by $10^{-4}$. Now this perturbed matrix has distinct eigenvalues, and can be diagonalized. The eigenvector matrix is

[V,~] = eig(A)

V = 0.707177488329075 0.707036066972953 0.707036066972953 0.707177488329075

This eigenvector matrix `V` defines a change of variables that transforms the system into a decoupled system whose solutions involve `exp(t/10000)` and `exp(-t/10000)`. But we know the solution should be very close to a linear function of $t$. It is probably a lousy idea to try to represent the solution in terms of these two exponentials.

This is reflected in the fact that `V` is badly conditioned.

cond(V) = 1.0000e+04

This matrix `A` is an example of why we can't use the Jordan Canonical Form for practical computation. Technically, `A` has a full set of linearly independent eigenvectors and a diagonal JCF. But any attempt to use this fact will be plagued by the condition number of `V` and a resulting loss of accuracy. On the other hand, if we try to say that `A` does not have a full set of eigenvectors and use a nondiagonal JCF, we are also making possibly significant errors.

Another example can be generated starting with MATLAB's test matrix,

A = gallery(5)

A = -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

This matrix was constructed so that its characteristic polynomial is

$$\lambda^5 = 0$$

All five eigenvalues are zero. The Jordan Canonical Form can be obtained from the Symbolic Math Toolbox which, in this case, does exact integer calculations with *no* round off error.

J = jordan(A)

J = 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0

With a well-conditioned change of variables, the equation

$$ \dot{x} = Ax$$

essentially becomes

$$ \frac{d^4u}{dt^4} = 0 $$

All solutions are cubic polynomials in $t$. So far, so good.

But now change the second diagonal element of `A` from -69 to

A(2,2) = -70

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

I want to regard that as a significant change; I am not willing to approximate this new matrix by the unaltered `gallery(5)`. The characteristic polynomial is now

factor(charpoly(A,'lambda'))

$$(\lambda-1)(\lambda^4 + 2\lambda^3 - 67\lambda^2 + 234\lambda - 168)$$

One of the eigenvalues is $1$. The other four are the distinct, irrational numbers which are the roots of an irreducible quartic. The Jordan Canonical Form is diagonal. We could ask the Symbolic Math Toolbox to compute it, but the result is not very usable.

Instead we can simply ask for

[V, lambda] = eig(A); lambda = diag(lambda)

lambda = -10.5726 + 0.0000i 3.7953 + 1.3331i 3.7953 - 1.3331i 1.0000 + 0.0000i 0.9820 + 0.0000i

The form of the solution to the differential equation has changed dramatically, from a cubic polynomial in $t$, to a function involving terms like $\exp(-10.5726 t)$ and $\exp(3.7953 t) \sin(1.3331 t)$. Moreover, the change of variables required to obtain this representation is badly conditioned.

cond(V) = 2.3890e+04

The fundamental difficulty is with matrices, like this last `A`, which are close to, but not exactly equal to, matrices with nondiagonal JCFs. If we try to diagonalize such matrices, then errors are magnified by the condition number of the eigenvector matrix, which can be arbitrarily large. On the other hand , to use the "nearby" nondiagonal JCF may also represent an unacceptably large error. We're damned if we do and damned if we don't.

The numerically reliable approach to all this is to avoid the Jordan Canonical Form altogether. The portions of MATLAB and Simulink that deal with matrix differential equations use other solution methods, including the Schur Form (which is triangular, instead of bidiagonal) and the matrix exponential (which is *not* computed from eigenvalues and eigenvectors).

I can't resist finishing up with a personal story. Prof. G. W. "Pete" Stewart of the University of Maryland introduces speakers at conferences by composing original limericks for them. Several years ago, he honored me with:

Said Professor Cleve Moler one day, I'd like to get e to the A. It's as easy as heck, To get e to the x, But I can't get e to the J.

Thanks, Pete.

Get
the MATLAB code

Published with MATLAB® R2017a

The ACM Special Interest Group on Programming Languages, SIGPLAN, expects to hold the fourth in a series of conferences on the History of Programming Languages in 2020, see HOPL-IV. The first drafts of papers are to be submitted by August 2018. That long lead time gives me the opportunity to write a detailed history of MATLAB. I plan to write the paper in sections, which I'll post in this blog as they are available. This is the seventh, and final, installment.... read more >>

]]>The ACM Special Interest Group on Programming Languages, SIGPLAN, expects to hold the fourth in a series of conferences on the History of Programming Languages in 2020, see HOPL-IV. The first drafts of papers are to be submitted by August 2018. That long lead time gives me the opportunity to write a detailed history of MATLAB. I plan to write the paper in sections, which I'll post in this blog as they are available. This is the seventh, and final, installment.

MATLAB has evolved over the past 35 years into a rich technical computing environment. Here is an overview of some of the tools available.

The MATLAB desktop was introduced in 2000. Here is the default layout of the desktop. This snapshot was taken while I was working on my previous post. Four panels are visible, the current folder viewer on the left, the workspace viewer on the right, the editor/debugger in the top center, and the traditional command window in the bottom center. A file viewer and a command history window are not in the default layout, but can be included in personalized layouts.

Any of the panels can be closed, or undocked into a standalone window. I have two screens available on the desk in my office. I usually put the command window on one screen and the editor on the other.

Here are the contents of the workspace when I execute the source script for my previous post. If I click on any of the variables, it will be opened in an appropriate viewer

When I work on this blog, I spend most of my time using this MATLAB editor. The source for each post is an executable script, with descriptive text in comments and section headings in comments beginning with double percent signs, `%%`. The final production of a post is initiated with the "Publish" tab. I have pretty much abandoned my old Unix text editor, "vi", after many years of faithful service.

The Live Editor was introduced in 2016 and is still under active development. MATLAB input, output and graphics are combined in a single interactive document. The document may be exported to HTML, PDF, or LaTeX.

Here are a few examples. The table generated by the small gradebook from my previous post is nicely formatted.

Expressions from the Symbolic Toolbox are typeset.

And graphical output may be included in the document.

The Parallel Computing Toolbox (PCT) was introduced at the Super Computing conference in 2004. The next year, at SC05, Bill Gates gave the keynote talk, using MATLAB to demonstrate Microsoft's entry into High Performance Computing.

The toolbox supports coarse grained, distributed memory parallelism by running many MATLAB workers on many machines in a cluster, or on many cores in a single machine. Originally MATLAB used MPI for the message passing, but in more recent versions of the toolbox our own utilities have replaced MPI.

By far the most popular feature of the PCT is the parallel for loop command, `parfor`.

Support for Graphics Processing Units was added to the Parallel Computing Toolbox in 2010. Eight years later, in release R2018a, the `gpuarray` has grown to have 385 associated methods, including all of the familiar matrix computations, `lu`, `eig`, `svd` and `mldivide` (backslash).

Much of modern MATLAB's power derives from the toolboxes available for specialized applications. In release 2018a there are 63 of them. Here is the list.

**Parallel Computing**

- MATLAB Distributed Computing Server
- Parallel Computing Toolbox

**Math, Statistics, and Optimization**

- Curve Fitting Toolbox
- Global Optimization Toolbox
- Model-Based Calibration Toolbox
- Neural Network Toolbox
- Optimization Toolbox
- Partial Differential Equation Toolbox
- Statistics and Machine Learning Toolbox
- Symbolic Math Toolbox
- Text Analytics Toolbox

**Control Systems**

- Aerospace Toolbox
- Automated Driving System Toolbox
- Control System Toolbox
- Fuzzy Logic Toolbox
- Model Predictive Control Toolbox
- Robotics System Toolbox
- Robust Control Toolbox
- System Identification Toolbox

**Signal Processing and Wireless Communications**

- Antenna Toolbox
- Audio System Toolbox
- Communications System Toolbox
- DSP System Toolbox
- LTE HDL Toolbox
- LTE System Toolbox
- Phased Array System Toolbox
- RF Toolbox
- Signal Processing Toolbox
- Wavelet Toolbox
- WLAN System Toolbox

**Image Processing and Computer Vision**

- Automated Driving System Toolbox
- Computer Vision System Toolbox
- Image Acquisition Toolbox
- Image Processing Toolbox
- Mapping Toolbox
- Vision HDL Toolbox

**Test and Measurement**

- Data Acquisition Toolbox
- Image Acquisition Toolbox
- Instrument Control Toolbox
- OPC Toolbox
- Vehicle Network Toolbox

**Computational Finance**

- Database Toolbox
- Datafeed Toolbox
- Econometrics Toolbox
- Financial Instruments Toolbox
- Financial Toolbox
- Risk Management Toolbox
- Spreadsheet Link
- Trading Toolbox

**Computational Biology**

- Bioinformatics Toolbox
- SimBiology

**Code Generation**

- Filter Design HDL Coder
- Fixed-Point Designer
- GPU Coder
- HDL Coder
- HDL Verifier
- MATLAB Coder
- Vision HDL Toolbox

**Application Deployment**

- MATLAB Compiler
- MATLAB Compiler SDK
- Spreadsheet Link

**Database Access and Reporting**

- Database Toolbox
- MATLAB Report Generator

Get
the MATLAB code

Published with MATLAB® R2018a