Application Center - Maplesoft

App Preview:

Intro to programming in Maple - logic, loops, and procedures

You can switch back to the summary page by clicking here.

Learn about Maple
Download Application



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


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;

f := proc (x) options operator, arrow; if x < 0 the...
f := proc (x) options operator, arrow; if x < 0 the...
f := proc (x) options operator, arrow; if x < 0 the...
f := proc (x) options operator, arrow; if x < 0 the...

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;

g := proc (x) options operator, arrow; if x < 0 the...
g := proc (x) options operator, arrow; if x < 0 the...
g := proc (x) options operator, arrow; if x < 0 the...
g := proc (x) options operator, arrow; if x < 0 the...

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);

[Maple Plot]

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));

g := proc (x) options operator, arrow; piecewise(x ...

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

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

[Maple Plot]

integrate it like a function

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


and differentiate it like a function

> diff(g(x),x);

PIECEWISE([-2/(1+x^2)^2*x, x <= 0],[-sin(x), 0 < 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:


(2+t)^3/6 if t is in the range -2..-1


(1+(1+t)*(3+(1+t)*(3-3*(1+t))))/6 if t is in the range -1..0

B(t) = ______________________________________________________________

(1+(1-t)*(3+(1-t)*(3-3*(1-t))))/6 if t is in the range 0..1


- (t-2)^3/6 if t is in the range 1..2


For t outside of the range -2..2 , B(t) = 0 .

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 -2*Pi .. 2*Pi , 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


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;

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 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;

answer := 1

answer := 5

answer := 14

answer := 30

answer := 55

answer := 91

answer := 140

answer := 204

answer := 285

answer := 385

answer := 506

answer := 650

answer := 819

answer := 1015

answer := 1240

answer := 1496

answer := 1785

answer := 2109

answer := 2470

answer := 2870

answer := 3311

answer := 3795

answer := 4324

answer := 4900

answer := 5525

answer := 6201

answer := 6930

answer := 7714

answer := 8555

answer := 9455

answer := 10416

answer := 11440

answer := 12529

answer := 13685

answer := 14910

answer := 16206

answer := 17575

answer := 19019

answer := 20540

answer := 22140

answer := 23821

answer := 25585

answer := 27434

answer := 29370

answer := 31395

answer := 33511

answer := 35720

answer := 38024

answer := 40425

answer := 42925

answer := 45526

answer := 48230

answer := 51039

answer := 53955

answer := 56980

answer := 60116

answer := 63365

answer := 66729

answer := 70210

answer := 73810

answer := 77531

answer := 81375

answer := 85344

answer := 89440

answer := 93665

answer := 98021

answer := 102510

answer := 107134

answer := 111895

answer := 116795

answer := 121836

answer := 127020

answer := 132349

answer := 137825

answer := 143450

answer := 149226

answer := 155155

answer := 161239

answer := 167480

answer := 173880

answer := 180441

answer := 187165

answer := 194054

answer := 201110

answer := 208335

answer := 215731

answer := 223300

answer := 231044

answer := 238965

answer := 247065

answer := 255346

answer := 263810

answer := 272459

answer := 281295

answer := 290320

answer := 299536

answer := 308945

answer := 318549

answer := 328350

answer := 338350

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;

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;

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.;

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;

x := .2836621855

x := .9600369303

x := .5734897327

x := .8400126809

x := .6674533830

x := .7854005360

x := .7071051035

x := .7602456870

x := .7246667298

x := .7487203837

x := .7325605057

x := .7434644380

x := .7361281031

x := .7410737901

x := .7377440894

x := .7399878116

x := .7384767772

x := .7394947924

x := .7388091199

x := .7392710309

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 cos(x) = x . (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;

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;

tempx := .2836621855

xold := 5.

x := .2836621855

tempx := .9600369303

xold := .2836621855

x := .9600369303

tempx := .5734897327

xold := .9600369303

x := .5734897327

tempx := .8400126809

xold := .5734897327

x := .8400126809

tempx := .6674533830

xold := .8400126809

x := .6674533830

tempx := .7854005360

xold := .6674533830

x := .7854005360

tempx := .7071051035

xold := .7854005360

x := .7071051035

tempx := .7602456870

xold := .7071051035

x := .7602456870

tempx := .7246667298

xold := .7602456870

x := .7246667298

tempx := .7487203837

xold := .7246667298

x := .7487203837

tempx := .7325605057

xold := .7487203837

x := .7325605057

tempx := .7434644380

xold := .7325605057

x := .7434644380

tempx := .7361281031

xold := .7434644380

x := .7361281031

tempx := .7410737901

xold := .7361281031

x := .7410737901

tempx := .7377440894

xold := .7410737901

x := .7377440894

tempx := .7399878116

xold := .7377440894

x := .7399878116

tempx := .7384767772

xold := .7399878116

x := .7384767772

tempx := .7394947924

xold := .7384767772

x := .7394947924

tempx := .7388091199

xold := .7394947924

x := .7388091199

tempx := .7392710309

xold := .7388091199

x := .7392710309

tempx := .7389598975

xold := .7392710309

x := .7389598975

tempx := .7391694877

xold := .7389598975

x := .7391694877

tempx := .7390283084

xold := .7391694877

x := .7390283084

tempx := .7391234099

xold := .7390283084

x := .7391234099

tempx := .7390593490

xold := .7391234099

x := .7390593490

tempx := .7391025015

xold := .7390593490

x := .7391025015

tempx := .7390734336

xold := .7391025015

x := .7390734336

tempx := .7390930142

xold := .7390734336

x := .7390930142

tempx := .7390798245

xold := .7390930142

x := .7390798245

tempx := .7390887092

xold := .7390798245

x := .7390887092

tempx := .7390827244

xold := .7390887092

x := .7390827244

tempx := .7390867558

xold := .7390827244

x := .7390867558

tempx := .7390840402

xold := .7390867558

x := .7390840402

tempx := .7390858695

xold := .7390840402

x := .7390858695

tempx := .7390846372

xold := .7390858695

x := .7390846372

tempx := .7390854673

xold := .7390846372

x := .7390854673

tempx := .7390849082

xold := .7390854673

x := .7390849082

tempx := .7390852848

xold := .7390849082

x := .7390852848

tempx := .7390850311

xold := .7390852848

x := .7390850311

tempx := .7390852020

xold := .7390850311

x := .7390852020

tempx := .7390850869

xold := .7390852020

x := .7390850869

tempx := .7390851644

xold := .7390850869

x := .7390851644

tempx := .7390851122

xold := .7390851644

x := .7390851122

tempx := .7390851474

xold := .7390851122

x := .7390851474

tempx := .7390851237

xold := .7390851474

x := .7390851237

tempx := .7390851396

xold := .7390851237

x := .7390851396

tempx := .7390851289

xold := .7390851396

x := .7390851289

tempx := .7390851361

xold := .7390851289

x := .7390851361

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, cos(x) = sum((-1)^(n/2)*x^n/n!,n = 0 .. infinity) for n even . Give x 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 f(theta) = sum((-1)^n*sin((2*n+1)*theta)/(2*n+1),n ... . Give theta 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 1/(2*n+1) . (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 x = exp(-x) 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 f(x) = 0 , for x . The secant method starts by asking you for two reasonably close guesses x[1] and x[2] . Then you evaluate the function at each value of x to get f[1] = f(x[1]) and f[2] = f(x[2]) . You learned in high school algebra what to do with two ordered pairs ( x[1], f[1] ) and ( x[2], f[2] ) like this: use the two-point line formula to put a straight line through them: y-f[2] = (f[2]-f[1])*(x-x[1])/(x[2]-x[1]) . But what we have just built is a straight-line approximation to the function f(x) , so to get an approximation to the value of x that solves our hard equation, we just solve the line formula above for the value of x that gives y = 0 and call it x[3] , our next best guess at where the solution is.

> restart;

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

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

> x3:=solve(Eqline,x);

x3 := (f2*x1-x2*f1)/(f2-f1)


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:

x[3] = x[2]-f[2]*(x[2]-x[1])/(f[2]-f[1]) .

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 ( x[2]-x[1] ) and ( f[2]-f[1] ) are small, so we are dividing zero by zero again, but the values of f[1] and f[2] in the numerator are extra small because we are close to f(x) = 0 , 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 tan(x)-2*x = 0 for a root near x = 1 . (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 f[2] = f[1] 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


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;

f := proc (x, y) local a, b, z; global d; a := 3.; ...
f := proc (x, y) local a, b, z; global d; a := 3.; ...
f := proc (x, y) local a, b, z; global d; a := 3.; ...
f := proc (x, y) local a, b, z; global d; a := 3.; ...
f := proc (x, y) local a, b, z; global d; a := 3.; ...

> d:=31.;

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;

cosdeg := proc (theta) local phi, answer; phi := .5...
cosdeg := proc (theta) local phi, answer; phi := .5...
cosdeg := proc (theta) local phi, answer; phi := .5...
cosdeg := proc (theta) local phi, answer; phi := .5...

Now test it

> yy:=cosdeg(45.);

yy := cos(.2500000000*Pi)

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

phi := .2500000000*Pi

answer := cos(.2500000000*Pi)

<-- 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 theta and returns the value of the Fourier series f(theta) 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 f(theta) vs. theta 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 a , then solves the differential equation

diff(f(x),`$`(x,2)) = -f(x)^3 with f(0) = 0 and diff(f,x)(0) = a . Let the procedure return f(1) as a single number (you will likely need to use trace to get this right). Then try different values of a and see if you can find at least 2 values of a (other than 0) for which f(1) = 0 . Find these 2 values of a 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);

f := proc (rkf45_x) local i, comp_soln_data, odepro...

> g:=f(1);

g := [t = 1., y(t) = 1.0806046976628, diff(y(t),t) ...
g := [t = 1., y(t) = 1.0806046976628, diff(y(t),t) ...

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


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 k in it, like this

cos(x) = k*x

where k varies from 0 up to 20. I want to build a function F(k) that returns the solution x of this equation for any k that I give it, and then I want to plot this function F(k) . 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);


> end;

F := proc (k) local s, x; s := fsolve(cos(x)-k*x,x,...
F := proc (k) local s, x; s := fsolve(cos(x)-k*x,x,...

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);

[Maple Plot]


Problem 8.13

Consider the polynomial x^4+a*x^3+x^2+x+1 with a a variable parameter. You can give this polynomial to fsolve and it will find all four roots as long as a has a value.

> a:=2.5;

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

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

a := 2.5

P := x^4+2.5*x^3+x^2+x+1

s := -2.150706634, -.7857634353, .2182350347-.73763...
s := -2.150706634, -.7857634353, .2182350347-.73763...

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 a 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 t(x) defined as the solution of the implicit equation

t = x*exp(-t)/(1+t^2) .

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 t(x) , 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 diff(t(x),x) . Warning: do not just take the x -derivative of the equation above treating t like a constant-- t is really t(x) on both sides of the equation . So rewrite the equation with t replaced by t(x) , then use diff to differentiate the whole equation with respect to x . Then have Maple solve for diff(t(x),x) and evaluate the resulting expression by modifying your procedure for t(x) . When you have built this procedure, use it to plot diff(t(x),x) 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