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

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.... 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 sixth such installment and the second of a three-part post about modern MATLAB.

MATLAB is not just a Matrix Laboratory any more. It has evolved over the past 35 years into a rich technical computing environment. Here are some samples of the features now available.

Handle Graphics has evolved over the years into a rich, tree-structured, object-oriented environment for producing publication quality plots and graphs. Let's see just one example.

sinc = @(x) sin(x)./x; x = (-5:.03:5)*pi; px = plot(x,sinc( x));

This is not ready yet for publication. We need to fine tune the annotation.

ax = px.Parent; ax.XLim = pi*[-5 5]; ax.XTick = pi*(-4:2:4); ax.XTickLabel = {'-4\pi','-2\pi','0','2\pi','4\pi'}; ax.YLim = [-0.3 1.1]; ax.Title.String = 'sinc(x)';

This plot is now ready to be included in a textbook or expository paper. We created two Handle Graphics objects, `px` and `ax`. `px` is a `line` object with 38 properties. `ax` is an `axes` object with 137 properties. Setting a few of these properties to values other than their defaults improves the scaling and annotation.

Structures, and associated "dot notation", were introduced in 1996. As we just saw, Handle Graphics now uses this notation. For another example, let's create a grade book for a small class.

Math101.name = {'Alice Jones';'Bob Smith';'Charlie Brown'}; Math101.grade = {'A-'; 'B+'; 'C'}; Math101.year = [4; 2; 3]

Math101 = struct with fields: name: {3×1 cell} grade: {3×1 cell} year: [3×1 double]

To call the roll, we need the list of names.

disp(Math101.name)

'Alice Jones' 'Bob Smith' 'Charlie Brown'

Changing Charlie's grade involves both structure and array notation.

```
Math101.grade(3) = {'W'};
disp(Math101.grade)
```

'A-' 'B+' 'W'

Tables, which were introduced in 2013, allow us to print the grade book in a more readable form.

T = struct2table(Math101)

T = 3×3 table name grade year _______________ _____ __________ 'Alice Jones' 'A-' 4.0000e+00 'Bob Smith' 'B+' 2.0000e+00 'Charlie Brown' 'W' 3.0000e+00

Major enhancements to the MATLAB object-oriented programming capabilities were made in 2008. Creating classes can simplify programming tasks that involve specialized data structures or large numbers of functions that interact with particular kinds of data. MATLAB classes support function and operator overloading, controlled access to properties and methods, reference and value semantics, and events and listeners.

Handle Graphics is one large, complex example of the object-oriented approach to MATLAB programming.

The `datetime` object, introduced in 2014, is another example. For many years, computations involving dates and time were done with using a collection of a dozen functions including `clock`, `datenum` and `datestr`. The `datetime` object combines and extends the capabilities of all those functions.

Here is an unexpected use of `datetime`. How long is a *microcentury*? Let's start with the date that I am publishing this post.

```
t1 = datetime('today')
```

t1 = datetime 27-Apr-2018

Add 100 years.

[y,m,d] = ymd(t1); t2 = datetime(y+100,m,d)

t2 = datetime 27-Apr-2118

A century is the difference between these two times. This difference is a `duration` object, expressed in hours:minutes:seconds.

c = t2 - t1

c = duration 876576:00:00

Here is a microcentury, in the same units.

microc = 1.e-6*c

microc = duration 00:52:35

So, a microcentury is 52 minutes and 35 seconds. The optimum time for an academic lecture.

But wait -- some centuries are longer than most, because they have four leap years, not the usual three. The year 2100, which we have just included in our century, will not be a leap year, On the other hand, the year 2000, which is divisible by 400, was a leap year. We need more precision.

microcentury_short = duration(microc, 'Format','hh:mm:ss.SSS')

microcentury_short = duration 00:52:35.673

Subtract, rather than add, 100 years.

t0 = datetime(y-100,m,d); microcentury_long = 1.e-6*duration(t1-t0, 'Format','hh:mm:ss.SSS')

microcentury_long = duration 00:52:35.760

The numerical solution of ordinary differential equations has been a vital part of MATLAB since its commercial beginning. ODEs are the heart of Simulink®, MATLAB's companion product.

The Van der Pol oscillator is a classical ODE example.

$$ \ddot{y} = \mu (1-y^2) \dot{y} - y$$

The parameter $\mu$ is the strength of the nonlinear damping term. When $\mu = 0$, we have the basic harmonic oscillator.

The MATLAB code expresses the oscillator as a pair of first order equations.

mu = 5; vdp = @(t,y) [0 1; -1 mu*(1-y(1)^2)]*y; tspan = [0 30]; y0 = [0 0.01]';

[t,y] = ode23s(vdp,tspan,y0);

plot(t,y,'.-') legend({'y','dy/dt'}) xlabel('t')

The Van der Pol oscillator, with the parameter `mu` set to a value of 5, is a mildly *stiff* differential equation. In anticipation, I used the `ode23s` solver; the `'s'` in the name indicates it is for stiff equations. In the plot you can see some clustering of steps where the solution is varying rapidly. A nonstiff solver would have taken many more steps.

A stiff ode solver uses an *implicit* method requiring the solution of a set of simultaneous linear equations at each step. The iconic MATLAB backslash operator is quietly at work here.

Support for *graph theory* was added in 2015. In this context, a *graph* is a set of nodes or vertices, together with a set of edges or lines between the nodes. There are two kinds of graphs. A `graph` is *undirected* if an edge between two nodes means that the nodes are connected to each other. A `digraph` is *directed* if the connections are one-way. Edges can have different weights, or all the weights can be the same.

G = digraph(1,2:5); G = addedge(G,2,6:15); G = addedge(G,15,16:20)

G = digraph with properties: Edges: [19×1 table] Nodes: [20×0 table]

Graphs have no inherent geometric structure. The tricky part of plotting a graph is locating the nodes in two- or three-dimensional space. The graph plot methods offer several different layouts.

plot(G,'layout','force','nodecolor',burnt_orange) no_ticks

The Image Processing Toolbox has been available since the 1990's. It is now at Version 10.2. The command

```
help images
```

lists over 500 functions. As expected, the toolbox makes extensive use of array functionality.

