Lesson 4: Numerical Computations; Newton's method restart;

One case where roundoff error can be severe is if you subtract two numbers that are very close together: the relative error of the result can be much larger than the relative errors of the inputs. This is called "catastrophic cancellation". Here is an example of catastrophic cancellation affecting the solution of a quadratic equation. This is the solution formula from high school. s:= solve(a*x^2+b*x+c=0, x); Now suppose b is large compared to a and c. r:= solve(x^2 + 444*x + 1=0, x); Use evalf to get a floating-point value, with Digits = 3. Digits:= 3; evalf(r[2]); No problem here, even though we're using Digits = 3. evalf(r[2],20); But what about the first solution? evalf(r[1]); The relative error is 100%. evalf(r[1],20); Even at Digits = 10 the relative error is not insignificant: we only get 5 significant digits. evalf(r[1],10); But we can actually get better relative error at the same setting of "Digits" by using a different formula that won't involve this catastrophic cancellation. Recall that the product of the roots is the constant term c divided by the leading term a. expand(s[1]*s[2]); So the first root should be s[1] = c/a/s[2]; which won't produce a catastrophic cancellation (when b is positive). In our case: evalf(1/r[2]);

Another situation where roundoff error can be significant is in adding small numbers to large numbers. What happens when you add 1.23 and 456 using Digits = 3? 1.23 + 456; Note that the last two digits of the smaller number were basically lost: they had almost no influence on the final result. In a more extreme case, where the two numbers differ by at least Digits orders of magnitude, the smaller number would be lost altogether. 0.123 + 456; This can have serious consequences if you are trying to add up many numbers that are each very small (e.g. in numerical integration). Suppose we start with 0 and add 1000 numbers, each equal to .00444. The correct answer would be 4.44. And it would be, if we used Digits = 10. Digits:= 10: total := 0: for count from 1 to 1000 do total := total + .00444 end do: total; This is what's called a "for loop" (we'll talk about them in more detail later). It starts out with total = 0, and repeats 1000 times the command total := total + .00444, i.e. adding .00444 to the current total and making that the new value of total. But now try it with Digits = 3. Digits:= 3: total := 0: for count from 1 to 1000 do total := total + .00444 end do: total; What's happening? After a certain number of iterations total reaches the value 1.00. After that, adding .00444 doesn't change the result. 1.00+.00444; Even with Digits = 10, you could still have a problem if the numbers being added were small enough. There are some rather clever ways to get around the problem (which we won't go into), but I just wanted to point it out. The main moral of the story is that, especially in a situation where you can see that one of these problems may occur, you should use a large enough setting of Digits. People sometimes say "I only need this result to 3 significant digits" or "my data come from measurements that are only correct to 3 significant digits, so the results should only be meaningful to 3 significant digits" and are tempted to reduce Digits to 3. That's usually a bad idea. Never use fewer than the default Digits = 10. When you come to the end of your calculation, you can always round the final answer to 3 significant digits if you want. If you wish you can also go to Options, Precision, and tell Maple to "Round screen display to" a certain number of digits. This will affect the numbers you see in Maple output, without changing the numbers that Maple uses internally. For serious numerical work, use Digits = 14 or 15, or even more if roundoff error appears to be a problem. The only real problem with large values of Digits is that the computations will require more time and the data may take up more memory. evalf(2.3456787, 3); Digits:= 10;

Up to now we have been working with expressions in Maple. A mathematical expression is something that might be constructed from constants and variables using the algebraic operators +-*/^ and mathematical functions. For example, here we assign an expression as the value of A: A := sin(x) + x^3 * y/Pi; If we want to substitute values for the variables x and y, we can do so using eval. eval(A, {x = 2, y = 1/4}); And this expression, which has no variables left in it, represents a number. We can use evalf to get its approximate floating-point value. evalf(%); A function, on the other hand, can be thought of as a procedure that takes, let's say, a number (or an expression that might represent a number), does something to it, and returns another number (or expression representing a number). Maple's built-in mathematical functions include, for example, sin and exp. But you can also define your own functions. The simplest way to do so is with -> (that's "minus" and "greater than", which together make something that looks like an arrow pointing to the right. Thus: f := x -> sin(x) + x/Pi; You can then evaluate the function at a particular value of x by putting the value in parentheses after the name f, just as you would for sin or exp. f(3); f(x); f(t); You can also have a function of several variables: g := (x,y) -> sin(x) + x^3 * y/Pi; g(2, 1/4); g(u,v-t); An expression is tied to the particular names of the variables in it. If you assigned a value to x or y, that will affect the value of A. x := 5; A; But the x in the definition of the function f or g is just a "dummy variable", and these functions are not affected. f(t); g(u,v); If you just look at the function f itself, you won't see the definition. f; In order to see what the definition is, you can use eval. eval(f); It's often more convenient to use functions rather than expressions, especially when you might want to look at the function at many different values of the variables.

restart; Suppose you have an expression such as B := sin(x) - x/Pi; and you would like to make it into a function. If it's a complicated expression, you probably don't want to type it in again. It's tempting to try f:= x -> B; But that doesn't work. f(3); Why doesn't it work? Because there are really two different variables called x here: the one in B, which is an ordinary variable (called a global variable in Maple terminology), that you can assign values to, while the one in the definition of f is a dummy variable or "formal parameter" that would be unaffected by assigning values to the global x. When you enter the command f(3); Maple substitutes 3 for the formal parameter x in the definition of f, but doesn't do anything to the global variable x. The right way to make an expression into a function, is using the unapply command. Thus: f := unapply(B, x); f(3); You can also use unapply to make a function of several variables. A:= sin(x+y) - exp(y+z); g:=unapply(A,x,y,z); Caution: you want to define f, not f(x). Here is a common mistake: h(x) := sin(x) - x/Pi; It looks OK, Maple doesn't complain, and you can even enter h(x); But what this definition did was literally to define h(x), not h(anything else). h(y); h(3); Maple doesn't know what h(y) is supposed to be, it just knows h(x). What you should have done was h := x -> sin(x) - x/Pi; h(y);

There are two methods of differentiation in Maple: D, which applies to functions, and diff, which applies to expressions. Suppose we have an expression, such as B, containing the variable x: B := sin(x) - x/Pi; f:= unapply(B,x); We can use diff to take its derivative with respect to x. diff(B,x); If the expression contains more than one variable, what diff gives you is the partial derivative with respect to this particular variable. A:= sin(x+y) - exp(y+z); g:= unapply(A,x,y,z); diff(A, x); diff(A, y); On the other hand, D takes a function of one variable and gives you another function of one variable, its derivative: D(f); You can treat D(f) like any other function. For example, you can evaluate it at a point. D(f)(Pi); Here I'm going to look at both types of derivative. On the left I'm going to take the function D(f) and evaluate it at x. On the right I'll take the derivative of the expression f(x) with respect to x. The results should be the same. D(f)(x) = diff(f(x),x); To evaluate a diff expression at a certain value of x, you need to use eval: D(f)(Pi) = eval(diff(f(x),x), x=Pi); D on a function of more than one variable gives an error. D(g); Instead, use D[i] for the partial derivative with respect to the i'th variable. D[1](g)(x,y,z) = diff(g(x,y,z), x); D[2](g)(x,y,z) = diff(g(x,y,z), y); D[4](g)(x,y,z);

When fsolve finds an approximate solution to an equation, it's using a numerical method. The prototype of these methods for solving equations is Newton's method. The methods fsolve uses are sometimes slightly more sophisticated, but still based on Newton's method, so it's worthwhile trying to understand Newton's method. Newton's method is a method of approximately solving an equation, say LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2J0YrLUYjNiYtRiw2JVEiZkYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi1RMCZBcHBseUZ1bmN0aW9uO0YnL0Y6USdub3JtYWxGJy8lJmZlbmNlR1EmZmFsc2VGJy8lKnNlcGFyYXRvckdGRC8lKXN0cmV0Y2h5R0ZELyUqc3ltbWV0cmljR0ZELyUobGFyZ2VvcEdGRC8lLm1vdmFibGVsaW1pdHNHRkQvJSdhY2NlbnRHRkQvJSdsc3BhY2VHUSYwLjBlbUYnLyUncnNwYWNlR0ZTLUkobWZlbmNlZEdGJDYkLUYjNiQtRiw2JVEieEYnRjZGOUZARkBGQC1GPTYtUSI9RidGQEZCRkVGR0ZJRktGTUZPL0ZSUSwwLjI3Nzc3NzhlbUYnL0ZVRlxvLUkjbW5HRiQ2JFEiMEYnRkBGQEYrRkA=. We start with an initial guess LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEieEYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYkLUkjbW5HRiQ2JFEiMEYnL0Y2USdub3JtYWxGJ0Y+LyUvc3Vic2NyaXB0c2hpZnRHRj1GPg==, and Newton's method produces a sequence of numbers LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEieEYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYkLUkjbW5HRiQ2JFEiMUYnL0Y2USdub3JtYWxGJ0Y+LyUvc3Vic2NyaXB0c2hpZnRHUSIwRidGPg==, LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYkLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEieEYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYkLUkjbW5HRiQ2JFEiMkYnL0Y2USdub3JtYWxGJ0Y+LyUvc3Vic2NyaXB0c2hpZnRHUSIwRidGPg==, ... that converges very rapidly to a solution (when things are working well). As we'll see, it doesn't always work well, and we'll investigate what can happen. The formula for Newton's method is LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2JkYrLUYjNictSSVtc3ViR0YkNiUtRiw2JVEieEYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYmRistRiM2Ji1GLDYlUSJuRidGOUY8LUkjbW9HRiQ2LVEiK0YnL0Y9USdub3JtYWxGJy8lJmZlbmNlR1EmZmFsc2VGJy8lKnNlcGFyYXRvckdGTi8lKXN0cmV0Y2h5R0ZOLyUqc3ltbWV0cmljR0ZOLyUobGFyZ2VvcEdGTi8lLm1vdmFibGVsaW1pdHNHRk4vJSdhY2NlbnRHRk4vJSdsc3BhY2VHUSwwLjIyMjIyMjJlbUYnLyUncnNwYWNlR0Znbi1JI21uR0YkNiRRIjFGJ0ZKRkpGK0ZKLyUvc3Vic2NyaXB0c2hpZnRHUSIwRictRkc2LVEiPUYnRkpGTEZPRlFGU0ZVRldGWS9GZm5RLDAuMjc3Nzc3OGVtRicvRmluRmVvLUYjNiYtRjQ2JUY2LUYjNiRGQ0ZKRl5vLUZHNi1RKCZtaW51cztGJ0ZKRkxGT0ZRRlNGVUZXRllGZW5GaG4tSSZtZnJhY0dGJDYoLUYjNiZGKy1GIzYmLUYsNiVRImZGJ0Y5RjwtRkc2LVEwJkFwcGx5RnVuY3Rpb247RidGSkZMRk9GUUZTRlVGV0ZZL0ZmblEmMC4wZW1GJy9GaW5GXnEtSShtZmVuY2VkR0YkNiQtRiM2JEZpb0ZKRkpGSkYrRkotRiM2JkYrLUYjNidGKy1GIzYmLUYsNiVRIkRGJy9GOkZORkpGanAtRmFxNiQtRiM2JEZncEZKRkpGSkZqcEZgcUZKRitGSi8lLmxpbmV0aGlja25lc3NHRl1vLyUrZGVub21hbGlnbkdRJ2NlbnRlckYnLyUpbnVtYWxpZ25HRmdyLyUpYmV2ZWxsZWRHRk5GSkYrRkpGK0ZKRitGSg== It's convenient to define a function I'll call newt to do this, so that this is LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYqLUklbXN1YkdGJDYlLUkjbWlHRiQ2JVEieEYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYmLUYvNiVRIm5GJ0YyRjUtSSNtb0dGJDYtUSIrRicvRjZRJ25vcm1hbEYnLyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR0ZFLyUpc3RyZXRjaHlHRkUvJSpzeW1tZXRyaWNHRkUvJShsYXJnZW9wR0ZFLyUubW92YWJsZWxpbWl0c0dGRS8lJ2FjY2VudEdGRS8lJ2xzcGFjZUdRLDAuMjIyMjIyMmVtRicvJSdyc3BhY2VHRlQtSSNtbkdGJDYkUSIxRidGQUZBLyUvc3Vic2NyaXB0c2hpZnRHUSIwRictRj42LVEifkYnRkFGQ0ZGRkhGSkZMRk5GUC9GU1EmMC4wZW1GJy9GVkZcby1GPjYtUSI9RidGQUZDRkZGSEZKRkxGTkZQL0ZTUSwwLjI3Nzc3NzhlbUYnL0ZWRmJvRmhuLUYvNiVRJW5ld3RGJ0YyRjUtSShtZmVuY2VkR0YkNiQtRiM2JC1GLDYlRi4tRiM2JEY6RkFGZW5GQUZBLUY+Ni1RIi5GJ0ZBRkNGRkZIRkpGTEZORlBGW29GXW9GQQ== newt := x -> evalf(x - f(x)/D(f)(x)); There are a couple of things to notice here. I used evalf, because I will want decimal numbers here. I don't want complicated exact expressions. 3 - f(3)/D(f)(3); % - f(%)/D(f)(%); The definition of newt uses f itself, not the definition of f, so in the definition, even if f has already been defined as, say, LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYpLUkjbWlHRiQ2JVEieEYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi1RJiZtYXA7RicvRjNRJ25vcm1hbEYnLyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR0Y9LyUpc3RyZXRjaHlHRjEvJSpzeW1tZXRyaWNHRj0vJShsYXJnZW9wR0Y9LyUubW92YWJsZWxpbWl0c0dGPS8lJ2FjY2VudEdGPS8lJ2xzcGFjZUdRLDAuMjc3Nzc3OGVtRicvJSdyc3BhY2VHRkwtRiw2JVEkc2luRicvRjBGPUY5LUkobWZlbmNlZEdGJDYkRitGOS1GNjYtUSomdW1pbnVzMDtGJ0Y5RjtGPi9GQUY9RkJGREZGRkgvRktRLDAuMjIyMjIyMmVtRicvRk5GZW4tSSZtZnJhY0dGJDYoRistRiw2JVElJnBpO0YnRlJGOS8lLmxpbmV0aGlja25lc3NHUSIxRicvJStkZW5vbWFsaWduR1EnY2VudGVyRicvJSludW1hbGlnbkdGYm8vJSliZXZlbGxlZEdGPUY5, you see f(x) rather than LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYnLUkjbWlHRiQ2JVEkc2luRicvJSdpdGFsaWNHUSZmYWxzZUYnLyUsbWF0aHZhcmlhbnRHUSdub3JtYWxGJy1JKG1mZW5jZWRHRiQ2JC1GIzYkLUYsNiVRInhGJy9GMFEldHJ1ZUYnL0YzUSdpdGFsaWNGJ0YyRjItSSNtb0dGJDYtUSomdW1pbnVzMDtGJ0YyLyUmZmVuY2VHRjEvJSpzZXBhcmF0b3JHRjEvJSlzdHJldGNoeUdGMS8lKnN5bW1ldHJpY0dGMS8lKGxhcmdlb3BHRjEvJS5tb3ZhYmxlbGltaXRzR0YxLyUnYWNjZW50R0YxLyUnbHNwYWNlR1EsMC4yMjIyMjIyZW1GJy8lJ3JzcGFjZUdGVS1JJm1mcmFjR0YkNihGOC1GIzYkLUYsNiVRJSZwaTtGJ0YvRjJGMi8lLmxpbmV0aGlja25lc3NHUSIxRicvJStkZW5vbWFsaWduR1EnY2VudGVyRicvJSludW1hbGlnbkdGX28vJSliZXZlbGxlZEdGMUYy. When you use newt, it will use whatever is the current definition of f. If you changed f, newt would still be good for doing Newton's method with the new version of f. Let's start with a typical initial guess, and see where Newton's method takes us. f:= x -> sin(x) - x/Pi; x[0] := 3; This is actually an entry in what Maple calls a "table''. If you've done some computer programming, you may be pleasantly surprised to see that you don't need to tell Maple that you're making a table, or how big it will be: the table is created automatically whenever you first use "x[...]'', and will hold whatever entries you put in it. The index inside the brackets is shown as a subscript in output. x[1] := newt(x[0]); x[2] := newt(x[1]); This could take a while if I had to type this in each time I want to get another LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2JC1JJW1zdWJHRiQ2JS1GLDYlUSJ4RicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUYjNiQtRiw2JVEibkYnRjdGOi9GO1Enbm9ybWFsRicvJS9zdWJzY3JpcHRzaGlmdEdRIjBGJ0ZCRitGQg==. Fortunately, we can automate the process using a "for loop". for count from 2 to 6 do x[count + 1] := newt(x[count]) end do;

for ... from ... to ... do ... end do -> unapply diff D