Chapter 1 ntroduction 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, an on-line manual from The Math Works is available. Select Help in the toolbar atop the MATLAB command window, then select MATLAB Help and Getting Started. A PDF version is vailable under Printable versions The document is also available from the Math Works Web site 10. Many other manuals produced by The Math Works are available on line and from the Web site A list of over 600 MATLAB-based books by other authors and publishers, in sev- eral languages, is available at [11. Three introductions to MATLAB are of particular interest here: a relatively short primer by Sigmon and Davis [8, a medium-sized mathematically oriented text by Higham and Higham 3, and a large, compreher sive 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 beer collected in a directory (or folder)named (The directory name is the initials of the book title. You can either start MATLAB in this directory or use patht to add the directory to the MATLAB path December 26 2005
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 investigate 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, an on-line manual from The MathWorks is available. Select Help in the toolbar atop the Matlab command window, then select MATLAB Help and Getting Started. A PDF version is available under Printable versions. The document is also available from The MathWorks Web site [10]. Many other manuals produced by The MathWorks are available on line and from the Web site. A list of over 600 Matlab-based books by other authors and publishers, in several languages, is available at [11]. Three introductions to Matlab are of particular interest here: a relatively short primer by Sigmon and Davis [8], a medium-sized, mathematically oriented text by Higham and Higham [3], and a large, comprehensive 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. December 26, 2005 1
Chapter 1. Introduction to MATLAB 1.1 The golden ratio What is the world's most interesting number? Perhaps you like T, or e, or 17. Some people might vote for o, the golden ratio, computed here by our first MATLAB phi =(1 sqrt(5))/2 Thi hi= 1.6180 Let's see more digits. format long phi 1.61803398874989 This didn't recompute it just displayed 15 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 haller rectangle with the same shap Figure l.l. The golden rectangle Equating the aspect ratios of the rectangles gives a defining equation for o This equation says that you can compute the reciprocal of o by simply subtracting one. How many numbers have that property? Multiplying the aspect ratio equation by o produces the polynomial equation
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))/2 This produces phi = 1.6180 Let’s see more digits. format long phi phi = 1.61803398874989 This didn’t recompute φ, it just displayed 15 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
1.1. The Golden Ratio The roots of this equation are given by the quadratic formula: 1士 The positive root is the golden ratio. If you have forgotten the quadratic formula, you can ask MaTLab to find the roots of the polynomial. MATLAB represents a polynomial by the vector of its coefficients, in descending order. So the vector p(a) The roots are computed by the roots function. r roots(p) 803398874989 1.61803398874989 ers are the only numbers whose reciprocal can be computed by subtracting You can use the Symbolic Toolbox, which connects MATLAB to Maple, to solve the aspect ratio equation without converting it to a polynomial. The equation is represented by a character string. The solve function finds two solutions r= solve(1/x=x-12) produce 1/2*5^(1/2)+1/2] [1/2-1/2*5(1/2)] e pretty function displays the results in a way that resembles typeset mathe- matics pretty (r) 1/2 [1/25+1/2] 1/2] [1/2-1/
1.1. The Golden Ratio 3 The roots of this equation are given by the quadratic formula: φ = 1 ± √ 5 2 . The positive root is the golden ratio. If you have forgotten the quadratic formula, you can ask Matlab to find the roots of the polynomial. Matlab represents a polynomial by the vector of its coefficients, in descending order. So the vector p = [1 -1 -1] represents the polynomial p(x) = x 2 − x − 1. The roots are computed by the roots function. r = roots(p) produces r = -0.61803398874989 1.61803398874989 These two numbers are the only numbers whose reciprocal can be computed by subtracting one. You can use the Symbolic Toolbox, which connects Matlab to Maple, to solve the aspect ratio equation without converting it to a polynomial. The equation is represented by a character string. The solve function finds two solutions. r = solve(’1/x = x-1’) produces r = [ 1/2*5^(1/2)+1/2] [ 1/2-1/2*5^(1/2)] The pretty function displays the results in a way that resembles typeset mathematics. pretty(r) produces [ 1/2 ] [1/2 5 + 1/2] [ ] [ 1/2] [1/2 - 1/2 5 ]
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 ph 1/2*5~(1/2)+1/2 TH ion can be converted to a numerical value in two different ways. It can to any number of digits using variable-precision arithmetic with the vpa(phi, 50) produces 50 digits 1.6180339887498948482045868343656381177203091798058 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) 1.61803398874989 The aspect ratio equation is simple enough to have closed-form symbolic so- lutions. More complicated equations have to be solved approximately. The inline function is a quick way to convert character strings to objects that can be arguments to the MATLAB functions that operate on other functions f inline(1/x -(x-1)2) defines f(a)=1/-(a-1)and produces Inline function f(x)=1/x-(x-1) (Beginning with MATLAB 7, inline objects will be superceeded by a more powerful construction known as anonymous functions. Inline objects are still allowed in MATLAB 7, but anonymous functions are preferred because they produce more efficient code. The graph of f(a)over the interval 0 <a<4 shown in Figure 1.2 is obtained ezplot(f, 0,4
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 = 1/2*5^(1/2)+1/2 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. 1.6180339887498948482045868343656381177203091798058 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 = 1.61803398874989 The aspect ratio equation is simple enough to have closed-form symbolic solutions. More complicated equations have to be solved approximately. The inline function is a quick way to convert character strings to objects that can be arguments to the Matlab functions that operate on other functions. f = inline(’1/x - (x-1)’); defines f(x) = 1/x − (x − 1) and produces f = Inline function: f(x) = 1/x - (x-1) (Beginning with Matlab 7, inline objects will be superceeded by a more powerful construction known as anonymous functions. Inline objects are still allowed in Matlab 7, but anonymous functions are preferred because they produce more efficient code.) The graph of f(x) over the interval 0 ≤ x ≤ 4 shown in Figure 1.2 is obtained with ezplot(f,0,4)
1.1. The Golden Ratio (x1) Figure 1.2. f(o)=0 The name ezplot stands for "easy plot, although some of the English-speaking world would pronounce it"e-zed plot " Even though f(a) becomes infinite as a-0, ezplot automatically picks a reasonable vertical scale The statement phi fzero(f, 1) looks for a zero of f(a)near a =1. It produces an approximation to o 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 goldrect 7 GOLDRECT Plot the golden rectangle hi=(1+sqrt(5))/2 to phi phi 0 oJ y [00110] u=[11]; [01] plot(x,y, 'b,,u,v,'b--) text (phi/2, 1.05, '\phi ')
1.1. The Golden Ratio 5 0 0.5 1 1.5 2 2.5 3 3.5 4 −3 −2 −1 0 1 2 3 4 5 6 7 x 1/x − (x−1) Figure 1.2. f(φ) = 0. 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]; plot(x,y,’b’,u,v,’b--’) text(phi/2,1.05,’\phi’)
Chapter 1. Introduction to MATLAB text(1+phi)/2,-.05,phi-1) text(-.05,5,1) text(.5,-.05,1) axis equa axis off set(gcf, 'color','white') The vectors x and y each contain five elements. Connecting consecutive (k, yk)pairs with straight lines produces the outside rectangle. The vectors u and v each contain two elements. The line connecting(u1, v1)with(u2, u2) sepa- ates the rectangle into the square and the smaller rectangle. The plot command draws these lines-the x-y lines in solid blue and the u-v line in dashed blue The next four statements place text at various points; the string"phi' denotes the Greek letter. The two axis statements cause the scaling in the z and y directions to be equal and then turn off the display of the axes. The last statement sets the background color of gcf, which stands for get current figure, to white A continued fraction is an infinite expression of the form a1+a2+a3+… If all the aks are equal to l, the continued fraction is another representation of the The following MATLAB function generates and evaluates truncated continued frac tion approximations to o. The code is stored in an M-file named goldfract m. function goldfract(n) ZGOLDFRACT Golden ratio continued fraction GOLDFRACT (n) displays n terms P 8b=[1+/(P")1 P end p= sprintf('%d/d',P, q
6 Chapter 1. Introduction to MATLAB text((1+phi)/2,-.05,’\phi - 1’) text(-.05,.5,’1’) text(.5,-.05,’1’) axis equal axis off set(gcf,’color’,’white’) The vectors x and y each contain five elements. Connecting consecutive (xk, yk) pairs with straight lines produces the outside rectangle. The vectors u and v each contain two elements. The line connecting (u1, v1) with (u2, v2) separates the rectangle into the square and the smaller rectangle. The plot command draws these lines—the x − y lines in solid blue and the u − v line in dashed blue. The next four statements place text at various points; the string ’\phi’ denotes the Greek letter. The two axis statements cause the scaling in the x and y directions to be equal and then turn off the display of the axes. The last statement sets the background color of gcf, which stands for get current figure, to white. A continued fraction is an infinite expression of the form a0 + 1 a1 + 1 a2+ 1 a3+··· . If all the ak’s are equal to 1, the continued fraction is another representation of the golden ratio: φ = 1 + 1 1 + 1 1+ 1 1+··· . The following Matlab function generates and evaluates truncated continued fraction approximations to φ. The code is stored in an M-file named goldfract.m. function goldfract(n) %GOLDFRACT Golden ratio continued fraction. % GOLDFRACT(n) displays n terms. p = ’1’; for k = 1:n p = [’1+1/(’ p ’)’]; end p p = 1; q = 1; for k = 1:n s = p; p = p + q; q = s; end p = sprintf(’%d/%d’,p,q)
11. The golden ratio format long p= eval(p) format short err =(1+sqrt(5))/2-p The statement oldfract(6) 1+1/(1+1/(1+1/(1+1/(1+1/(1+1/(1)))))) 21/13 1.61538461538462 e The three p's are all different representations of the same approximation to o 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'(tha goldfract(o) 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 1 So the iteration starts with and repeatedly replaces the fraction The statement p= sprintf('7d/7d',p, q)
1.1. The Golden Ratio 7 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/13 p = 1.61538461538462 err = 0.0026 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 1 + 1 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 p = sprintf(’%d/%d’,p,q)
Chapter 1. Introduction to MATLAB 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 ctually do the division expressed in the second p The final quantity err is the difference between p and o. 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 generate by goldfract(n) theoretically approaches o. But limitations on the size of the integers in the numerator and denominator, as well as roundoff error in the actual Hoating-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 abaca 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 would not have a sequence named after him. The number of pairs would simply double each month. After n months there would be 2" pairs of rabbits. Thats 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: The initial conditions are that in the first month there is one pair of rabbits and in he second there are two pairs f1=1,f2=2
8 Chapter 1. Introduction to MATLAB 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 transcribed 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: f1 = 1, f2 = 2
1.2. Fibonacci Numbers The following MATLAB function, stored in the M-file fibonacci. m, produces a vector containing the first n Fibonacci numbers function f fibonacci(n) Z FIBONACCI Fibonacci sequence f FIBONACCI(n) generates the first n Fibonacci numbers f zeros(n, 1) f(k)=f(k-1)+f(k-2) With these initial conditions, the answer to Fibonacci's original question about the size of the rabbit popu ter one year is given I fibonacci(12) 8 The answer is 233 pairs of rabbits.(It would be 4096 pairs if the number doubled every month for 12 months. Let's look carefully at fibonacci. m. It's a good example of how to create a MATLAB function. The first line is The first word on the first line says this is a function M-file, not a script. The remainder of the first line says this particular function produces one output result f, and takes one input argument, n. The name of the function specified on the first line is not actually used, because MATLAB looks for the name of the M-file, but it is common practice to have the two match. The next two lines are comments that provide the text displayed when you ask for help help fibonacci
1.2. Fibonacci Numbers 9 The following Matlab function, stored in the M-file fibonacci.m, produces a vector containing the first n Fibonacci numbers. function f = fibonacci(n) % FIBONACCI Fibonacci sequence % f = FIBONACCI(n) generates the first n Fibonacci numbers. f = zeros(n,1); f(1) = 1; f(2) = 2; for k = 3:n f(k) = f(k-1) + f(k-2); end With these initial conditions, the answer to Fibonacci’s original question about the size of the rabbit population after one year is given by fibonacci(12) This produces 1 2 3 5 8 13 21 34 55 89 144 233 The answer is 233 pairs of rabbits. (It would be 4096 pairs if the number doubled every month for 12 months.) Let’s look carefully at fibonacci.m. It’s a good example of how to create a Matlab function. The first line is function f = fibonacci(n) The first word on the first line says this is a function M-file, not a script. The remainder of the first line says this particular function produces one output result, f, and takes one input argument, n. The name of the function specified on the first line is not actually used, because Matlab looks for the name of the M-file, but it is common practice to have the two match. The next two lines are comments that provide the text displayed when you ask for help. help fibonacci produces
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 he first line is not displayed when you ask for help on the function he next line f zeros(n, 1) by-1 matrix all zeros and assigns it to f. In mAtlAB matrix with only one column is a column vector and a matrix with only one row is t two lir f(2)=2 provide the initial conditions The last three lines are the for statement that does all the work for k= 3: n (k-1)+f(k-2) 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 i FIBNUM(n) generates the nth Fibonacci number if n<=1 else f fibnum(n-1)+fibnum(n-2) end The statement fibnum(12) produces
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 languages. 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