This image, "Psychedelic Pluto", combines image processing and numerical linear algebra. The NASA press release says "New Horizons scientists made this false color image of Pluto using a technique called **principal component analysis** (emphasis mine) to highlight the many subtle color differences between Pluto's distinct regions." The image was presented by Will Grundy of the New Horizons’ surface composition team on Nov. 9, 2015.

*Image Credit: NASA/JHUAPL/SwRI*

Images of Pluto are of particular interest because MATLAB and Simulink were used extensively in both navigation and the flight control system for the New Horizons flyby mission.

The Symbolic Math Toolbox has also been available since the 1990's. It is now at version 8.1. The command

```
help symbolic
```

lists almost 500 functions.

The toolbox provides functions for solving, plotting, and manipulating symbolic math equations. The operations include differentiation, integration, simplification, transforms, and equation solving. Computations can be performed either analytically or using variable-precision arithmetic.

The toolbox has recently been integrated with the Live Editor, which I will discuss in my next blog post. This example does not yet use the Live Editor.

The statement

```
x = sym('x')
```

x = x

creates a symbolic variable that just represents itself. The variable does not have a numeric value. The statement

f = 1/(4*cos(x)+5)

f = 1/(4*cos(x) + 5)

creates a symbolic function of that variable.

Without the Live Editor, the closest we can come to LaTeX style typeset mathematics in the command window is

pretty(f)

1 ------------ 4 cos(x) + 5

Let's differentiate `f` twice.

f2 = diff(f,2)

f2 = (4*cos(x))/(4*cos(x) + 5)^2 + (32*sin(x)^2)/(4*cos(x) + 5)^3

And then integrate the result twice. Do we get back to where we started?

g = int(int(f2))

g = -8/(tan(x/2)^2 + 9)

No, `f` and `g` are not the same expression. Let's plot both of them.

subplot(211) ezplot(f) subplot(212) ezplot(g)

The two plots have the same shape, but the y axes are not the same. Let's compute the difference; this is the "error" made by differentiating twice and then integrating twice.

e = f - g

e = 8/(tan(x/2)^2 + 9) + 1/(4*cos(x) + 5)

The error isn't zero. What is it? First a plot.

half_fig ezplot(e)

Now let's do the most difficult computation of this entire example.

s = simplify(e)

s = 1

So, `e` is an elaborate way of writing the constant function, 1. That's the "constant of integration", the infamous +C that we always forget to include in indefinite integration.

Here's your homework for today. Where did the limits on the y axis for `ezplot(e)` come from?

Get
the MATLAB code

Published with MATLAB® R2018a

Today is Friday, the 13th. In many parts of the world, today is regarded as *unlucky*. But I want to revisit an old question: is today *unlikely*? What are the chances that the 13th of any month falls on a Friday? Computing the answer makes use of a new MATLAB® feature, the `datetime` method.... read more >>

Today is Friday, the 13th. In many parts of the world, today is regarded as *unlucky*. But I want to revisit an old question: is today *unlikely*? What are the chances that the 13th of any month falls on a Friday? Computing the answer makes use of a new MATLAB® feature, the `datetime` method.

I wrote about Friday the 13th in a blog post several years ago, Cleve's Corner, Friday the 13th. I am reusing a portion of that post today.

Leap years make our calendar a nontrivial mathematical object. The leap year rule can be implemented by this anonymous function.

leapyear = @(y) mod(y,4)==0 & mod(y,100)~=0 | mod(y,400)==0;

This says that leap years happen every four years, except for the turn of centuries not divisible by 400. Try a few year numbers.

y = [2018 2020 2000 2100]'; isleap = [y leapyear(y)]

isleap = 2018 0 2020 1 2000 1 2100 0

So, this year is not a leap year, the next one is in 2020. The turn of the century 2000 was a leap year, but 2100 will not be.

The leap year rule implies that our calendar has a period of 400 years. The calendar for years 2000 to 2399 will be reused for 2400 to 2799. In a 400 year period, there are 97 leap years, 4800 months, 20871 weeks, and 146097 days. So the average number of days in a calendar year is not 365.25, but

dpy = 365+97/400

dpy = 365.2425

We can compute the probability that the 13th of a month is a Friday by counting how many times that happens in 4800 months. The correct probability is then that count divided by 4800. Since 4800 is not a multiple of 7, the probability does not reduce to 1/7.

For many years, computations involving dates and time have been done using a whole bunch of functions.

`clock``datenum``datevec``datestr``dateticks``now``weekday`

Those functions have served us well, but they have some significant shortcomings. They are not particularly easy to use. The display formatting is limited. The floating point precision of `datenum` effectively limits the calculation of durations. There are no immediate plans to retire `datenum` and its friends, but something preferable is now available.

The `datetime` object, introduced in R2014b, combines and extends these functions into one.

`datetime`

The documentation for `datetime` is available online, or, from within MATLAB, with the command

doc datetime

For example, here is the date and time when I am publishing this post.

```
d = datetime('now')
```

d = datetime 13-Apr-2018 21:21:19

I can specify a custom display format with

d = datetime(d,'Format','eeee, MMMM d, yyyy HH:mm:ss')

d = datetime Friday, April 13, 2018 21:21:19

The available display formats include support for an international standard, ISO 8601.

`datetime`'s are much more precise than `datenum`'s can ever hope to be. Try this

hms = datetime(d,'Format','HH:mm:ss.SSS')

hms = datetime 21:21:19.502

We see that `d` is carrying fractional seconds to at least three places beyond the decimal point. How about nanoseconds?

secperday = 24*60*60 nano = 1.e-9/secperday hms = datetime(d,'Format','HH:mm:ss.SSSSSSSSS') hmsplusnano = hms + nano

secperday = 86400 nano = 1.1574e-14 hms = datetime 21:21:19.502000000 hmsplusnano = datetime 21:21:19.502000001

A nanosecond is four orders of magnitude smaller than roundoff error of `datenum`'s in this era. Compare

nano epsofdatenum = eps(datenum(d))

nano = 1.1574e-14 epsofdatenum = 1.1642e-10

I am pleased to reveal some details about the arithmetic in `datetime`. The object is using a form of quadruple precision known as double-double. In fact, the second half of the 112-bit format is stored as the imaginary part of the `data` field, as you see in the following.

disp(struct(hmsplusnano))

