Baixe Numerical Computing with MATLAB - Cleve Moler e outras Notas de estudo em PDF para Engenharia Elétrica, somente na Docsity!
Preface
Numerical Computing with MATLAB is a textbook for an introductory course in numerical methods, Matlab, and technical computing. The emphasis is on informed use of mathematical software. We want you to learn enough about the mathematical functions in MATLAB that you will be able to use them correctly, appreciate their limitations, and modify them when necessary to suit your own needs. The topics include:
- introduction to MATLAB
- linear equations
- interpolation
- zero finding
- least squares
- quadrature
- ordinary differential equations
- random numbers
- Fourier analysis
- eigenvalues and singular values
- partial differential equations
George Forsythe initiated a software-based numerical methods course at Stan- ford University in the late 1960s. The textbooks by Forsythe, Malcolm, and Moler [1] and Kahaner, Moler, and Nash [2] that evolved from the Stanford course were based upon libraries of Fortran subroutines. This textbook is based upon MATLAB. NCM, a collection of over 70 M- files, forms an essential part of the book. Many of the over 200 exercises involve modifying and extending the programs in NCM. The book also makes extensive use of computer graphics, including interactive graphical expositions of numerical algorithms. The prerequisites for the course, and the book, include:
2 Preface
- calculus
- some familiarity with ordinary differential equations
- some familiarity with matrices
- some computer programming experience
If you’ve never used Matlab before, the first chapter will help you get started. If you’re already familiar with Matlab you can glance over most of the first chapter quickly. Everyone should read the section in the first chapter about floating point arithmetic. There is probably too much material here for a one-quarter or one-semester course. Plan to cover the first several chapters and then choose the portions of the last four chapters that interest you. Make sure that the NCM collection is installed on your network or your per- sonal computer as you read the book. The software is available from a Web site devoted to the book,
http://www.mathworks.com/moler
There are three types of NCM files:
- gui files. Interactive graphical demonstrations.
- tx files. Textbook implementations of built-in Matlab functions.
- Others. Miscellaneous files, primarily associated with exercises.
When you have NCM available,
ncmgui
produces the figure shown on the next page. Each thumbnail plot is actually a push button that launches the corresponding gui. This book would not have been possible without the staff of The MathWorks. They are a terrific group of people and have been especially supportive of this book project. Out of the many friends and colleagues who have made specific contributions, I want to mention five in particular. Kathryn Ann Moler has used early drafts of the book several times in courses at Stanford and has been my best critic. Tim Davis and Charlie Van Loan wrote especially helpful reviews. Lisl Urban did an immaculate editing job. My wife Patsy has lived with my work habits and my laptop and loves me anyway. Thanks, everyone.
- Cleve Moler, January 5, 2004
4 Preface
Bibliography
[1] G. Forsythe, M. Malcolm, and C. Moler, Computer Methods for Math- ematical Computations, Prentice Hall, Englewood Cliffs, 1977.
[2] D. Kahaner, C. Moler, and S. Nash, Numerical Methods and Software, Prentice Hall, Englewood Cliffs, 1989.
[3] The MathWorks, Inc., Numerical Computing with MATLAB, http://www.mathworks.com/moler
2 Chapter 1. Introduction to 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 15 significant digits instead of five. 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 a polynomial equation
φ^2 − φ − 1 = 0
The roots of this equation are given by the quadratic formula.
φ =
1.1. The Golden Ratio 3
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.
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/25^(1/2)+1/2] [ 1/2-1/25^(1/2)]
The pretty function displays the results in a way that resembles typeset mathe- matics.
pretty(r)
produces
[ 1/2 ] [1/2 5 + 1/2] [ ] [ 1/2] [1/2 - 1/2 5 ]
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)
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
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’) text((1+phi)/2,-.05,’\phi - 1’) text(-.05,.5,’1’) text(.5,-.05,’1’) axis equal axis off set(gcf,’color’,’white’)
6 Chapter 1. Introduction to MATLAB
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 (u 1 , v 1 ) with (u 2 , v 2 ) sepa- rates 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
a 0 +
a 1 + (^) a 2 + 11 a 3 +...
If all the ak’s are equal to 1, the continued fraction is another representation of the golden ratio.
φ = 1 +
1+...
The following Matlab function generates and evaluates truncated continued frac- tion 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)
format long p = eval(p)
format short err = (1+sqrt(5))/2 - p
8 Chapter 1. Introduction to MATLAB
The final quantity err is the difference between p and φ. With only six terms, the approximation is accurate to less than three digits. How many terms does it take to get 10 digits of accuracy? As the number of terms n increases, the truncated continued fraction gener- ated 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. One of the exercises 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
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
1.2. Fibonacci Numbers 9
% 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 a small 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
FIBONACCI Fibonacci sequence f = FIBONACCI(n) generates the first n Fibonacci numbers.
1.2. Fibonacci Numbers 11
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
12 Chapter 1. Introduction to MATLAB
and ends with
Do you see why we chose n = 40? Use the up arrow key on your keyboard to bring back the previous expression. Change it to
f(2:n)./f(1:n-1) - phi
and then press the Enter key. What is the value of the last element? The population of Fibonacci’s rabbit pen doesn’t double every month; it is multiplied by the golden ratio every month. It is possible to find a closed-form solution to the Fibonacci number recurrence relation. The key is to look for solutions of the form
fn = cρn
for some constants c and ρ. The recurrence relation
fn = fn− 1 + fn− 2
becomes
ρ^2 = ρ + 1
We’ve seen this equation before. There are two possible values of ρ, namely φ and 1 − φ. The general solution to the recurrence is
fn = c 1 φn^ + c 2 (1 − φ)n
The constants c 1 and c 2 are determined by initial conditions, which are now conveniently written
f 0 = c 1 + c 2 = 1 f 1 = c 1 φ + c 2 (1 − φ) = 1
An exercise will ask you to use the Matlab backslash operator to solve this 2-by- system of simultaneous linear equations, but it is actually easier to solve the system by hand.
c 1 =
φ 2 φ − 1
c 2 = −
(1 − φ) 2 φ − 1
Inserting these in the general solution gives
fn =
2 φ − 1
(φn+1^ − (1 − φ)n+1)
14 Chapter 1. Introduction to MATLAB
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 0’s and 1’s for use with sparse matrix and image processing functions.
Figure 1.3. Fractal fern
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 also view the file with
F = imread(’fern.png’); image(F)
1.3. Fractal Fern 15
If you like the image, you might even choose to make it your computer desktop background. However, you should really run fern on your own computer to see the dynamics of the emerging fern in high resolution. The fern is generated by repeated transformations of a point in the plane. Let x be a vector with two components, x 1 and x 2 , representing the point. There are four different transformations, all of them of the form
x → Ax + b
with different matrices A and vectors b. These are known as affine transformations. The most frequently used transformation has
A =
, b =
This transformation shortens and rotates x a little bit, then adds 1.6 to its second component. Repeated application of this transformation moves the point up and to the right, heading towards the upper tip of the fern. Every once in a while, one of the other three transformations is picked at random. These transformations move the point into the lower subfern on the right, the lower subfern on the left, or the stem. Here is the complete fractal fern program.
function fern %FERN MATLAB implementation of the Fractal Fern % Michael Barnsley, Fractals Everywhere, Academic Press, % This version runs forever, or until stop is toggled. % See also: FINITEFERN.
shg clf reset set(gcf,’color’,’white’,’menubar’,’none’, ... ’numbertitle’,’off’,’name’,’Fractal Fern’) x = [.5; .5]; h = plot(x(1),x(2),’.’); darkgreen = [0 2/3 0]; set(h,’markersize’,1,’color’,darkgreen,’erasemode’,’none’); axis([-3 3 0 10]) axis off stop = uicontrol(’style’,’toggle’,’string’,’stop’, ... ’background’,’white’); drawnow
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];
