This week’s post is by Reece Teramoto.... read more >>

]]>This week’s post is by Reece Teramoto.

MathWorks recently took over sponsorship of the Math Modeling Challenge, a contest for high school juniors and seniors in the U.S organized by SIAM – Society for Industrial and Applied Mathematics. https://m3challenge.siam.org. Students work in teams to solve real-world problems that professional mathematicians face. Teams are given 14 hours to complete a three-part problem using math modeling and put together a paper showing their results. Teams can choose to use any software or language (or none) to do their mathematical computations.

We decided to put together a set of sample solutions to this past year’s problem done with MATLAB.

This past year’s challenge focused on the problem of wasted food. Students were challenged to perform computations to analyze the amount of wasted food and apply the results to come up with solutions of feeding the millions of food-insecure individuals in the country. You can view the full problem here and the winning solutions here.

In this post, I’ll be walking through a brief overview of each part of the problem as well as a summary of a sample solution put together in MATLAB. The corresponding live script for each example will demonstrate how to do the following in MATLAB:

- Import data
- Process and model data
- Visualize the results (charts, graphs, curve fitting, plots)
- Make meaningful conclusions from the results

While all three parts of the problem demonstrate the steps above, the three parts will each focus on a different tool in MATLAB. We will focus on the Import Tool, Curve Fitting app, and Mapping Toolbox.

Before reading through the rest of this post, please download the full MATLAB examples from File Exchange so that you may refer to them.

- Just Eat It! (Import Tool)
- Food Foolish? (Curve Fitting app)
- Hunger Game Plan? (Mapping Toolbox)
- Want to participate in the challenge?

You can find the full example for this problem by opening the file `Just_Eat_It.mlx` from the archive linked to at the beginning of this blog post.

For this part of the problem, students are to create a mathematical model that a state could use to determine if it could feed its food-insecure population using the wasted food generated in that state. Students are provided with sample data for food waste and insecurity in the state of Texas.

This is one way to approach this problem, employing technical computing skills:

- Calculate how many dollars of food the state consumers waste.
- Calculate how many dollars are needed to feed all the food-insecure individuals in the state.
- Compare the two values to determine how much of the state’s food-insecure population can be fed by using the wasted food of the state’s consumers.

First, we can use the MATLAB Import Tool to import the provided spreadsheet data and save everything to tables. This import can be performed with a few clicks:

Once we have our data in MATLAB, we can use many built-in MATLAB functions to organize and summarize our data.

From the data, we can calculate the dollar value of wasted food in Texas for different categories, which we can visualize:

bar(foodWasteTable.TypeofFood, foodWasteTable.WastedDollars); ytickformat('usd') xlabel('Type of Food') ylabel('Amount Wasted') title('Dollars of Wasted Food in Texas (2016)') % change the y-axis to be in the form of n x 10^6 (millions) ax = gca; ax.YAxis.Exponent = 6;

% sort from smallest to largest food waste foodWasteTable = sortrows(foodWasteTable,'WastedDollars'); % combine the 3 smallest categories and label them as "other" p = pie([sum(foodWasteTable.WastedDollars(1:3)); foodWasteTable.WastedDollars(4:end)]); legend(["other"; foodWasteTable.TypeofFood(4:end)], 'Location', 'eastoutside') title('Categories of Wasted Food')

As one might expect, the two largest types of wasted food are meat and produce. After further calculations, we find that the value of wasted food in Texas is about $2.3 billion annually.

The second step is finding the dollars are needed to feed all the food-insecure individuals. Feeding America performs this calculation here. We can easily translate this into MATLAB to find that the annual cost to feed food insecure individuals in Texas is about $2.2 billion.

For the final comparison, and based on our assumptions noted in the file `Just_Eat_It.mlx`, we determined that we can feed all the food-insecure individuals in Texas that by repurposing wasted food in the state.

You can find the full example for this problem by opening the file `Food_Foolish.mlx` from the archive linked to at the beginning of this blog post.

Personal choices when it comes to food consumption primarily occur at the grocery store, school cafeteria, restaurants, and at home. For this part of the problem, students are to create a mathematical model that can be used to determine the amount of food waste a household generates in a year based on their traits and habits. The model should then be applied to the following households:* *

- Single parent with a toddler, annual income of $20,500
- Family of four (two parents, two teenage children), annual income of $135,000
- Elderly couple, living on retirement, annual income of $55,000
- Single 23-year-old, annual income of $45,000

We will use the Curve Fitting App in MATLAB to find the relationship between household income and wasted food. Our approach is as follows:

- Gather data for household incomes and wasted food. We can think of each data point as a point on an XY-plane, where the X-axis is household income and the Y-axis is wasted food.
- Use the Curve Fitting App to generate a function that models the trend.
- Visualize this relationship.
- Evaluate the curve at the income levels of the sample family data.