Warning: Calling STRUCT on an object prevents the object from hiding its implementation details and should thus be avoided. Use DISP or DISPLAY to see the visible public details of an object. See 'help struct' for more information. UTCZoneID: 'UTC' UTCLeapSecsZoneID: 'UTCLeapSeconds' ISO8601Format: 'uuuu-MM-dd'T'HH:mm:ss.SSS'Z'' epochDN: 719529 MonthsOfYear: [1×1 struct] DaysOfWeek: [1×1 struct] dateFields: [1×1 struct] noConstructorParamsSupplied: [1×1 struct] data: 1.5237e+12 + 1.0000e-06i fmt: 'HH:mm:ss.SSSSSSSSS' tz: '' isDateOnly: 0

The summary description in the `datetime` documentation describes some important features of the object, which has nearly 100 methods defined.

`datetime` arrays represent points in time using the proleptic
ISO calendar. `datetime` values have flexible display formats
up to nanosecond precision and can account for time zones,
daylight saving time, and leap seconds.

I had to look up the definition of "proleptic". In this context it means the calendar can be extended backwards in time to apply to dates even before it was adopted.

Let's get on with Friday the 13th. Like most things in MATLAB `datetime` works with vectors. The statements

y = 2000; m = 1:4800; t = 13; v = datetime(y,m,t);

produce a row vector of 4800 `datetime`'s for the 13th of every month. The first few and last few are

fprintf(' %s',v(1:4)) fprintf(' ...\n') fprintf(' %s',v(end-3:end))

13-Jan-2000 13-Feb-2000 13-Mar-2000 13-Apr-2000 ... 13-Sep-2399 13-Oct-2399 13-Nov-2399 13-Dec-2399

The statement

w = weekday(v);

produces a row vector of 4800 flints between 1 (for Sunday) and 7 (for Saturday). The first few and last few are

fprintf('%3d',w(1:4)) fprintf(' ...') fprintf('%3d',w(end-3:end))

5 1 2 5 ... 2 4 7 2

Now to get the counts of weekdays, all we need is

c = histcounts(w)

c = 687 685 685 687 684 688 684

There you have it. The count for Friday, 688, is higher than for any other day of the week. The 13th of any month is slightly more likely to fall on a Friday.

The probabilities are

p = c/4800

p = 0.1431 0.1427 0.1427 0.1431 0.1425 0.1433 0.1425

Compare these values with their mean, `1/7 = 0.1429`. Four probabilities are below the mean, three are above.

Let's pair the counts with the days of the week, using a `categorical` variable, introduced in R2013b.

c = categorical(w,1:7,{'Sun' 'Mon' 'Tue' 'Wed' 'Thu' 'Fri' 'Sat'}); summary(c)

Sun Mon Tue Wed Thu Fri Sat 687 685 685 687 684 688 684

When we plot a histogram the appropriate labels are provided on the x-axis. Friday the 13th is the winner.

histogram(c) set(gca,'ylim',[680 690]) title('Weekday counts, 400 years')

Thanks to Peter Perkins at MathWorks. I learned a lot from him while working on this post.

Get
the MATLAB code

Published with MATLAB® R2018a

Let me tell you about a beautiful, fractal curve, the Dragon Curve. Download my new `dragon` program from the File Exchange and follow along.... read more >>

Let me tell you about a beautiful, fractal curve, the Dragon Curve. Download my new `dragon` program from the File Exchange and follow along.

Start with a long, narrow strip of paper, like the one in the following photograph. Fold it in half once, end over end. Then fold it a second time, being careful to fold in the same direction, right over left or left over right, as the first. Then a third fold, and a fourth. Unfold the strip so that all the folds form the right angles pictured. You will have created a Dragon Curve of degree four.

Image credit: http://www.cutoutfoldup.com/216-dragon-curve.php

How many times do you imagine you could fold the paper before it is too thick to fold again? Probably five, maybe six.

The degree one curve has one right angle; let's denote that by R. The degree two curve has two right angles, followed by a left; that's RRL. The degree three signature is RRLRRLL.

Can you read off the signature of degree four from the photo? Do you see the pattern? Probably not yet.

The number of segments between folds doubles with each increase in degree. Here's code to display the signatures of degrees one through six. That's all I can display in one column of this blog post.

The key operation in the code is `LR-fliplr(s)`. That flips the string `s` end to end and replaces all the R's by L's and the L's by R's.

s = 'R'; disp(s) LR = 'L'+'R'; for n = 2:6 s = [s 'R' LR-fliplr(s)]; disp(s) end

R RRL RRLRRLL RRLRRLLRRRLLRLL RRLRRLLRRRLLRLLRRRLRRLLLRRLLRLL RRLRRLLRRRLLRLLRRRLRRLLLRRLLRLLRRRLRRLLRRRLLRLLLRRLRRLLLRRLLRLL

The following figures are all from my new MATLAB program, `dragon`. The computations are based on a growing and shrinking vector `z` of points in the complex plane.

Start with

z = [0 1]

This isn't complex yet, but it soon will be. When we plot `z`, we get a single black line, representing an unfolded strip of paper.

Let's fold the paper once. This is done with the help of the complex number

w = (1+i)/2

In polar form, $w$ is $e^{i \pi/4}/\sqrt{2}$. So complex multiplication by $w$ rotates by $45^\circ$ and scales by $1/\sqrt{2}$.

In MATLAB, `w'` is the complex conjugate of `w`. This is used in the "fold" operation for any vector `z`.

zleft = w*z zright = 1 - w'*fliplr(z) z = [zleft zright]

Folding turns `[0 1]` into

[0 w w 1]

When we plot any vector, we break it in half, then plot the left half in black and the right half in blue.

The black line has been rotated $45^\circ$ to the left and a blue line, rotated $45^\circ$ to the right, appended.

Fold it again.

The black line becomes a black right angle and the blue line becomes a blue "left angle".

Two more folds produce a degree four Dragon Curve.

A total of eight folds gives degree eight.

And after a total of eighteen folds, we have our dragon.

We could keep on folding, but we now have $2^{18}$ short line segments and have reached the limits of screen resolution. More folds do not change any more pixels.

Let's do a different kind of rotation. Rotate the entire dragon by $90^\circ$, superimpose it on the original, and plot the four halves with four colors.

Do it twice more.

We could now do translations in all four compass directions and superimpose all the dragons, but I won't show that. It turns out that as both the degrees and the translates go to infinity, the curves, which are intricately intertwined, never intersect, and completely cover and tile the complex plane.

