Download Introduction to Numerical Methods - Notes | AMATH 352 and more Study notes Mathematics in PDF only on Docsity!
Applied Math 352 Notes by R. J. LeVeque
1 Introduction to Numerical Methods
Computers are quite stupid. They really don’t “know” how to do much of anything, in terms of what is built into the hardware. As far as mathematics is concerned, they are basically capable only of simple arithmetic: addition, subtraction, multiplication and division. Anything else that we want to do mathematically using a computer has to be eventually reduced to a (typically very long!) sequence of these simple operations. Remember that everything on a computer is represented as a string of zeros and ones. Any data, such as numbers, photographs, music, etc., are stored as a large set of zeros and ones, binary digits or bits for short. Computer storage consists of some sort of medium that has a huge number of locations, each of which can be set to a 0 value or a 1 value. How these are stored depends on what sort of medium it is. On a CD or DVD, for example, there are trillions of individual locations that can each be etched with a laser making a tiny dent (representing 1) or not (representing 0). Operations such as addition on a computer are performed by electronic components that take two electrical signals, each representing a sequence of 0’s and 1’s, and producing as output a third electrical signal that is somehow produced from the two inputs. If the two input signals represent two numbers expressed in binary notation, then it’s possible to design an “adder” that takes these two signals and produces an output signal that is the binary representation of the sum of the two numbers. It seems reasonable that such a device can be built, since there is a fairly simple algorithm for combining the bits in the summands to get the bits in the sum (corresponding bits are added using the rules 0 + 0 = 0, 0 + 1 = 1, and 1 + 1 = 0 with a 1 carried to the next place). With a bit more work, hardware implementing multiplication can also be designed. Any other mathematical operations that we want to perform on a computer must somehow be reduced to these operations. Consider the simple process of computing a square root, for example. You may think a computer or even the simplest calculator knows how to do this directly since there is a square root button on most calculators and every computer language has some built in function with a name like sqrt that returns a square root. But think about it: how would you design an electrical circuit that takes the bits representing a number and massages them directly to get the bits representing its square root? There is no simple bit-by-bit algorithm the way there is for addition or multiplication. When you press the square root button on your calculator or use sqrt in a computer program, you are actually invoking some other computer program the reduces the problem of computing a square root to a sequence of simple arithmetic operations (additions, multiplica- tions) that are then done by the hardware. There are many ways that this might be done and we’ll look at a couple in a moment. This is the sort of problem that we will be studying in this course: how do we take compli- cated mathematical problems and reduce them ultimately to a series of arithmetic operations that a binary computer can perform in a finite amount of time that will produce some reason- able approximation to the true solution of the problem. Of course we will mostly look at more challenging problems than computing the square root, and when we write programs we will
freely use built-in functions such as sqrt and hence build upon algorithms that are already available, but it’s useful to consider the square root problem as an illustration of some of the issues that arise when trying to develop an approximation algorithm.
1.1 Computing a square root
How would you compute
19, for example, if you were only allowed to use +, −, ∗, /? Here are two possible approaches....
First approach. We know that 4^2 = 16 and 5^2 = 25 so
19 must be between 4 and 5. Let’s try squaring the midpoint of this interval, 4.5. We find that (4.5)^2 = 20.25, so this is too big. Hence
19 must be between 4 and 4.5. Again try the midpoint: (4.25)^2 = 18.0625. This is too small, so
19 must lie between 4.25 and 4.5. Next we can try the midpoint of this interval and find that 4. 3752 = 19.1406, so this is too big. And so on... This is rather tedious but if we keep doing this long enough we will get down to a very small interval that contains√
This is an iterative method: We have some way of taking an approximation (an interval containing
19 in this case) and improving it at the next iteration to get a smaller interval. If we repeat this process many more times and print out the midpoint of the interval we have in each iteration as a single number approximating
19, we would find the following:
iteration midpoint 0 4. 1 4. 2 4. 3 4. 4 4. 5 4. 6 4. 7 4. 8 4. 9 4. 10 4. 11 4. 12 4. 13 4. 14 4. 15 4. 16 4. 17 4. 18 4. 19 4. 20 4.
For comparison, the correct value of
19 is
4.358898943540...
6 8 10 12 14 16 18 20
−0.
0
1
Figure 1.1: The definite integral (1.1) interpreted as an area.
is the matlab prompt, so anything following this represents something you can type into matlab at the prompt.
- A simple function f (x) can be specified using the syntax
f = @(x) 3*x + sin(x); for example. This defines f (x) = 3x + sin(x). This function can then be evaluated at a point or passed into another function, as will be done in the examples below.
1.3 Computing a definite integral.
Suppose we need to compute (^) ∫ 17
9
cos(
x + log(x)) dx. (1.1)
In calculus you learned how to compute definite integrals by finding an anti-derivative of the integrand and evaluating it at the endpoints of the interval. However, this only works if you can find an anti-derivative. The function cos(
x + log(x)) does not have an anti-derivative that can can be written down in terms of standard elementary functions, so this approach fails. However, the definite integral (1.1) does have a well-defined value. It is the area under the curve y = cos(
x + log(x)) for x between 9 and 17. This is sketched in Figure 1.1. The area of the shaded region is the integral, and the desired value is approximately
∫ (^17)
9
cos(
x + log(x)) dx = 6. 952058905443586... (1.2)
How might we approximate this? One approach is to split the interval [9, 17] into pieces of length h and construct a set of rectangles as shown in Figure 1.2(a) whose base has length h and whose height is given by evaluating f (x) at one edge of the piece. The sum of the areas of all these rectangles will approximate the area under the curve. In Figure 1.2(a) only four rectangles have been used and the approximation will not be very good, but by using smaller h and many more rectangles, an arbitrarily good approximation can be computed — in fact you probably saw the definite integral defined in calculus as the limit as h → 0 of this approximation.
6 8 10 12 14 16 18 20
−0.
0
1
6 8 10 12 14 16 18 20
−0.
0
1
Figure 1.2: (a) The definite integral (1.1) approximated by the areas of 4 rectangles. (b) The definite integral (1.1) approximated by the areas of 4 trapezoids.
Making h smaller (using more rectangles) requires more work: to compute the area of each rectangle we need to evaluate f (x) at one point, and so for every rectangle we require one function evaluation of the integrand. For some problems it can be quite expensive (i.e., require quite a bit of computer time) to evaluate f at each point. So we would like to develop methods that compute more accurate approximations with fewer function evaluations. Figure 1.2(b) shows a better approximation, obtained by approximating the area by a set of trapezoids instead of rectangles. The number of function evaluations is nearly the same (one more, at the right end of the interval), but the approximation obtained is typically much better. For the very coarse approximations shown in Figure 1.2, the rectangle and trapezoid approximations give:
rectangle = 6.516622749765546e+ trapezoid = 6.832431699945039e+
Increasing the number of intervals improves the approximations. For example, if we subdivide [9, 17] into 1000 subintervals instead of only 4, the rectangle and trapezoid approximations are
rectangle = 6.950792491769874e+ trapezoid = 6.952056992070894e+
In matlab there are a number of different built-in functions to approximate the values of integrals by numerical methods. This procedure is often called quadrature (motivated by the idea of replacing the area by a set of quadrilaterals such as rectangles or trapezoids), and the matlab routines have names like quad, quadl, quadv, etc. These implement various more sophisticated approaches to estimating definite integrals and we’ll look more closely at them later. They are very simple to use:
f = @(x) cos(sqrt(x) + log(x)); integral = quad(f, 9, 17) integral =
This is a linear system of 3 equations in 3 unknowns. If you have studied linear algebra before, you know that this system can be written in matrix-vector form as Ax = b where
A =
(^) , x =
x 1 x 2 x 3
(^) , b =
and that the solution can be computed exactly by performing Gaussian elimination on this system. Doing so, you should find that
x 1 = − 13 / 15 , x 2 = 7/ 5 , x 3 = − 14 / 15.
This problem is easily solved in matlab by using the backslash operator:
A = [3 4 0; 2 1 5; 7 8 -2] A = 3 4 0 2 1 5 7 8 - b = [3; -5; 7] b = 3
x = A\b x = -8.666666666666666e- 1.400000000000000e+ -9.333333333333335e-
This is a good approximation to the true solution. The errors, in the 16th digit, are due to the fact that the computer uses finite precision arithmetic, each number is stored using only a finite number of bits and there is no way to exactly store numbers like − 13 /15 or 14/15. For such a small system of equations, Gaussian elimination can be done by hand and the exact solution obtained. However, in many practical problems it is necessary to solve linear systems of equations with thousands or millions of equations or unknowns, and then it is necessary to use a computer implementation of Gaussian elimination. Solving such a system requires doing billions or trillions of arithmetic operations and one might worry about several things, such as:
- How long will this take on a computer? Will I live long enough to see the answer?
- Since computers store only approximations of real numbers and computer arithmetic is not exact, do the errors made in the process of doing Gaussian elimination accumulate to give such a large error that the result obtained is meaningless?
We will consider both these questions and others in detail later.
1.6 Finding the eigenvalues of a matrix
If you don’t know what eigenvalues are, don’t worry — we will learn what they are and how they are used in various contexts. A 3 × 3 matrix like the matrix A in (1.4) has 3 eigenvalues. In general these can not be found exactly, but there are many approaches to approximating them numerically. In matlab, we can find good approximations using the eig function:
lambda = eig(A) lambda = 8.582575694955835e+ -5.825756949558420e- -5.999999999999998e+
1.7 Solving differential equations
If you have taken a course on ordinary differential equations like Math 307 or AMath 351 then you learned many techniques for solving particular differential equations analytically. Already in Calculus you should have seen some differential equations, and know for example that the solution to y′(t) = 3y(t), y(0) = 2
is y(t) = 2e^3 t. Many differential equations that arise in practice cannot be solved exactly. Instead we can compute approximations to the solution at a large set of points using a numerical method. For example, the differential equation
y′(t) = cos(ty(t))
1 + y(t), y(0) = 2 (1.5)
cannot be solved exactly, but we can produce a plot of an approximate solution over the interval 0 ≤ t ≤ 4, for example, by:
% plotODEsolution.m f = @(t,y) cos(t*y) * sqrt(1+y); y0 = 2; ODEsolution = ode23(f,[0 4],y0); t = linspace(0, 4, 100); y = deval(ODEsolution, t); plot(t,y)
This gives the plot shown in Figure 1.3. These commands could all be typed in at the prompt. However, for things like this is it much better to create an m-file script. This is simply a text file with extension .m that contains the commands to be executed. Typing the name of this file (without the .m) at the prompt causes this sequence of commands to be executed. So if the above commands were in the file plotODEsolution.m then typing
plotODEsolution
step in gaining trust in something we have computed. So we will not only learn how to derive and implement numerical methods, we will also learn how to analyze them and estimate the error in the results they produce.
- For more complicated problems than the ones considered in this section it is often nec- essary to develop a new numerical method, or at least adapt a well-known method to a new context. In order to do this effectively, it is necessary to know a lot about existing methods for standard problems.
- Finally, it’s lots of fun to invent, analyze, and test new methods!
New numerical methods are constantly being invented and analyzed, there are dozens of technical journals and thousands of books devoted to this topic. Numerical methods are invented and studied by mathematicians who want to analyze the properties of methods (prove they converge, provide error estimates, etc.) and by scientists and engineers who want to solve real world problems in many different fields. Historically there have been two main sides to science: theory and experiments. Recently a third major discipline has been added to this: computational science, in which computations are used to test out theories and compare to experiments. To do so the theory is usually expressed in the form of mathematical equations that usually cannot be solved exactly. The rapid increase in computing power and speed has made it possible to do amazing computations that were inconceivable even a few years ago. But computers are still as stupid as ever and can only do basic arithmetic (incredibly quickly). Numerical methods for solving more complicated mathematical problems lie at the heart of this exciting development.
