Download Golden Ratio and Fibonacci Numbers: Mathematical Concepts in MATLAB and more Exams Matlab skills in PDF only on Docsity!
Chapter 1
Introduction to MATLAB
This book is an introduction to two subjects: Matlab and numerical computing. This first chapter introduces Matlab by presenting several programs that inves- tigate elementary, but interesting, mathematical problems. If you already have some experience programming in another language, we hope that you can see how Matlab works by simply studying these programs. If you want a more comprehensive introduction, there are many resources available. You can select the Help tab in the toolstrip atop the Matlab com- mand window, then select Documentation, MATLAB and Getting Started. A MathWorks Web site, MATLAB Tutorials and Learning Resources [11], offers a number of introductory videos and a PDF manual entitled Getting Started with MATLAB. An introduction to MATLAB through a collection of mathematical and com- putational projects is provided by Moler’s free online Experiments with MATLAB [6]. A list of over 1500 Matlab-based books by other authors and publishers, in several languages, is available at [12]. Three introductions to Matlab are of par- ticular interest here: a relatively short primer by Sigmon and Davis [9], a medium- sized, mathematically oriented text by Higham and Higham [3], and a large, com- prehensive manual by Hanselman and Littlefield [2]. You should have a copy of Matlab close at hand so you can run our sample programs as you read about them. All of the programs used in this book have been collected in a directory (or folder) named
NCM
(The directory name is the initials of the book title.) You can either start Matlab in this directory or use
pathtool
to add the directory to the Matlab path.
September 18, 2013
2 Chapter 1. Introduction to MATLAB
1.1 The Golden Ratio
What is the world’s most interesting number? Perhaps you like π, or e, or 17. Some people might vote for ϕ, the golden ratio, computed here by our first Matlab statement.
phi = (1 + sqrt(5))/
This produces
phi =
Let’s see more digits.
format long phi
phi =
This didn’t recompute ϕ, it just displayed 16 significant digits instead of 5. The golden ratio shows up in many places in mathematics; we’ll see several in this book. The golden ratio gets its name from the golden rectangle, shown in Figure 1.1. The golden rectangle has the property that removing a square leaves a smaller rectangle with the same shape.
φ
φ − 1
1
1
Figure 1.1. The golden rectangle.
Equating the aspect ratios of the rectangles gives a defining equation for ϕ:
1 ϕ
ϕ − 1 1
This equation says that you can compute the reciprocal of ϕ by simply subtracting one. How many numbers have that property? Multiplying the aspect ratio equation by ϕ produces the polynomial equation
ϕ^2 − ϕ − 1 = 0.
4 Chapter 1. Introduction to MATLAB
The variable r is a vector with two components, the symbolic forms of the two solutions. You can pick off the first component with
phi = r(1)
which produces
phi = 5^(1/2)/2 + 1/ This expression can be converted to a numerical value in two different ways. It can be evaluated to any number of digits using variable-precision arithmetic with the vpa function.
vpa(phi,50)
produces 50 digits.
It can also be converted to double-precision floating point, which is the principal way that Matlab represents numbers, with the double function.
phi = double(phi)
produces
phi =
The aspect ratio equation is simple enough to have closed-form symbolic so- lutions. More complicated equations have to be solved approximately. In Matlab an anonymous function is a convenient way to define an object that can be used as an argument to other functions. The statement
f = @(x) 1./x-(x-1)
defines f (x) = 1/x − (x − 1) and produces
f = @(x) 1./x-(x-1)
The graph of f (x) over the interval 0 ≤ x ≤ 4 shown in Figure 1.2 is obtained with
1.1. The Golden Ratio 5
0 0.5 1 1.5 2 2.5 3 3.5 4
−
−
−
0
1
2
3
4
5
6
7
x
1/x − (x−1)
Figure 1.2. f (ϕ) = 0.
ezplot(f,0,4)
The name ezplot stands for “easy plot,” although some of the English-speaking world would pronounce it “e-zed plot.” Even though f (x) becomes infinite as x → 0, ezplot automatically picks a reasonable vertical scale. The statement
phi = fzero(f,1)
looks for a zero of f (x) near x = 1. It produces an approximation to ϕ that is accurate to almost full precision. The result can be inserted in Figure 1.2 with
hold on plot(phi,0,’o’)
The following Matlab program produces the picture of the golden rectangle shown in Figure 1.1. The program is contained in an M-file named goldrect.m, so issuing the command
goldrect runs the script and creates the picture.
% GOLDRECT Plot the golden rectangle
phi = (1+sqrt(5))/2; x = [0 phi phi 0 0]; y = [0 0 1 1 0]; u = [1 1]; v = [0 1];
1.1. The Golden Ratio 7
end p = sprintf(’%d/%d’,p,q)
format long p = eval(p)
format short err = (1+sqrt(5))/2 - p
The statement
goldfract(6)
produces
p = 1+1/(1+1/(1+1/(1+1/(1+1/(1+1/(1))))))
p = 21/
p =
err =
The three p’s are all different representations of the same approximation to ϕ. The first p is the continued fraction truncated to six terms. There are six right parentheses. This p is a string generated by starting with a single ‘ 1 ’ (that’s goldfract(0)) and repeatedly inserting the string ‘1+1/(’ in front and the string ‘)’ in back. No matter how long this string becomes, it is a valid Matlab expression. The second p is an “ordinary” fraction with a single integer numerator and denominator obtained by collapsing the first p. The basis for the reformulation is
p q
p + q p
So the iteration starts with 1 1 and repeatedly replaces the fraction
p q
with p + q p
The statement
8 Chapter 1. Introduction to MATLAB
p = sprintf(’%d/%d’,p,q)
prints the final fraction by formatting p and q as decimal integers and placing a ‘/’ between them. The third p is the same number as the first two p’s, but is represented as a conventional decimal expansion, obtained by having the Matlab eval function actually do the division expressed in the second p. The final quantity err is the difference between p and ϕ. With only 6 terms, the approximation is accurate to less than 3 digits. How many terms does it take to get 10 digits of accuracy? As the number of terms n increases, the truncated continued fraction generated by goldfract(n) theoretically approaches ϕ. But limitations on the size of the integers in the numerator and denominator, as well as roundoff error in the actual floating-point division, eventually intervene. Exercise 1.3 asks you to investigate the limiting accuracy of goldfract(n).
1.2 Fibonacci Numbers
Leonardo Pisano Fibonacci was born around 1170 and died around 1250 in Pisa in what is now Italy. He traveled extensively in Europe and Northern Africa. He wrote several mathematical texts that, among other things, introduced Europe to the Hindu-Arabic notation for numbers. Even though his books had to be tran- scribed by hand, they were widely circulated. In his best known book, Liber Abaci, published in 1202, he posed the following problem:
A man put a pair of rabbits in a place surrounded on all sides by a wall. How many pairs of rabbits can be produced from that pair in a year if it is supposed that every month each pair begets a new pair which from the second month on becomes productive? Today the solution to this problem is known as the Fibonacci sequence, or Fibonacci numbers. There is a small mathematical industry based on Fibonacci numbers. A search of the Internet for “Fibonacci” will find dozens of Web sites and hundreds of pages of material. There is even a Fibonacci Association that publishes a scholarly journal, the Fibonacci Quarterly. If Fibonacci had not specified a month for the newborn pair to mature, he would not have a sequence named after him. The number of pairs would simply double each month. After n months there would be 2n^ pairs of rabbits. That’s a lot of rabbits, but not distinctive mathematics. Let fn denote the number of pairs of rabbits after n months. The key fact is that the number of rabbits at the end of a month is the number at the beginning of the month plus the number of births produced by the mature pairs:
fn = fn− 1 + fn− 2.
The initial conditions are that in the first month there is one pair of rabbits and in the second there are two pairs:
f 1 = 1, f 2 = 2.
10 Chapter 1. Introduction to MATLAB
FIBONACCI Fibonacci sequence f = FIBONACCI(n) generates the first n Fibonacci numbers.
The name of the function is in uppercase because historically Matlab was case insensitive and ran on terminals with only a single font. The use of capital letters may be confusing to some first-time Matlab users, but the convention persists. It is important to repeat the input and output arguments in these comments because the first line is not displayed when you ask for help on the function. The next line
f = zeros(n,1);
creates an n-by- 1 matrix containing all zeros and assigns it to f. In Matlab, a matrix with only one column is a column vector and a matrix with only one row is a row vector. The next two lines,
f(1) = 1; f(2) = 2;
provide the initial conditions. The last three lines are the for statement that does all the work.
for k = 3:n f(k) = f(k-1) + f(k-2); end
We like to use three spaces to indent the body of for and if statements, but other people prefer two or four spaces, or a tab. You can also put the entire construction on one line if you provide a comma after the first clause. This particular function looks a lot like functions in other programming lan- guages. It produces a vector, but it does not use any of the Matlab vector or matrix operations. We will see some of these operations soon. Here is another Fibonacci function, fibnum.m. Its output is simply the nth Fibonacci number.
function f = fibnum(n) % FIBNUM Fibonacci number. % FIBNUM(n) generates the nth Fibonacci number. if n <= 1 f = 1; else f = fibnum(n-1) + fibnum(n-2); end
The statement
fibnum(12)
produces
1.2. Fibonacci Numbers 11
ans = 233 The fibnum function is recursive. In fact, the term recursive is used in both a mathematical and a computer science sense. The relationship fn = fn− 1 + fn− 2 is known as a recursion relation and a function that calls itself is a recursive function. A recursive program is elegant, but expensive. You can measure execution time with tic and toc. Try
tic, fibnum(24), toc
Do not try
tic, fibnum(50), toc Now compare the results produced by goldfract(6) and fibonacci(7). The first contains the fraction 21/13 while the second ends with 13 and 21. This is not just a coincidence. The continued fraction is collapsed by repeating the statement
p = p + q;
while the Fibonacci numbers are generated by
f(k) = f(k-1) + f(k-2);
In fact, if we let ϕn denote the golden ratio continued fraction truncated at n terms, then fn+ fn
= ϕn.
In the infinite limit, the ratio of successive Fibonacci numbers approaches the golden ratio:
lim n→∞
fn+ fn
= ϕ.
To see this, compute 40 Fibonacci numbers. n = 40; f = fibonacci(n);
Then compute their ratios.
f(2:n)./f(1:n-1)
This takes the vector containing f(2) through f(n) and divides it, element by element, by the vector containing f(1) through f(n-1). The output begins with
1.3. Fractal Fern 13
This is an amazing equation. The right-hand side involves powers and quo- tients of irrational numbers, but the result is a sequence of integers. You can check this with Matlab, displaying the results in scientific notation.
format long e n = (1:40)’; f = (phi.^(n+1) - (1-phi).^(n+1))/(2*phi-1)
The .^ operator is an element-by-element power operator. It is not necessary to use ./ for the final division because (2*phi-1) is a scalar quantity. The computed result starts with
f = 1.000000000000000e+ 2.000000000000000e+ 3.000000000000000e+ 5.000000000000001e+ 8.000000000000002e+ 1.300000000000000e+ 2.100000000000000e+ 3.400000000000001e+
and ends with
5.702887000000007e+ 9.227465000000011e+ 1.493035200000002e+ 2.415781700000003e+ 3.908816900000005e+ 6.324598600000007e+ 1.023341550000001e+ 1.655801410000002e+
Roundoff error prevents the results from being exact integers, but
f = round(f)
finishes the job.
1.3 Fractal Fern
The M-files fern.m and finitefern.m produce the “Fractal Fern” described by Michael Barnsley in Fractals Everywhere [1]. They generate and plot a potentially infinite sequence of random, but carefully choreographed, points in the plane. The command
fern
runs forever, producing an increasingly dense plot. The command
14 Chapter 1. Introduction to MATLAB
Figure 1.3. Fractal fern.
finitefern(n)
generates n points and a plot like Figure 1.3. The command
finitefern(n,’s’)
shows the generation of the points one at a time. The command
F = finitefern(n);
generates, but does not plot, n points and returns an array of zeros and ones for use with sparse matrix and image-processing functions. The NCM collection also includes fern.png, a 768-by-1024 color image with half a million points that you can view with a browser or a paint program. You can
16 Chapter 1. Introduction to MATLAB
p = [ .85 .92 .99 1.00]; A1 = [ .85 .04; -.04 .85]; b1 = [0; 1.6]; A2 = [ .20 -.26; .23 .22]; b2 = [0; 1.6]; A3 = [-.15 .28; .26 .24]; b3 = [0; .44]; A4 = [ 0 0 ; 0 .16];
cnt = 1; tic while ~get(stop,’value’) r = rand; if r < p(1) x = A1x + b1; elseif r < p(2) x = A2x + b2; elseif r < p(3) x = A3x + b3; else x = A4x; end set(h,’xdata’,x(1),’ydata’,x(2)); cnt = cnt + 1; drawnow end t = toc; s = sprintf(’%8.0f points in %6.3f seconds’,cnt,t); text(-1.5,-0.5,s,’fontweight’,’bold’); set(stop,’style’,’pushbutton’,’string’,’close’, ... ’callback’,’close(gcf)’)
Let’s examine this program a few statements at a time.
shg
stands for “show graph window.” It brings an existing graphics window forward, or creates a new one if necessary.
clf reset
resets most of the figure properties to their default values.
set(gcf,’color’,’white’,’menubar’,’none’, ... ’numbertitle’,’off’,’name’,’Fractal Fern’)
changes the background color of the figure window from the default gray to white and provides a customized title for the window.
x = [.5; .5];
provides the initial coordinates of the point.
1.3. Fractal Fern 17
h = plot(x(1),x(2),’.’);
plots a single dot in the plane and saves a handle, h, so that we can later modify the properties of the plot.
darkgreen = [0 2/3 0];
defines a color where the red and blue components are zero and the green component is two-thirds of its full intensity.
set(h,’markersize’,1,’color’,darkgreen,’erasemode’,’none’);
makes the dot referenced by h smaller, changes its color, and specifies that the image of the dot on the screen should not be erased when its coordinates are changed. A record of these old points is kept by the computer’s graphics hardware (until the figure is reset), but Matlab itself does not remember them.
axis([-3 3 0 10]) axis off
specifies that the plot should cover the region
− 3 ≤ x 1 ≤ 3 , 0 ≤ x 2 ≤ 10 ,
but that the axes should not be drawn.
stop = uicontrol(’style’,’toggle’,’string’,’stop’, ... ’background’,’white’);
creates a toggle user interface control, labeled ’stop’ and colored white, in the default position near the lower left corner of the figure. The handle for the control is saved in the variable stop.
drawnow
causes the initial figure, including the initial point, to actually be plotted on the computer screen. The statement
p = [ .85 .92 .99 1.00];
sets up a vector of probabilities. The statements
A1 = [ .85 .04; -.04 .85]; b1 = [0; 1.6]; A2 = [ .20 -.26; .23 .22]; b2 = [0; 1.6]; A3 = [-.15 .28; .26 .24]; b3 = [0; .44]; A4 = [ 0 0 ; 0 .16];
define the four affine transformations. The statement
cnt = 1;
initializes a counter that keeps track of the number of points plotted. The statement
1.4. Magic Squares 19
displays the elapsed time since tic was called and the final count of the number of points plotted. Finally,
set(stop,’style’,’pushbutton’,’string’,’close’, ... ’callback’,’close(gcf)’)
changes the control to a push button that closes the window.
1.4 Magic Squares
Matlab stands for Matrix Laboratory. Over the years, Matlab has evolved into a general-purpose technical computing environment, but operations involving vectors, matrices, and linear algebra continue to be its most distinguishing feature. Magic squares provide an interesting set of sample matrices. The command help magic tells us the following:
MAGIC(N) is an N-by-N matrix constructed from the integers 1 through N^2 with equal row, column, and diagonal sums. Produces valid magic squares for all N > 0 except N = 2.
Magic squares were known in China over 2,000 years before the birth of Christ. The 3-by-3 magic square is known as Lo Shu. Legend has it that Lo Shu was discovered on the shell of a turtle that crawled out of the Lo River in the 23rd century b.c. Lo Shu provides a mathematical basis for feng shui, the ancient Chinese philosophy of balance and harmony. Matlab can generate Lo Shu with
A = magic(3)
which produces
A = 8 1 6 3 5 7 4 9 2
The command
sum(A)
sums the elements in each column to produce
15 15 15
The command
sum(A’)’
transposes the matrix, sums the columns of the transpose, and then transposes the results to produce the row sums
15 15 15
20 Chapter 1. Introduction to MATLAB
The command
sum(diag(A))
sums the main diagonal of A, which runs from upper left to lower right, to produce
15
The opposite diagonal, which runs from upper right to lower left, is less important in linear algebra, so finding its sum is a little trickier. One way to do it makes use of the function that “flips” a matrix “upside-down.”
sum(diag(flipud(A)))
produces
15
This verifies that A has equal row, column, and diagonal sums. Why is the magic sum equal to 15? The command
sum(1:9)
tells us that the sum of the integers from 1 to 9 is 45. If these integers are allocated to 3 columns with equal sums, that sum must be
sum(1:9)/
which is 15. There are eight possible ways to place a transparency on an overhead projec- tor. Similarly, there are eight magic squares of order three that are rotations and reflections of A. The statements
for k = 0: rot90(A,k) rot90(A’,k) end
display all eight of them.
8 1 6 8 3 4 3 5 7 1 5 9 4 9 2 6 7 2
6 7 2 4 9 2 1 5 9 3 5 7 8 3 4 8 1 6
2 9 4 2 7 6 7 5 3 9 5 1 6 1 8 4 3 8