I found this at an on-line shopping center named Etsy, in the "Fractal and Geeky Jewelry" section. It's from an outfit named DragonNerd in Winnipeg. Count the number of bends to determine the degree.

The last time I blogged about jewelry was in my piece about the Pentium Affair.

Image credit: http://www.etsy.com/shop/DragonNerd

All the figures above are from my new program, `dragon`, available from MATLAB Central and also included in Cleve's Laboratory.

`dragon` provides four pushbuttons

These buttons

- Unwind the rotations
- Reduce the degree
- Increase the degree
- Superimpose rotations

The Web is full of stuff about the Dragon Curve. Just ask Google and Wikipedia.

I highly recommend two videos from Numberphile. In one, Don Knuth tells about a wrong turn he made in the Dragon Curve ceramic displayed in his home, https://www.youtube.com/watch?v=v678Em6qyzk.

In the other, the popular British author, Rob Eastaway, describes the appearance of the Curve in Michael Crichton's *Jurassic Park*, https://www.youtube.com/watch?v=wCyC-K_PnRY.

Get
the MATLAB code

Published with MATLAB® R2018a

Today is Easter Sunday. Why? How is the date of Easter determined?... read more >>

]]>Today is Easter Sunday. Why? How is the date of Easter determined?

This is a revision of a blog post from 2013.

Easter Day is one of the most important events in the Christian calendar. It is also one of the most mathematically elusive. In fact, regularization of the observance of Easter was one of the primary motivations for calendar reform centuries ago.

Easter is linked to the Jewish Passover. The informal rule is that Easter Day is the first Sunday after the first full moon after the vernal equinox. But the ecclesiastical full moon and equinox involved in this rule are not always the same as the corresponding astronomical events, which, after all, depend upon the location of the observer on the earth.

I am using a revision of the `easter` program in my blog post from 2013 and *Experiments with MATLAB*. Let's check this year.

easter_2018 = datestr(easter(2018))

easter_2018 = '01-Apr-2018'

My MATLAB® program is based on the algorithm presented in the first volume of the classic series by Donald Knuth, The Art of Computer Programming. Knuth has used it in several publications to illustrate different programming languages. The task has often been the topic of an exercise in computer programming courses.

Knuth says that the algorithm is due to the Neapolitan astronomer Aloysius Lilius and the German Jesuit mathematician Christopher Clavious in the late 16th century and that it is used by most Western churches to determine the date of Easter Sunday for any year after 1582.

The date varies between March 22 and April 25.

The earth's orbit around the sun and the moon's orbit around the earth are not in sync. It takes the earth about 365.2425 days to orbit the sun. This is known as a tropical year. The moon's orbit around the earth is complicated, but an average orbit takes about 29.53 days. This is known as a synodic month. The fraction

```
year = 365.2425;
month = 29.53;
format rat
ratio = year/month
```

ratio = 6444/521

is not the ratio of small integers. However, in the 5th century BC, an astronomer from Athens named Meton observed that the ratio is very close to 235/19.

```
format short
ratio
meton = 235/19
```

ratio = 12.3685 meton = 12.3684

In other words, 19 tropical years is close to 235 synodic months. This Metonic cycle was the basis for the Greek calendar and is the key to the algorithm for determining Easter.

Here is the revised MATLAB program. Try other years.

Can you spot the change made to my old program in the blog post from 2013?

```
type easter
```

function dn = easter(y) % EASTER Date of Easter. % EASTER(y) is the datenum of Easter in year y. % Example: % datestr(easter(2020)) % % Ref: Donald Knuth, The Art of Computer Programming, % Fundamental Algorithms, pp. 155-156. % Copyright 2014-18 Cleve Moler % Copyright 2014-18 The MathWorks, Inc. % Golden number in 19-year Metonic cycle. g = mod(y,19) + 1; % Century number. c = floor(y/100) + 1; % Corrections for leap years and moon's orbit. x = floor(3*c/4) - 12; z = floor((8*c+5)/25) - 5; % Epact. e = mod(11*g+20+z-x,30); if (e==25 && g>11 || e==24), e = e + 1; end % Full moon. n = 44 - e; if n < 21, n = n + 30; end % Find a Sunday. d = floor(5*y/4) - x - 10; % Easter is a Sunday in March or April. d = n + 7 - mod(d+n,7); dn = datenum(y,3,32);

[1] Donald E. Knuth, The Art of Computer Programming, Volume 1: Fundamental Algorithms (3rd edition), pp. 159-160, Addison-Wesley, 1997, ISBN 0-201-89683-4, PDF available.

[2] Wikipedia, Primary article on Easter. <http://en.wikipedia.org/wiki/Easter>

[3] Wikipedia, Computus, details on calculation of Easter. <http://en.wikipedia.org/wiki/Computus>

[4] Wikipedia, Metonic cycle. <http://en.wikipedia.org/wiki/Metonic_cycle>

Get
the MATLAB code

Published with MATLAB® R2018a

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.... 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 fifth such installment and the first of a multipart post about modern MATLAB.

Despite its name, MATLAB is not just a Matrix Laboratory any more. It has evolved over the last 35 years into a rich technical computing environment.

PC-MATLAB had only about a dozen graphics commands, including `plot` for two-dimensional work and `mesh` for three dimensions. The results were rendered on low resolution displays with limited memory and few, if any, colors. But as the equipment evolved, so did MATLAB graphics capabilities.

In 1986 we released a version of MATLAB named PRO-MATLAB, for UNIX workstations including the Sun-3. The displays and the window systems on these workstations supported more powerful graphics.

In 1992 we released MATLAB 4 for both PCs and workstations. This release included significant new graphics features, including color. I remember struggling with conversion of the RGB color represention that was the basis for our software to the CMYK representation that was expected by the printers of our paper documention. The problem is under specified. We have three values, R, G and B. They are weighted average of four values, C, M, Y and K. What are those four values? I was immensely pleased when I saw the first successful page proofs for the color insert.

This is not the place to describe all the graphics features we have in MATLAB today. Here's a sample, four different views of one of our venerable test surfaces.

[x,y,z] = peaks;

$$ z = 3(1-x) e^{-x^2-(y+1)^2} -10(\textstyle\frac{x}{5}-x^3-y^5) e^{-x^2-y^2} -\textstyle\frac{1}{3} e^{-(x+1)^2-y^2} $$

