** Next:** About this document ...
**Up:** No Title
** Previous:** No Title

Solving Equations

Numerical techniques play a key role in computations, and Maple offers a variety of means for their performing. The present lab is composed in order to concentrate on the computational tools rather than on mathematical issues. Trying to make your knowledge of Maple more systematic, the lab pursues the following goals:

- (i)
- to familiarize you with various Maple tools for numerical solving of equations and evaluating of definite integrals;
- (ii)
- to provide you with basic experience in working with Maple as programming language.

This lab consists of two major parts - Background and Exercises. In the first, the relevant theoretical notes and description of the use of appropriate Maple commands and procedures are given. The second part includes four problems devoted to different numerical methods and various tools of Maple.

We know simple formulas for solving linear and quadratic equations, and there are more complicated formulas for cubic and quartic equations (equations of degree three and four). At one time it was hoped that similar formulas might be found for quintic and higher degree equations, but Niels Abel (a great Norwegian mathematician of the beginning of 19th century) showed that no formulas like these are possible for polynomial equations of degree greater than 4.

When exact formulas for solving an equation *f*(*x*) = 0 are not available, we
can turn to numerical techniques from Calculus to approximate the solutions we
seek. Two techniques considered in the course (the Bisection Method and
Newton's method) are the examples of an idea of successive approximations
allowing to localize the point (i.e., solutions, or roots) where the graph of
*y* = *f*(*x*) passes through zero.

The idea of the Bisection Method is simple. If we know that a zero exists in
the closed interval [*a*, *b*], the zero must lie in the sub-interval
[*a*, (*a*+*b*)/2], or [(*a*+*b*)/2, *b*]. From the sign of *f*[(*a*+*b*)/2], we can determine which interval contains the
zero. By repeatedly bisecting the interval, we can ``close in" on the zero of
the function.

Newton's Method is based on the idea of using tangent lines to replace the
graph of *y* = *f*(*x*) near the points where *f* is zero. The initial estimate *x _{0}* may
be found by graphing or just plain guessing. The method then uses the tangent
to the curve

Convergence of the Bisection Method is an almost obvious issue. The required
accuracy can be initially chosen, and the iterations stop at step *n* when the
error becomes less than (or equal to) the interval's half-length, i.e.:

where *a*_{n} and *b*_{n} are the left and right points of the
interval at the nth step respectively.

For Newton's approach, the convergence is not so obvious; moreover, the method
does not always lead to convergence and requires an appropriate analysis (see
Theorem A in Section 10.4 of the text). In practice, Newton's Method usually
converges with impressive speed, but since this is not guaranteed one must
test that convergence is actually taking place. The relevant tool for this
testing is a sufficient condition of convergence establishing that the method
converges to zero of *f* if

on an open interval containing the root.

Another important application of numerical techniques is integration. We may be unable to evaluate an integral (because of the absence of an antiderivative in terms of elementary functions) or desire to compose a procedure for a computer or calculator to evaluate definite integrals. Whatever the reason, we turn to numerical methods such as the Trapezoidal Rule and Simpson's Rule.

These methods actually estimate the area under a curve *y* = *f*(*x*) on the
interval [*a*,*b*]. The accuracy of the entire procedure depends on how well the
upper boundary of each approximating sub-region fits the shape of the given
curve. It is clear, therefore, that trapezoids often result in a better
approximation than rectangles (recall the geometric interpretation of
the Riemann Sum), and Simpson's Rule using parabolic arcs for its
upper boundaries provides a better approximation than trapezoids.

To approximate with the use of the
Trapezoidal Rule, we apply the formula:

There is also an error term that can be used to estimate the error. More
precisely, we have:

where

for some point *c* between *a* and *b*, provided exists on [*a*,*b*].

To approximate the definite integral by Simpson's Rule, the following formula
is used:

An error control in Simpson's Rule term is established through the analysis of
the term *E*_{n} in (2) that now is specified as

where, again, .

There is a number of ways to implement the numerical techniques noted above working with Maple. In the present lab, you have to check out some of them; particularly, the special Maple commands directly implementing these methods and the procedures composed with the element of Maple programming means.

Maple is a very powerful computation system, so it is nothing surprising that many popular numerical procedures can be run just by a single command. In fact, some tools for solving equations were already in use in your previous labs. Here is a brief generalized review of these tools.

Maple has two commands for solving equations, `solve` and
`fsolve`, providing exact
and approximate solutions respectively. For example, to solve the quadratic
equation

we may either use

> sol1 := solve(x^2+2*x-1=0,x);

or label the left side of the equation by:

> f := x^2=2*x-1:

> sol2 := solve(f=0,x);

or assign a variable name to the entire equation entering:

> eq1 := x^2+2*x-1=0:

> solve(eq1,x);

More than one equation can be solved by using curly brackets { }. For example:

> eq_a := 3*x+2*y=1;

> eq_b := x+2*y=3;

> solve({eq_a,eq_b}, {x,y});

However, many equations cannot be solved exactly. If you forget Abel's Theorem and try finding the exact solutions to any polynomials of 5th degree or higher, you will encounter something like that:

> eq2 := x^7+3*x^4+2*x-1=0:

> solve(eq2,x);

Maple's response indicates that it does not know how to solve this equation
exactly. In this situation, use `fsolve` command:

> fsolve(eq2,x);

It is wise to use the **fsolve** command in conjunction with a plot in order to
catch other solutions which may exist. Consider an equation:

> eq3 := x^2+1/x-1/x^2=0:

> fsolve(eq,x);

However, the command:

> plot(x^2+1/x-1/x^2, x=-3..3, y=-20..20);shows that there is another branch of the graph of this function, and it crosses the

> fsolve(eq,x,x=0..1);

Generally, **fsolve** find all the roots, however, some more complicated equations
should be handled by a **plot** together with specifying a range with the
**fsolve**
command.

Maple can be used to plot an equation, but in this case the **implicitplot**
command must be used instead of plot. In order to use this command, first load
the **plots** package with the command

> with(plots);

Try the **implicitplot** command to plot a circle using an equation of a circle.

In order to conduct computation using Eqs. (1) and (3), it might be useful to enter the general equation of the function and then find its value in various points by substitution (command subs).

It is clear from the conceptual description of the methods above (see the text, Sections 10.3 and 10.4), that they can provide reasonable closeness to
the exact value if the number of approximating strips is large enough.
Obviously, such processes as partition and successive computations in the
Trapezoidal and Simpson's Method and iterations in the Bisection and Newton's
Methods could be algorithmized and ordered to a computer. Actually,
the **solve**
and **fsolve** Maple commands are nothing but the built-in procedures involving
sets of simple operations/commands and implementing such algorithms.

Maple possesses appropriate means which allow one to repeat some particular commands or set of commands many times with slightly different values; a number of commands could be gathered in an automating procedure. In this capacity, Maple reveals its features common for programming languages. The relative simplicity of the algorithms related to the numerical methods studying in this lab suggests the idea to introduce you to the elements of programming in Maple on the example of these methods.

To repeat a command you can compose a looping structure using the `
for` and `do`
commands. For example, the structure

> for k from 1 by 0.5 to 3

do print (k, k^2); od;prints out 5 pairs of numbers. Notice that the index moves from 1 to 3 by step 0.5. You may drop a portion

To type this structure, press **Shift + Enter** in order to move the cursor to next
line without executing the line. When you press enter at the end of the last
line, all the lines will be performed, causing 5 (or 3) pairs of numbers to be
displayed.

You can use this looping structure as part of a procedure:

> maketable :=

proc(n) for k from 1 by 0.5 to n do print (k, k^2); od; end;

Remember to press **Shift + Enter** at the end of every line but the last as you
enter this procedure. The procedure **maketable** uses *n* as its argument. Press
**Enter** to define this procedure. Only a compact version of the procedure
definition is displayed. Do not be concerned if a warning appears on the
screen.

To use this procedure later, just type:

> maketable(7);

As you might expect, a table of squares of 7 numbers is displayed.

The loop is one of the two most basic constructs in programming. The other
one is the **if** or conditional statement. It arises in many contexts. For
example, you can use the **if** statement to implement an absolute value
function:

Maple executes the **if** statement as follows: if *x* < 0, then Maple calculates
-*x*; otherwise, it calculates *x*.

The procedure **Abs** to compute the absolute value function is shown below. The
closing word **fi** (**if** in reverse) completes the **if** statement. In either case,
the absolute value of *x* is the last result that Maple computes and so is the
value that **Abs** returns.

> Abs :=

proc(x) if x < 0 then -x; else x; fi; end;

> Abs(3); Abs(-2.3);

3 2.3

Notice that the portion **else x** can be dropped.

Some specific commands usually appear in connection with procedures. Command
**local** specifies variables whose values Maple knows only inside the procedure;
these are called local variables. While Maple executes a procedure, a global
variable by the same name remains unchanged. Command **type** **(a, *)** checks the
type of a variable *a*; some key word is used as *****; for example,
**integer, even,
equation,** etc.

The syntax of all these commands is quite clear; additional information is available through Maple Help.

Find the equation of the parabola that passes through the points (-1,2), (1,-1.5), and (4,7).

Here is a Maple source code for a non-standard procedure implementing Newton's
Method.

**Note:** You may avoid typing the procedure by using the Getting Started
Worksheet attached tothis lab. The text for the Newton Procedure is
there.

> Newton := proc (f, a, b, iter)

local first, second, x, i, d1; if not type (f, procedure) then ERROR ('the first argument must be a function name.') fi; if nargs <> 4 then ERROR ('You need four arguments.') fi; if evalf(f(a)*f(b)) >= 0 then print ('f(a) and f(b) do not have opposite signs.'); fi; first := evalf((a+b)/2,15); print ('Initial estimate', evalf(first)); d1 := diff(f(x), x); second := first - evalf(f(first)/subs(x=first, d1), 15); for i from 2 to iter do first := second; print (i - 1, evalf(first(f(first)/subs(x=first, d1), 15); od; print ('Final estimate', evalf(second)); end;

- (a)
- Analyze the code. Explain how the iteration scheme of Newton's Method is
implemented in it. Express plausible hypotheses about the nature of the
particular cases handled by
**if**statements with text outputs. - (b)
- Consider function and find its root(s) using the
commands
**solve**and**fsolve**. Is the answer the same? Why? - (c)
- For the function in (b), apply the Newton procedure. Compare the result with (b). How many iterations did you need to get the root with nine digits of accuracy?
- (d)
- Repeat points (b) and (c) for the function
*f*(*x*) =*x*- 3^{3}*x*+ 1. Illustrate your solution by plotting the function.

Find the area under the curve *f*(*x*) =
*e*^{-x2} between 0 and 2 using the
formula of the Trapezoidal Rule (1) with *n* = 4 and 8. Evaluate the
corresponding integral by Maple tools. Compare the results and specify the
number of decimal places that are correct in the approximate
results.

Here is a Maple source code for a non-standard procedure implementing
Simpson's Method.

**Note:** You may avoid typing the procedure by using the Getting Started
Worksheet attached tothis lab. The text for the Newton Procedure is
there.

> Simpson := proc(f, a, b, n)

local i, j, h, s, v; if not type(f, procedure) then ERROR ('the first argument must be a function name.') fi; if nargs <> 4 then ERROR ('You need four arguments.') fi; if (n mod 2 <> 0) then ERROR ('# of subdivisions must be even.') fi; h := (b-a)/n; for i from 1 to (n-1) do v[i] := a + i*h; od; s := f(a) + f(b); for i from 1 by 2 to (n-1) do s := s + 4*f(v[i]); od; for j from 2 by 2 to (n-2) do s := s + 2*f(v[j]); od; s := (h/3) * s; evalf(s, 15); end;

[(a)] Analyze the code. Explain how the algorithm of the Simpson's Method is
implemented in it. Express plausible hypotheses about the nature of the
particular cases handled by **if** statements with text outputs.

[(b)] We know that

Use the Simpson procedure with 4 different *n* (the highest value *n* = 40) to get
4 approximations for . Use your watch to estimate time required to
perform each computation. Compare the results with the value of provided
by Maple.

2/2/1999