The U.S. Bureau of Labor Statistics has collected data on annual food expenditures for individuals of varying income. You can find the data in this spreadsheet. We can import this data into MATLAB using the Import Tool and store the data in a table.

Once we have the data points, we can use the Curve Fitting App to generate the curve fit. The app makes it easy to quickly try out different fits of a curve to the data points and see the immediate result of each fit:

Once we have our fitted curve, we can pass in the sample family income data:

sampleIncome = [20500 135000 55000 45000]'; sampleWaste = foodWasteCurve(sampleIncome); families = {'Single parent with a toddler', 'Family of four (two parents, two teenage children)', ... 'Elderly couple, living on retirement','Single 23-year-old'}; solution = table(sampleIncome,sampleWaste, 'RowNames', families)

We can also plot both a bar graph and line graph of the trend of food waste for various incomes. Then, we can plot points for the 4 sample families given in the problem statement. This graph makes it easy to visualize the relationship between household income and annual dollars of wasted food. Since we assumed that food waste depends *only* on household income, we have a solution for this part of the problem.

% bar graph bar(incomeRange, wasteDemo) hold on % line graph plot(incomeRange, wasteDemo, 'LineWidth', 3, 'Color', 'green') % plot points of the 4 families from the problem statement scatter(solution.sampleIncome, solution.sampleWaste, 'Marker', '*', 'LineWidth', 5, 'MarkerEdgeColor', 'red') hold off xlabel('Income (USD)') ylabel('Annual Food Waste (USD)') title ('Food Waste by Income')

** **You can find the full example for this problem by opening the file `Hunger_Game_Plan.mlx` from the archive linked to at the File Exchange <link>.

For this problem, students were to use mathematical modeling to provide insight on a strategy that a local community could use to repurpose wasted food with minimal cost.

Our approach was as follows:

- Use Mapping Toolbox to display location of two dozen grocery stores, including all the Stop and Shops, in Boston.
- Display the location of The Greater Boston Food Bank.
- Find the shortest path that starts at the food bank, visits all the grocery stores to collect otherwise-wasted produce from the past week, and ends at the food bank.
- Calculate the trucking costs and dollars of saved produce from the past week.

Since Stop and Shop is the most popular grocery store chain in New England and The Greater Boston Food Bank is an organization that MathWorks supports, our approach for this problem was to find the shortest path that a food collection truck could start at The Greater Boston Food Bank, visit each Stop and Shop and a few other stores in the Boston area to collect otherwise-wasted produce, and return to The Greater Boston Food Bank. We want the shortest path since this would minimize the cost of operating the truck.

This particular scenario is known in computer science as the Travelling Salesman Problem (TSP), in which a salesman starts from their home city and wants to visit multiple cities and return home in the shortest possible path.

Mapping Toolbox lets us plot points on a map given their latitude and longitude values. We start by displaying all the Stop and Shop locations in the Boston area (identified by a green marker). The coordinates were obtained from Google Maps.

This TSP with only 13 locations is not much of a challenge, so we add 11 more locations (shown in dark red) chosen at random, as well as The Greater Boston Food Bank (indicated in blue).

wmmarker(Lats(1:13), Lons(1:13), 'Icon', 'stopandshoplogo.jpg') wmmarker(Lats(14:24), Lons(14:24), 'color', darkred) wmmarker(Lats(25), Lons(25), 'color', ‘blue’)

The previous Cleve’s Corner post described a near-optimal solution of the Travelling Salesman Problem.. We can pass in the coordinates to traveler to obtain the shortest path between them. We can view this path using the Mapping Toolbox as well:

Lats = Lats(path); Lons = Lons(path); wmline(Lats, Lons)

Although the path returned does not consider the exact path of the roads, we still get a good estimation of the shortest path and the distance traveled. The expense of operating an 18-wheeler is a little over \$1 per mile, we can calculate the total cost of gas for this trip—around \$30.

Additionally, we can use some data from the National Resources Defense Council to calculate that the average grocery store in the U.S. wastes about \$7500 per week in fruits and vegetables. This means that by completing this round trip, we will salvage almost \$200,000 worth of otherwise-wasted fruits and vegetables from just these locations in the Boston area from the past week. Pretty good!

Hopefully this has given you a glimpse at how easy MATLAB makes it to analyze, visualize, and perform computations on data. To view the full MATLAB solutions for the above problems, please download the files from File Exchange.

Additionally, if you know of any high school juniors and seniors who would be interested in the challenge, encourage them to participate! The challenge is free to enter and there are $100,000 in scholarships at stake. The challenge is done completely online and there are plenty of reference materials to help teams prepare. For more information, please visit the following link:

]]>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