% First, an unadorned surface plot. p1 = subplot(2,2,1); surf(x,y,z) axis tight % Next, a 3D bar chart, with some adjustment of the default annotation. p2 = subplot(2,2,2); k = 1:4:49; b = bar3(flipud(z(k,k))); fixup(p2,b) % Then, a 2D plot of the columns of the array. p3 = subplot(2,2,3); plot(z) % Finally, a filled contour plot with 20 levels and a modest color map. p4 = subplot(2,2,4); contourf(z,20); colormap(p4,pink) axis off

For many years MATLAB had only one numeric data type, IEEE standard 754 double precision floating point, stored in the 64-bit format.

```
clear
format long
phi = (1 + sqrt(5))/2
```

phi = 1.618033988749895

Over several releases, we added support for single precision. Requiring only 32 bits of storage, single precision cuts memory requirements for large arrays in half. It may or may not be faster.

MATLAB does not have declarations, so single precision variables are obtained from double precision ones by an executable conversion function.

p = single(phi)

p = single 1.6180340

In 2004 MATLAB 7 introduced three unsigned integer data types, `uint32`, `uint16` and `uint8`, three signed integer data types, `int32`, `int16`, and `int8`, and one logical data type, `logical`.

q = uint16(1000*phi) r = int8(-10*phi) s = logical(phi)

q = uint16 1618 r = int8 -16 s = logical 1

Let's see how much storage is requied for these quantities

whos

Name Size Bytes Class Attributes p 1x1 4 single phi 1x1 8 double q 1x1 2 uint16 r 1x1 1 int8 s 1x1 1 logical

In 1992 MATLAB 4 introduced sparse matrices. Only the nonzero elements are stored, along with row indices and pointers to the starts of columns. The only change to the outward appearance of MATLAB is a pair of functions, `sparse` and `full`. Nearly all the operations apply equally to full or sparse matrices. The sparse storage scheme represents a matrix in space proportional to the number of nonzero entries, and most of the operations compute sparse results in time proportional to the number of arithmetic operations on nonzeros.

For example, consider the classic finite difference approximation to the Laplacian differential operator. The function `numgrid` numbers the points in a two-dimensional grid, in this case n-by-n points in the interior of a square.

```
clear
n = 100;
S = numgrid('S',n+2);
```

The function `delsq` creates the five-point discrete Laplacian, stored as a sparse N-by-N matrix, where N = n^2. With five or fewer nonzeros per row, the total number of nonzeros is a little less than 5*N.

A = delsq(S); nz = nnz(A)

nz = 49600

For the sake of comparison, let's create the full version of the matrix and check the amount of storage required. Storage required for `A` is proportional to N, while for `F` it is proportion to N^2.

F = full(A); whos

Name Size Bytes Class Attributes A 10000x10000 873608 double sparse F 10000x10000 800000000 double S 102x102 83232 double n 1x1 8 double nz 1x1 8 double

Let's time the solution to a boundary value problem. For the sparse matrix the time is O(N^2). For n = 100 it's instanteous.

b = ones(n^2,1); tic u = A\b; toc

Elapsed time is 0.018341 seconds.

The full matrix time is O(N^3). It would require several seconds to compute the same solution.

tic u = F\b; toc

Elapsed time is 7.810682 seconds.

Cell arrays were introduced with MATLAB 5 in 1996. A cell array is an indexed, possibly inhomogeneous collection of MATLAB objects, including other cell arrays. Cell arrays are created by curly braces, {}.

```
c = {magic(3); uint8(1:10); 'hello world'}
```

c = 3×1 cell array {3×3 double } {1×10 uint8 } {'hello world'}

Cell arrays can be indexed by both curly braces and smooth parentheses. With braces, c{k} is the contents of the k-th cell. With parentheses, c(k) is another cell array containing the specified cells.

M = c{1} c2 = c(1)

M = 8 1 6 3 5 7 4 9 2 c2 = 1×1 cell array {3×3 double}

Think of a cell array as a collection of mailboxes. `box(k)` is the `k`-th mail box. `box{k}` is the mail in the `k`-th box.

For many years, text was a second class citizen in MATLAB.

With concerns about cross platform portability, Historic MATLAB had its own internal character set. Text delineated by single quotes was converted a character at a time into floating point numbers in the range 0:51. Lower case was converted to upper. This process was reversed by the DISP function.

Here is some output from Historic MATLAB.

<> H = 'hello world'

H = 17 14 21 21 24 36 32 24 27 21 13

<> disp(H)

HELLO WORLD

MathWorks versions of MATLAB have relied on the ASCII character set, including upper and lower case. Character vectors are still delineated by single quotes.

```
h = 'hello world'
```

h = 'hello world'

disp(h)

hello world

d = uint8(h)

d = 1×11 uint8 row vector 104 101 108 108 111 32 119 111 114 108 100

Short character strings are often used as optional parameters to functions.

[U,S,V] = svd(A,'econ') % Economy size, U is the same shape as A.

plot(x,y,'o-') % Plot lines with circles at the data points.

Multiple lines of text, or many words in an array, must be padded with blanks so that the character array is rectangular. The `char` function provides this service. For example, here is a 3-by-7 array.

cast = char('Alice','Bob','Charlie')

cast = 3×7 char array 'Alice ' 'Bob ' 'Charlie'

Or, you could use a cell array.

cast = {'Alice', 'Bob', 'Charlie'}'

cast = 3×1 cell array {'Alice' } {'Bob' } {'Charlie'}

In 2016, we began the process of providing more comprehensive support for text by introducing the double quote character and the `string` data structure.

cast = ["Alice", "Bob", "Charlie"]'

cast = 3×1 string array "Alice" "Bob" "Charlie"

There are some very convenient functions for the new string data type.

```
proverb = "A rolling stone gathers momentum"
words = split(proverb)
```

proverb = "A rolling stone gathers momentum" words = 5×1 string array "A" "rolling" "stone" "gathers" "momentum"

Addition of strings is concatenation.

merge = cast(1); plus = " + "; for k = 2:length(cast) merge = merge + plus + cast(k); end merge

merge = "Alice + Bob + Charlie"

The `regexp` function provides the *regular expression* pattern matching seen in Unix and many other programming languages.

