Chapter8.mws
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: ifthenelseend 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
fthenelifelseend 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 nonnumber 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
ifthenelseend 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;
ifthenelifend 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 socalled cubic Bspline 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 threedimensional 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 1e8, 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 xxold>1e8 (this is why we need the initial fake value of xold)
>
for i from 1 while abs(xxold)>1e8 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 1e8. 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 1e4 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 twopoint line formula to put a straight line through them:
. But what we have just built is a straightline 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:=0f2=(f2f1)/(x2x1)*(xx2);
>
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 1e10 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 selfcontained 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 builtin Maple commands.
Problem 8.9
Write a procedure that computes the cubic Bspline 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 onedimensional function you would use
B:=proc(x)
, for 2d 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 righthand 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 notationthe 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
>
>