Chapter8.mws

• Chapter 8: Programming -- logic, loops, and procedures
  • Logic: if-then-else-end if
  • Loops: for, do, end do, while, and break
  • Procedures: putting loops and logic to work
  • Example: plotting results from fsolve

---

Chapter 8: Programming -- logic, loops, and procedures

Sometimes Maple won't hav e a command that does what you need to do. When this happens you have to do your own programming in Maple using logic, loops, and procedures. We will start with logic

Logic: if-then-else-end if

Debugging

The basic logic command in Maple is the if statement. To illustrate it, suppose that you wanted a Maple function f(x) that would accept a value x and return -1 if x<0, 0 if x=0, and +1 if x>0. You would write it this way:

> restart;

> f:=x-> if x<0 then -1 elif x=0 then 0 else 1 end if;

Here's a test to see if it works

> f(-.5);f(0);f(.5);

(Note: this is the same as Maple's signum function.)

(As you will see in a minute, this i f-then-elif-else-end if way of building functions with special cases doesn't work well in Maple. The piecewise syntax introduced below is the best syntax to use. I am only using this more awkward form to introduce you to logic statements.)

>

The statement defining the function f almost reads like English if you just translate elif as "else if" and end if as "we're done". As you can see there are three elements that go into this statement.

The first are the logical syntax commands if ... then, elif ... then, else, and end if . They have to be in this order.

The second element is the conditional statement (logical test statement). These are things like x<0, etc.. Maple allows the following logical tests

< less than

<= less than or equal

> greater than

>= greater than or equal

= equal

<> not equal

When Maple encounters a test like this it evaluates it as a Boolean expression, i.e., as an expression whose value can be true or false. To help you debug things Maple has a a Boolean evaluator command evalb . It works like this

> evalb(7>5);evalb(3=4);evalb(3<>4);

Stare at the results of executing these commands until they make sense.

Note: and , or , and not can also be used to combine logic statements, like this

-2<=x and x<=-1 .

There are other conjuctions as well; see ?if and ?boolean for more details.

The third element of the if statement is the set of statements that it executes. In the example above these are just the numbers -1, 0, and 1, but any valid expression could go in these positions. For instance, we could use this construct to define a function which is one thing for negative numbers and another for positive numbers, like this:

> g:=x-> if x<0 then 1/(1+x^2) else cos(x) end if;

Test it to make sure there are no errors

> g(-1.5);g(2.5);

Ok, we seem to be getting numbers out. But watch what happens if you do this

> g(Pi);

Error, (in g) cannot evaluate boolean: Pi < 0

You need to able to interpret this error message because you will be getting it occasionally. It tells you that it tried to evaluate g, but that when it did, at least one of the Boolean expressions couldn't be evaluated. This happened because Pi is not really a number in Maple, so we just asked for a comparison between the number 0 and the non-number Pi. (Pi is more of an idea than a number in Maple.) You get the same error if you try to evaluate g(a), where a hasn't been given a number yet

> g(a);

Error, (in g) cannot evaluate boolean: a < 0

The problem with Pi can be fixed by remembering that Pi is tricky and always using evalf(Pi)

when you want a number:

> g(evalf(Pi));

Now let's plot our function and see what it looks like

> plot(g(x),x=-10..10);

Error, (in g) cannot evaluate boolean: x < 0

Whoa, bad booleans again. Now what went wrong? Well, with the if statement in g(x) Maple can't plot with this syntax. When x goes from the plot command into g, x is apparently not yet a number. A clue to solving this problem is found by looking at what Maple did to our statement g:=x->... above; it turned it into a procedure (more about these later). And Maple plots procedures using "operator notation", meaning that you don't say g(x),x=-10..10 in the plot command, you just leave x out, like this

> plot(g,-10..10);

But if you try to integrate the function or differentiate it, you are at a dead end:

> int(g(x),x=-2..2);

Error, (in g) cannot evaluate boolean: x < 0

> diff(g(x),x);

Error, (in g) cannot evaluate boolean: x < 0

This difficulty is the reason that the if-then-else-end if syntax is not very convenient for defining functions in Maple. But piecewise functions like this are pretty important, so it would be good if there were some clean way to handle them. Fortunately, there is. Maple has a piecewise command that is used for exactly this purpose. To redefine the function g(x) using piecewise we would do this

> g:=x->piecewise(x<0, 1/(1+x^2),x>=0, cos(x));

Now it really is a function and you can plot it like a function

> plot(g(x),x=-5..5);

integrate it like a function

> int(g(x),x=-5..5);

and differentiate it like a function

> diff(g(x),x);

So this is the best way to do piecewise functions; if-then-elif-end if is mostly used inside loops and procedures, which we will get to shortly.

>

Problem 8.1

Use the piecewise construct to define a Maple function which evaluates the so-called cubic B-spline function and plot it from t=-4..4. The function is defined this way:

_____________________________________________________________

if is in the range -2..-1

______________________________________________________________

if is in the range -1..0

B(t) = ______________________________________________________________

if is in the range 0..1

______________________________________________________________

- if is in the range 1..2

______________________________________________________________

For outside of the range -2..2 , .

You will discover that when you need to specify a range with a beginning and an end, like -1 to 0, you will need to combine two inequalities, i.e. b<a<c. Maple does this by using two different inequality conditions and combining them with and, like this:

b<a and a<c . When you get the function defined, plot it between -4 and 4.

----------------------------------------------------------------------------

Problem 8.2

Integrate B(t) from -3..3. Also plot its first three derivatives on the interval -3..3.

----------------------------------------------------------------------------

Problem 8.3

Here's an example of a problem where piecewise won't work cleanly and it is better to use the if...end if syntax.

Use this syntax to make a plot of the sine function from , but with the peaks and valleys all chopped off flat at height 0.9 or depth -0.9. The Maple command signum will be useful in doing this, as will your old friend evalf . Note that instead of doing logic horizontally (in x), as we do when using piecewise , you will be doing logic vertically (in y) instead.

Go to top of section

Loops: for, do, end do, while, and break

Debugging

Much of computing is done in loops: Maple's sum command has a loop at the bottom, as does numerical integration, fsolve , numerical differential equation solving, dsolve(...type=numeric) , and many more. You have seen an example of a homemade loop in the three-dimensional integration procedure in the calculus section[ Triple integrals] and in the differential equation section on Shooting. When Maple won't do what you want it to do, you will often find yourself writing your own loop. So that we have something concrete to talk about, let's write our own loop to sum the squares of the integers from 1 to 100. Maple's sum command does it this way

> sum(i^2,i=1..100);

Somewhere down in Maple's depths, what is actually being done by this command is a loop. A handmade one would go like this.

Sum Loop:

First you decide where you want the answer to be stored and set that variable to zero:

> answer:=0;

Then you make a loop that counts from 1 to 100, and inside that loop you square the counter and add its square into answer, like this using Maple's for..from..to..do..end do syntax. Note that do means "beginning of the loop" and end do means "end of the loop".

> for i from 1 to 100 do

> answer:=answer + i^2;

> end do;

Well, that was more information than you wanted. To fix this, end the end do that terminates the loop with a colon instead of a semicolon then look at answer at the end.

> answer:=0;

> for i from 1 to 100 do

> answer:=answer + i^2;

> end do:

> answer;

And what if you want to sum the odd numbers from 1 to 99? Well, you just want to start at 1, end at 99, and step up by 2. A loop that does this using Maple's by syntax is

> answer:=0;

> for i from 1 by 2 to 99 do

> answer:=answer + i^2;

> end do:

> answer;

Successive Substitution Loop:

If you know how far you want to count this is all you have to do, but sometimes you don't know how many times you want to loop. Instead, you want to loop until some condition is met. For instance, here is an interesting little trick that some of you may have noticed on your hand calculators. What happens if you just take the cosine of a number over and over and over? Let's let Maple do it 20 times and show us what it's doing as it goes along.

> x:=5.;

> for i from 1 to 20 do

Take the cosine of x and store the result in x so we can do it again and again

> x:=cos(x);

> end do;

It looks like the answer is converging to something. And in fact, if you think for just a minute, you will be able to see that the number it is converging to has to be a solution of the equation . (Explain to your partner or to the TA why this is true.) But you can also see that 20 times was not enough to get 10 digits right; what we would really like to do is loop until successive answers are within some tolerance of each other, say 1e-8, just to make the problem definite. Maple will let you end the loop that way too, with just a little bit of extra programming.

Choose a starting value of x and an initial fake value of xold, the "previous value of x"

> x:=5.;xold:=0;

Keep looping as long as |x-xold|>1e-8 (this is why we need the initial fake value of xold)

> for i from 1 while abs(x-xold)>1e-8 do

Put the new x from taking the cosine into a temporary variable tempx

> tempx:=cos(x):

Put the current x into xold so we will be able to tell when to quit

> xold:=x:

Put the new x (called t right now) into x

> x:=tempx;

Go back up and do it again

> end do;

This is called successive iteration, and it is used in computing a lot. It would be a good idea for you to look at the structure of this loop carefully to see how it is built. In particular, pay attention to when each variable is computed and where it is stored. You may also notice that I tried in the loop to have only x print by using colons on the t and xold lines; this didn't work and I don't know how to fix it.

>

Well, those are the two basic kinds of loops. For additional information see ?for.

Problem 8.4

Write a loop that evaluates the power series for the cosine function, for even . Give a value outside the loop, use the by syntax to get even numbers, and keep looping until the next term added to the sum is less than 1e-8. To do this you will need to know how to make a loop that adds. There is one at the beginning of this section, right here: Adding loop

--------------------------------------------------------------------------------

Problem 8.5

Write a loop that evaluates the Fourier series . Give a value outside the loop and keep looping until the magnitude of the terms being added drop below 1e-4 in magnitude. Since the sine function is never greater than 1 in magnitude, this means that you only need to apply the while test to . (It is bad to test on the whole term with the sine included because the sine might just happen to be really small and cause the loop to terminate too soon.)

--------------------------------------------------------------------------------

Problem 8.6

Write a loop to solve the equation by successive substitution. Get the answer to 8 significant figures. To see how to build a successive substitution loop go here (just above): Successive substitution loop

--------------------------------------------------------------------------------

Problem 8.7 (Secant Method)

Here's a loop that computational physicists use all the time. It's called the secant method and it's a way of solving hard equations in one variable. In fact, Maple almost certainly uses a fancy varition of this technique inside the fsolve command. I also used it in the shooting procedure for differential equations in Chapter 7. Here's the idea.

Suppose you want to solve a hard equation of the form , for . The secant method starts by asking you for two reasonably close guesses and . Then you evaluate the function at each value of to get and . You learned in high school algebra what to do with two ordered pairs ( ) and ( ) like this: use the two-point line formula to put a straight line through them: . But what we have just built is a straight-line approximation to the function , so to get an approximation to the value of that solves our hard equation, we just solve the line formula above for the value of that gives and call it , our next best guess at where the solution is.

> restart;

> Eqline:=0-f2=(f2-f1)/(x2-x1)*(x-x2);

> x3:=solve(Eqline,x);

>

This is a good algebraic expression, but not a good numerical one because as we get close to the root when both x1 and x2, and f1 and f2, are close to each other we get x3 by almost dividing zero by zero. This can cause floating point troubles, so it is better to rearrange this formula so that it looks like this:

.

With this form we get the next best guess by adding something small to the previous best guess and we get in less numerical trouble. It is true that both ( ) and ( ) are small, so we are dividing zero by zero again, but the values of and in the numerator are extra small because we are close to , so the sins we commit in dividing zero by zero don't affect the answer as much.

Ok, here is finally the problem: Write a loop that uses the secant method to solve the equation for a root near . (Make your loop similar to a Successive substitution loop.) In the loop you will find that it is desirable to jump out of the loop if to as many significant digits as Maple is keeping (typically 10; see online help about Digits ) because this will keep you from dividing by zero. This is hard to do with while , but fortunately there is a way to jump out using the break command.

Here is a stupid example: suppose you were looping from 1 to 100 but you wanted to jump out when i was 20. You would use

> for i from 1 to 100 do

> if i=20 then break end if;

> end do;

Notice that after the jump, i has the value it had when the jump occurred.

> i;

>

You will want to use this in the secant method procedure you are developing here. To test your use of it in your loop, keep repeating the calculation until the value of f2 is less than 1e-10 by using while ; this should make division by zero happen so you will need the break command to protect your loop.

You may discover, however, that getting Maple to divide by zero in this problem is difficult. If it won't cooperate, don't worry about it. The break statement is still a good idea in general, though, so I want you to practice putting it in.

Go to top of section

Procedures: putting loops and logic to work

Debugging

Well, you have probably noticed that loops and logic are messy. They take up a lot of space in your worksheet and to do them over and over you have to keep moving the cursor around. Maple has provided a way to put complicated stuff like this into something that works like a function so you can just pass values into it and get answers back out. The structure that does this is the procedure, and you have already encountered some of these in earlier sections. Now you get to learn how to write your own. Here is the basic structure.

> restart;

Assign the procedure to a name and define an argument list for it

> f:=proc(x,y)

Declare the names of local variables, those only used within the procedure

> local a,b,z;

Declare the names of global variables, those from outside to be used inside

> global d;

Now come the Maple statements that get the job done. In this case we use local a and b, global d, and passed in x and y to compute function value z, which is local. So how does Maple know what value to give us back? The returned value is always the last thing computed in the procedure.

> a:=3.;b:=17.;

> z:=a*b*x*y*d; # last value

> end;

> d:=31.;

> f(1,2);

>

To illustrate how a procedure works, here is one that simply evaluates the cosine of an angle in degrees.

> restart;

> cosdeg:=proc(theta)

> local phi,answer;

> phi:=theta*Pi/180.;

> answer:=cos(phi);

> end;

Now test it

> yy:=cosdeg(45.);

Because a procedure is a program that you write, it will usually fail. (I am not trying to insult you. This sentence is also true: "When scientists and engineers write procedures, they usually fail.") And when it fails you will need to see what it is doing to find out how to fix it. The Maple command to tell a procedure to reveal to you its inner workings is the command trace . To debug the procedure above use

> trace(cosdeg);

then execute it like this

> cosdeg(45);

{--> enter cosdeg, args = 45

<-- exit cosdeg (now at top level) = cos(.2500000000*Pi)}

And then when it is working right and you don't want to see all of this stuff every time you use it type

> untrace(cosdeg);

>

Problem 8.8

Notice also that you didn't get what you probably expected out of this cosdeg procedure. Fix it so that it gives a floating point number back. Make sure that your fix takes place inside the procedure and not in the call you make to it.

--------------------------------------------------------------------------------

If you look at how cosdeg is used you can see that a procedure is used just like a function. But because it is a self-contained little program it can be a really complicated function. And this is, in fact, the usual reason for writing a procedure. You have to be careful, though, because functions defined this way may not be compatible with other Maple commands like plot , int , sum , dsolve , etc.. In addition, procedures are usually slower than the built-in Maple commands.

Problem 8.9

Write a procedure that computes the cubic B-spline function described in Problem 8.1. Plot your procedure (function) from -3..3. Note: if you are going to plot the result of a procedure the answer has to be returned as the last result, and the only arguments that can be passed in are the appropriate independent variables, i.e. for a one-dimensional function you would use B:=proc(x) , for 2-d use C:=proc(x,y) , etc.. When you plot output from a procedure, you can treat it like a function using the plot(f(x),x=a..b) command, or you can use Maple's operator notation with the name of x left out: plot(B,-3..3) .

--------------------------------------------------------------------------------

Problem 8.10

You can also use procedures to perform a task that you want to do a whole bunch of times while you change some parameter. Write a procedure that takes as inputs the parameter Nterms and the angle and returns the value of the Fourier series given in Problem 8.5. Use Maple's sum command inside the procedure and keep terms up to n=Nterms. Then make three plots of vs. with Nterms=10, Nterms=100, and Nterms=400. This will give you an instructive way of seeing how the Fourier representation improves as you keep more terms. Because your procedure has more than one input variable, you will need to use a plot command of form plot(f(Nterms,theta),theta=0..Pi) .

--------------------------------------------------------------------------------

Problem 8.11

Write a procedure that takes as input the parameter , then solves the differential equation

with and . Let the procedure return as a single number (you will likely need to use trace to get this right). Then try different values of and see if you can find at least 2 values of (other than 0) for which . Find these 2 values of accurate to 5 significant figures.

This is a pretty hard problem, so here are some hints. You will want to put dsolve(...type=numeric) into the procedure and have it return a number. But the output from the numerical version of dsolve is a set of equations, like this

> restart;

> f:=dsolve({diff(y(t),t$2)=-y(t),y(0)=2,D(y)(0)=0},{y(t)},type=numeric);

> g:=f(1);

When we encountered this before we used the assign statement to get the number assigned to y(t), like this

> assign(g);

Error, (in assign/internal) invalid left hand side in assignment

This failed because the derivative symbol is not a valid variable to assign to. So now you have two choices. (i) Change the second order ode into a coupled set of ode's and use the assign statement to extract the number or (ii) learn how to use Maple's op and rhs commands to extract the number. The command op picks out elements of almost anything Maple generates: expressions, functions, lists, sets, etc.. Go back up to the restart above, execute the restart ,the dsolve , and the g:=f(1); but do not execute the assign statement. Then in the statement below change the 1 to 2 and then to 3, and see what you get back

> op(1,g);

Error, wrong number (or type) of parameters in function diff

The number f(1) is really close now; all we have to do is get it out of the equation, and that's just what rhs does: it extracts the right-hand side of an equation. So see what happens when you do this

> rhs(op(2,g));

Error, wrong number (or type) of parameters in function diff

You will need to use this construction inside your procedure so that it can return a number.

>

--------------------------------------------------------------------------------

Problem 8.12

Write a procedure that changes the Pnm(x) Associated Legendre Function discussed in Problem 3.4 into a real function that returns a number for any valid choice of n and m, including m=0. This means that you have to solve the "differentiate zero times" problem discussed in Chapter 3. Associated Legendre functions Test your procedure sufficiently that you know it is working correctly. (a) First write the procedure so that it returns an expression in x, then (b) change it so it returns a floating point number.

Go to top of section

Example: plotting results from fsolve

Debugging

Here is something that it is a little technical, but I run up against it all the time, so maybe you will too.

Suppose you have a hard equation to solve that has a variable parameter in it, like this

where varies from 0 up to 20. I want to build a function that returns the solution of this equation for any that I give it, and then I want to plot this function . The only way I know to do this is to write a procedure. After the Maple commands inside the procedure below have executed, F is assigned the last result calculated in the procedure, in this case the value of x from the fsolve , which is just what we want F to be.) Study the procedure below, then execute the plot.

> restart;

Declare F to be a procedure

> F:=proc(k)

Declare variables that will only be used inside the procedure

> local s,x;

Use fsolve to solve the equation using the value of k passed in through proc(k)

> s:=fsolve(cos(x)-k*x,x,0..2);

Return

> end;

Try a few values

> F(1);F(2);

Make a plot (use operator notation--the form plot(F(k),k=0..20) doesn't work)

> plot(F,0..20);

>

Problem 8.13

Consider the polynomial with a variable parameter. You can give this polynomial to fsolve and it will find all four roots as long as has a value.

> a:=2.5;

> P:=x^4+a*x^3+x^2+x+1;

> s:=fsolve(P,x,complex);

If you want to select one of the roots, say the second one, you just do this

> s[2];

Copy the procedure at the start of this section and modify it to plot the second root of this polynomial as varies between 2 and 10. (Just put the Maple code in this problem into the procedure.) Note: procedures return the last value calculated, so in your new procedure, right after the s:=fsolve(....) you will need the command s[2]; to return the value of the second root.

>

---------------------------------------------------------------------

Problem 8.14

(a) Write a procedure that returns the function defined as the solution of the implicit equation

.

In the procedure use Maple's fsolve command to do the solve and when you have built it, plot t(x) on the range x=0..10.

(b) Do the same thing as in (a), but this time use the function notation to define , i.e., t:=x->fsolve(...). Are the two ways of doing this equivalent?

(c) Write either a function or a procedure that will give you the derivative of this implicitly defined function, i.e., get a Maple function or procedure for . Warning: do not just take the -derivative of the equation above treating like a constant-- is really on both sides of the equation . So rewrite the equation with replaced by , then use diff to differentiate the whole equation with respect to . Then have Maple solve for and evaluate the resulting expression by modifying your procedure for . When you have built this procedure, use it to plot on the range 0..10 and verify visually that it looks like the derivative of the function you plotted in part (a).

--------------------------------------------------------------------------------

Go to top of chapter

>

>