The distinction between functions and commands was never clear in the early days of MATLAB. This was eventually resolved by *the command/function duality*: A command statement of the form

cmd arg1 arg2

is the same function statement with `char` arguments

```
cmd('arg1,'arg2')
```

For example

```
diary notes.txt
```

is the same as

```
f = 'notes.txt';
diary(f)
```

MATLAB has functions that take other functions as arguments. These are known as "function functions". For many years function arguments were specified by a string. For example, suppose we want to numerically solve the ordinary differential equation for the Van de Pol oscillator. Create a file `vanderpol.m` that evaluates the differential equation.

```
type vanderpol
```

function dydt = vanderpol(t,y) dydt = [y(2); 5*(1-y(1)^2)*y(2)-y(1)]; end

Then pass the function name as a string to the ODE solver `ode45`.

```
tspan = [0 150];
y0 = [1 0]';
[t,y] = ode45('vanderpol',tspan,y0);
```

It turns out that over half the time required by this computation is spent in repeatedly decoding the string argument. To improve performance we introduced the *function handle*, which is the function name preceeded by the "at sign", `@vanderpol`.

The "at sign" is also used to create a function handle that defines an *anonymous function*, which is the MATLAB instantiation of Church's lambda calculus.

@(x) sin(x)./x

Ironically anonymous functions are often assigned to variables, thereby negating the anonymity.

sinc = @(x) sin(x)./x; vdp = @(t,y) [y(2); 5*(1-y(1)^2)*y(2)-y(1)];

Let's compare the times required in our Van der pol integration with strings, function handles and anonymous functions.

```
tic, [t,y] = ode45('vanderpol',tspan,y0); toc
tic, [t,y] = ode45(@vanderpol,tspan,y0); toc
tic, [t,y] = ode45(vdp,tspan,y0); toc
```

Elapsed time is 0.173259 seconds. Elapsed time is 0.041710 seconds. Elapsed time is 0.032854 seconds.

We see that the times required for the later two are comparable, and are significantly faster than the first.

Get
the MATLAB code

Published with MATLAB® R2018a

Happy Pi Day, 3/14.... read more >>

]]>Happy Pi Day, 3/14.

Here are 10,000 digits of π. Notice the six consecutive 9's in the sixth row, digits 763 through 768.

n = 100; % Get pi from Symbolic Toolbox. p = vpa(pi,n^2); % Convert digits to characters, then to individual integers. p = char(p); p(2) = []; % Decimal point p = p - '0'; p(end+1:n^2) = 0; % Square array P = reshape(p,n,n)'; % Extra row and column to satisfy pcolor. P(n+1,n+1) = NaN; % Plot with pcolor. pcolor(P) axis square axis off % Title title('\pi','fontsize',24) % Make room for colorbar pos = get(gca,'pos') - [.05 0 .0 0]; set(gca,'ydir','rev','pos',pos) % Custom colorbar. a = axes('pos',[.80 .60 .02 .30]); y = (0:9)'; y(11,2) = NaN; pcolor(y) set(a,'xtick',[],'ydir','rev','ytick',(1.5:1:10.5)', ... 'yaxislocation','right', ... 'yticklabel',{' 0',' 1',' 2',' 3',' 4',' 5',' 6',' 7',' 8',' 9'})

Get
the MATLAB code

Published with MATLAB® R2018a

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.... 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 fourth such installment. MATLAB takes the first steps from a primitive matrix calculator to a full-fledged technical computing environment.

I spent the 1979-80 academic year at Stanford, as a Visiting Professor of Computer Science. In the fall I taught CS237a, the graduate course in Numerical Analysis. I introduced the class to the matrix calculator that I called Historic MATLAB in my previous post about MATLAB history.

There were maybe fifteen to twenty students in the class. About half of them were Math and CS students. As I remember, they were not impressed with MATLAB. It was not a particularly sophisticated programming language. It was not numerical analysis research.

But the other half of the students were from engineering at Stanford. They loved MATLAB. They were studying subjects like control theory and signal processing that, back then, I knew nothing about. Matrices were a central part of the mathematics in these subjects. The students had been doing small matrix problems by hand and larger problems by writing Fortran programs. MATLAB was immediately useful.

Jack Little studied electrical engineering at MIT in the late 1970's and then came west to Stanford for grad school. He didn't take my CS237 course, but a friend of his did. The friend showed him MATLAB and Little immediately adopted it for his own work in control systems and signal processing.

In 1983 Little suggested the creation of a commercial product based on MATLAB. I said I thought that was a good idea, but I didn't join him initially. The IBM PC had been introduced only two years earlier and was barely powerful enough to run something like MATLAB, but Little anticipated its evolution. He left his job, bought a Compaq PC clone at Sears, moved into the hills behind Stanford, and, with my encouragement, spent a year and a half creating a new and extended version of MATLAB written in C. A friend, Steve Bangert, joined the project and worked on the new MATLAB in his spare time.

The three of us -- Little, Bangert, and myself -- formed MathWorks in California in 1984. PC-MATLAB made its debut in December 1984 at the IEEE Conference on Decision and Control in Las Vegas.

Little and Bangert made many important modifications and improvements to Historic MATLAB when they created the new and extended version. The most significant were functions, toolboxes and graphics

The MATLAB function naming mechanism uses the underlying computer file system. A *script* or *function* is a body of MATLAB code stored in a file with the extension `.m`. If the file begins with the keyword `function`, then it is a function, otherwise it is a script. The name of the file, minus the `.m`, is the name of the script or function. For example, a file named `hello.m` containing the single line

disp('Hello World')

is a MATLAB script. Typing the file name without the `.m` extension at the command line produces the traditional greeting.

hello

Hello World

Functions usually have input and output arguments. Here is a simple example that uses a vector outer product to generate the multiplication table you memorized in elementary school

```
type mults
```

function T = mults(n) j = 1:n; T = j'*j; end

This statement produces a 10-by-10 multiplication table.

T = mults(10)

T = 1 2 3 4 5 6 7 8 9 10 2 4 6 8 10 12 14 16 18 20 3 6 9 12 15 18 21 24 27 30 4 8 12 16 20 24 28 32 36 40 5 10 15 20 25 30 35 40 45 50 6 12 18 24 30 36 42 48 54 60 7 14 21 28 35 42 49 56 63 70 8 16 24 32 40 48 56 64 72 80 9 18 27 36 45 54 63 72 81 90 10 20 30 40 50 60 70 80 90 100

Input arguments to MATLAB functions are "passed by value" or, more precisely, passed by "copy on write". It is not necessary to make copy of an input argument if the function does not change it. For example, the function

function C = commuator(A,B) C = A*B - B*A; end

does not change `A` or `B`, so it does not copy either.

Toolboxes are collections of functions written in MATLAB that can be read, examined, and possibly modified, by users. The `matlab` toolbox, which is an essential part of every installation, has several dozen functions that provide capabilities beyond the built-in functions in the core.

Other toolboxes are available from MathWorks and many other sources. The first two specialized toolboxes were control and signal.

Technical computer graphics have been an essential part of MATLAB since the first MathWorks version. Let's look at two examples from that first version. Notice that neither one has anything to do with numerical linear algebra, but they both rely on vector notation and array operations. One demo was called simply `wow`.

```
type wow
```

% Lets define a peculiar vector, and then plot % a growing cosine wave versus a growing sine: t = 0:(.99*pi/2):500; x = t.*cos(t); y = t.*sin(t); plot(x,y,'k') axis square

wow

The second example investigates the Gibbs phenomena by plotting the successive partial sums of the Fourier series for a square wave.

```
type square_wave
```

% For a finale, we can go from the fundamental to the 19th harmonic % and create vectors of successively more harmonics, saving all % intermediate steps as rows in a matrix. t = 0:pi/128:pi; y = zeros(10,max(size(t))); x = sin(t); y(1,:) = x; i = 1; for k = 3:2:19 i = i + 1; x = x + sin(k*t)/k; y(i,:) = x; end % We can plot these successively on an overplot, showing the transition % to a square wave. Note that Gibb's effect says that you will never % really get there. % Or we can plot this as a 3-d mesh surface. mesh(y) colormap([0 0 0]) % Black and white. axis tight axis off view([-27 39])

square_wave

I cannot find a disc for the very first version of PC-MATLAB. The oldest I can find is version 1.3. Here are all the functions, key words, and operators in that version.

**Information and assistance:**

help - help facility (help intro, help index) demo - run demonstrations browse - browse utility who - lists current variables in memory what - lists .M files on disk dir - directory list of files on disk casesen- turn off case sensitivity format - set output format pack - memory garbage collection and compaction edit - invoke editor

**Attributes and Array Manipulation:**

max - maximum value in vector min - minimum value in vector sum - sum of elements in vector prod - product of elements in vector cumsum - cumulative sum of elements in vector cumprod - cumulative product of elements in vector mean - mean value in a vector median - median value in a vector std - standard deviation in a vector size - row and column dimensions of an array length - length of a vector : - array subscripting, selection, and removing sort - sorting find - find array indices of logical values hist - histograms

**Array Building Functions:**

diag - creates diagonal arrays, removes diagonal vectors eye - creates identity matrices ones - generates an array of all ones zeros - generates an array of all zeros rand - random numbers and arrays magic - generates a magic square tril - lower triangular part triu - upper triangular part hilb - generates Hilbert matrix invhilb - generates inverse Hilbert matrix

**Matrix Properties:**

cond - condition number in 2-norm det - determinant norm - 1-norm, 2-norm, F-norm, infinity norm rank - rank rcond - condition estimate

**Matrix Functions:**

inv - inverse expm - matrix exponential logm - matrix logarithm sqrtm - matrix square root funm - arbitrary matrix functions pinv - pseudoinverse poly - characteristic polynomial kron - Kronecker tensor product

**Matrix Decompositions and Factorizations:**

chol - Cholesky factorization eig - eigenvalues and eigenvectors hess - Hessenberg form lu - factors from Gaussian elimination orth - orthogonalization qr - orthogonal-triangular decomposition schur - Schur decomposition svd - singular value decomposition

**Polynomial manipulations:**

poly - characteristic polynomial roots - find polynomial roots polyval - evaluate polynomial conv - multiplication

**Signal Processing:**

filter - digital filter fft - fast Fourier transform ifft - inverse fast Fourier transform conv - convolution length - length of a sequence

**Permanent Variables:**

ans - answer when expression is not assigned eps - machine epsilon, floating point relative accuracy pi - 3.14159265358979323846264338327950288419716939937511... inf - infinity nan - Not-a-Number nargin - number of input arguments in a function nargout - number of output arguments in a function

**Control Flow:**

if s conditionally execute statements elseif s else end for i=v,.. end repeat statements a specific number of times while s,..,end do while break break out of a for or while loop return return from function

**Graphics and Plotting:**

shg - show graphics screen cla - clear alpha screen clg - clear graphics screen plot - linear X-Y plot loglog - loglog X-Y plot semilogx - semi-log X-Y plot semilogy - semi-log X-Y plot polar - polar plot mesh - 3-dimensional mesh surface title - plot title xlabel - x-axis label ylabel - y-axis label grid - draw grid lines axis - manual axis scaling print - hardcopy disp - compact matrix display

**User-defined Functions, Procedures, and Programming:**

type - list a Function or file what - lists .M files on disk input - get a number from the user pause - pause until <cr> entered return - return from a user defined Function eval - interpret text in a variable setstr - set flag indicating matrix is a string echo - enables command echoing nargin - number of input arguments in a function nargout - number of output arguments in a function exist - checks if a variable or function exists

**Disk Files:**

diary - diary of the session in a disk file load - loads variables from a file save - saves variables on a file type - list a file dir - directory list of files on disk delete - delete files chdir - change working directory

**Elementary Math Functions:**

abs - absolute value or complex magnitude conj - complex conjugate fix - round towards zero imag - imaginary part real - real part round - round to nearest integer sign - signum function sqrt - square root rem - remainder sin, cos, tan asin, acos, atan atan2 - four quadrant arctangent exp - exponential base e log - natural logarithm log10 - log base 10 bessel - Bessel functions rat - rational number approximation

**Special Symbols:**

= - assignment statement [ - used to form vectors and matrices ] - see [ ( - arithmetic expression precedence ) - see ( . - Decimal point .. - continue a statement to the next line , - separates matrix subscripts and function arguments ; - Ends rows. Suppresses the display of a variable % - comment statement : - subscripting, vector generation, data selection ' - transpose, string delimiter

**Element-by-element arithmetic operators:**

+ - addition - - subtraction .* - multiplication ./ - division .\ - left division .^ - raise to a power .' - transpose

**Matrix operators:**

* - multiplication / - division \ - left division ^ - raise to a power ' - matrix conjugate transpose

**Relational operators:**

< - less than <= - less than or equals > - greater than >= - greater than or equals == - equals ~= - not equals

**Logical operators:**

+ - OR .* - AND ~ - NOT (complement)

Get
the MATLAB code

Published with MATLAB® R2018a

```
syms x
n = 7;
x7 = expand((x+1)^n)
```

x7 = x^7 + 7*x^6 + 21*x^5 + 35*x^4 + 35*x^3 + 21*x^2 + 7*x + 1Formally, the binomial coefficients are given by $${{n} \choose {k}} = \frac {n!} {k! (n-k)!}$$ But premature floating point overflow of the factorials makes this an unsatisfactory basis for computation. A better way employs the recursion $$ {{n} \choose {k}} = {{n-1} \choose {k}} + {{n-1} \choose {k-1}}$$ This is used by the MATLAB function

P = pascal(7)

P = 1 1 1 1 1 1 1 1 2 3 4 5 6 7 1 3 6 10 15 21 28 1 4 10 20 35 56 84 1 5 15 35 70 126 210 1 6 21 56 126 252 462 1 7 28 84 210 462 924The other is lower triangular, with the binomial coefficients in the rows. (We will see why the even numbered columns have minus signs in a moment.)

L = pascal(7,1)

L = 1 0 0 0 0 0 0 1 -1 0 0 0 0 0 1 -2 1 0 0 0 0 1 -3 3 -1 0 0 0 1 -4 6 -4 1 0 0 1 -5 10 -10 5 -1 0 1 -6 15 -20 15 -6 1The individual elements are

P(i,j) = P(j,i) = nchoosek(i+j-2,j-1)And (temporarily ignoring the minus signs) for

L(i,j) = nchoosek(i-1,j-1)The first fun fact is that

L = chol(P)'

L = 1 0 0 0 0 0 0 1 1 0 0 0 0 0 1 2 1 0 0 0 0 1 3 3 1 0 0 0 1 4 6 4 1 0 0 1 5 10 10 5 1 0 1 6 15 20 15 6 1So we can reconstruct

P = L*L'

P = 1 1 1 1 1 1 1 1 2 3 4 5 6 7 1 3 6 10 15 21 28 1 4 10 20 35 56 84 1 5 15 35 70 126 210 1 6 21 56 126 252 462 1 7 28 84 210 462 924

triprint(L)

1 1 1 1 2 1 1 3 3 1 1 4 6 4 1 1 5 10 10 5 1 1 6 15 20 15 6 1

L = pascal(n,1) L_squared = L^2

L = 1 0 0 0 0 0 0 1 -1 0 0 0 0 0 1 -2 1 0 0 0 0 1 -3 3 -1 0 0 0 1 -4 6 -4 1 0 0 1 -5 10 -10 5 -1 0 1 -6 15 -20 15 -6 1 L_squared = 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1Here is an exercise for you. What is

X = rot90(L,-1) X_cubed = X^3

X = 1 1 1 1 1 1 1 -6 -5 -4 -3 -2 -1 0 15 10 6 3 1 0 0 -20 -10 -4 -1 0 0 0 15 5 1 0 0 0 0 -6 -1 0 0 0 0 0 1 0 0 0 0 0 0 X_cubed = 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1

```
odd = @(x) mod(x,2)==1;
n = 56;
L = abs(pascal(n,1));
spy(odd(L))
title('odd(L)')
```

n = 12; A = fliplr(abs(pascal(n,1))) for k = 1:n F(k) = sum(diag(A,n-k)); end F

A = 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 1 2 1 0 0 0 0 0 0 0 0 1 3 3 1 0 0 0 0 0 0 0 1 4 6 4 1 0 0 0 0 0 0 1 5 10 10 5 1 0 0 0 0 0 1 6 15 20 15 6 1 0 0 0 0 1 7 21 35 35 21 7 1 0 0 0 1 8 28 56 70 56 28 8 1 0 0 1 9 36 84 126 126 84 36 9 1 0 1 10 45 120 210 252 210 120 45 10 1 1 11 55 165 330 462 462 330 165 55 11 1 F = 1 1 2 3 5 8 13 21 34 55 89 144

L = pascal(12,1); t = L(3:end,3)'

t = 1 3 6 10 15 21 28 36 45 55Here's an unusual series relating the triangle numbers to $\pi$. The signs go + + - - + + - - .

pi - 2 = 1 + 1/3 - 1/6 - 1/10 + 1/15 + 1/21 - 1/28 - 1/36 + 1/45 + 1/55 - ...

```
type pi_pascal
```

function pie = pi_pascal(n) tk = 1; s = 1; for k = 2:n tk = tk + k; if mod(k+1,4) > 1 s = s + 1/tk; else s = s - 1/tk; end end pie = 2 + s;Ten million terms gives $\pi$ to 14 decimal places.

```
format long
pie = pi_pascal(10000000)
err = pi - pie
```

pie = 3.141592653589817 err = -2.398081733190338e-14

D = diag(1:7,-1)

D = 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 2 0 0 0 0 0 0 0 0 3 0 0 0 0 0 0 0 0 4 0 0 0 0 0 0 0 0 5 0 0 0 0 0 0 0 0 6 0 0 0 0 0 0 0 0 7 0is

expm_D = round(expm(D))

expm_D = 1 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 1 2 1 0 0 0 0 0 1 3 3 1 0 0 0 0 1 4 6 4 1 0 0 0 1 5 10 10 5 1 0 0 1 6 15 20 15 6 1 0 1 7 21 35 35 21 7 1

Edward Scheinerman, The Mathematics Lover's Companion, Yale University Press, 2017 https://yalebooks.yale.edu/book/9780300223002/mathematics-lovers-companion... read more >>

]]>ezplot('(x^2+y^2-1)^3 = x^2*y^3', ... [-1.5,1.5,-1.2,1.5]) colormap([.5 .2 .2])

Edward Scheinerman, *The Mathematics Lover's Companion*, Yale University Press, 2017 https://yalebooks.yale.edu/book/9780300223002/mathematics-lovers-companion

Get
the MATLAB code

Published with MATLAB® R2018a