Simple equations such as 3x + 10 = 16 can be solved easily by rearranging the equation to the form x = 2. As we have seen with Maple there are many equations that can be solved this way, but there are also many that cannot be solved by any formula or rules. How do we find the solution to these equations?

Maple has the command fsolve and Matlab has a routine fzero that can be used to solve equations. In this section, we will examine how these methods work. From our discussion, you should get some appreciation of how these methods work. fsolve and fzero are not symbolic methods, but rather numerical methods for solving equations.

The Bisection Method, the Successive Substitution Method, Newton's Method, and the Secant Method are examples of these numerical methods. You may have heard of some of these methods in calculus or other math classes. Maple and Matlab use combinations of advanced numerical methods to solve equations. In this section, we will introduce you to the Bisection Method, and then show the Maple and Matlab routines that can be used to find roots and evaluate integrals.

Numerical methods find solutions by converging upon the answer. Convergence is the act of coming together. In this section, you will see how the process of convergence is used in root finding and evaluating integrals. In this lecture we discuss methods for solving equations in one variable. There are many engineering problems that require the solution of a nonlinear equation, that is equations that cannot be rearranged so the variable representing the unknown quantity is to the first power. Nonlinear equations may have multiple real roots, or no real roots at all. We will work with two specific problems to illustrate several numerical methods for root finding.

As an example of where nonlinear equations arise, consider a problem from structural engineering.

An important consideration in structural engineering is determining the maximum load that a structure can bear. In the structural analysis of the two-member frame shown in the figure, the critical load is given by the equation shown below. This equation is derived using methods from engineering statics. There are several assumptions about this frame that we will not dwell on.

This equation states that the critical load, P, is dependent on the length L, the stiffness EI, and a parameter z.

The desired value of z is the smallest positive root of this equation. We will refer to this equation for z as the two member frame problem.

Given the value of the parameter z, the length of the beams L, and the beam stiffness EI, then the critical load P can be determined. From a design standpoint, an engineer would want to make sure that the loads on a structure are much less than the critical load so that the structure does not buckle.

Notice that there is no way to manipulate the equation for z to obtain a nice formula to evaluate z. This is a nonlinear equation.

As another example, consider a cantilevered beam of length L that is fixed at the left side to be horizontal, as shown in the figure. If the beam has weight per unit length of w and stiffness EI, then the equation for the position of the center of the beam is that shown below.

Or, if we let a = x/L, so that α ranges from 0 to 1, then the equation becomes:

This equation gives the height y as a function of the fraction of the length to the end.

Notice that the maximum deflection is at the end of the beam, at α=1, and the deflection there is - wL^4/(8 EI).

We want to find the value of α where the beam has half the total deflection. We need to solve the equation where the function y(α) has the value one-half the maximum deflection. This is the equation shown here.

The solution of this equation gives the value of α so that at &alphalL the beam has half the total deflection. We will refer to this equation as the cantilevered beam problem.

These two engineering problems have equations that cannot be manipulated symbolically to get a formula for the solution. We must use numerical methods to solve for the root of these equations.

The bisection method is an example of an iterative numerical method for root finding. That is, it generates sequences of values that, hopefully, converge to the root of the equation.

The bisection method works with the standard equation form of f(x) = 0. The equations from the two member frame problem and the cantilevered beam problem can be rearranged to satisfy this form as follows.

The bisection method begins with an interval [l,r] with f(l) and f(r) having opposite signs, see figure below. If one of f(l) and f(r) is positive and the other is negative, and f is a continuous function, then somewhere between l and r the function f must take on the value of zero. This is the point we wish to find. The variables l and r are used to denote the left and right endpoints of the interval.

At each step of the bisection method, the interval is replaced by an interval with half the width. Each step consists of the following calculations:

Use the bisection method to find the root(s) of the equation for the buckling load of the two-member frame. The equation is rearranged to get it in the standard form

From a graph of the function we see that the desired root is between 2 and 4, and closer to 4. The first steps of this method are illustrated in the table below.

step	l	m	r	f(l)	f(m)	f(r)
0	2	3	4	-3.18504	-1.06562	0.357821
1	3	3.5	4	-1.06562	-0.48695	0.357821
2	3.5	3.75	4	-0.48695	-0.1339	0.357821

The first values of l and r are 2 and 4, respectively. Note that in this case f(l) is negative and f(r) is positive. The value of m is 3 and because f(m) is negative at step 0, the new interval is [3,4]. That is l = 3 and r remains at 4.

With the interval [l,r] being [3,4] the new midpoint is 3.5. The value of f(m) is now negative, so the next interval is [3.5,4].

The computation for 12 iterations is shown in the following table. The values of l, m, and r are converging to a value near 3.8288, while the values of f(l), f(m), and f(r) are approaching zero.

step	l	m	r	f(l)	f(m)	f(r)
0	2	3	4	-3.18504	-1.06562	0.357821
1	3	3.5	4	-1.06562	-0.48695	0.357821
2	3.5	3.75	4	-0.48695	-0.1339	0.357821
3	3.75	3.875	4	-0.1339	0.085953	0.357821
4	3.75	3.8125	3.875	-0.1339	-0.02903	0.085953
5	3.8125	3.84375	3.875	-0.02903	0.027049	0.085953
6	3.8125	3.828125	3.84375	-0.02903	-0.00132	0.027049
7	3.828125	3.835938	3.84375	-0.00132	0.012778	0.027049
8	3.828125	3.832031	3.835938	-0.00132	0.005706	0.012778
9	3.828125	3.830078	3.832031	-0.00132	0.002187	0.005706
10	3.828125	3.829102	3.830078	-0.00132	0.000431	0.002187
11	3.828125	3.828613	3.829102	-0.00132	-0.00045	0.000431
12	3.828613	3.828857	3.829102	-0.00045	-8E-06	0.000431

If we continue the computation by adding more iterations, the value of z converges to 3.8288.

The bisection method is not a very fast method. It takes many evaluations of the function f to get a good solution. Notice that the bisection method only uses the sign of the function to make its decision about the next interval. Methods that use information about the magnitude of f(x) can be much faster than the bisection method. Examples of such methods are Newton's method and the secant method.

Maple has a numerical routine, the fsolve command, for finding roots. The main difference between the solve and fsolve commands are that solve is for symbolic computations and fsolve is for numerical computations. 'f' in Maple commands stands for floating points aka decimals. So whenever you compute using fsolve or evalf, your answer will be returned in decimals. The routines fsolve in Maple and fzero in Matlab may use the bisection method as part of their solution strategy, but they also use methods that are much more efficient.

All of these methods are iterative. They have a way to improve on an estimate to the solution. The methods terminate when the solution has been determined to sufficient accuracy.

Both Maple and Matlab have routines that will perform the successive iterations necessary to solve non-linear equations.

In Maple, roots can be found by using the fsolve command, which is a bisection-like method. This command may have been introduced in the Symbolic Computation section. For a quick review, let's say you wanted to find the roots of the following equation:

With the fsolve command, you must put in an interval range. fsolve does use some symbolic information about the function. If the function is a polynomial, then fsolve will solve for all the roots in the specified interval. If the function is not a polynomial, then fsolve will return with only one root. If there is more than one root in the interval, it is not clear which root Maple will return with. For this reason, it is always best to use plotting to determine a good interval that contains the one root being sought. Remember that fsolve only finds the real roots of an equation.

Important Help: Use plotting before using fsolve. This way, you can approximate the roots and enter smaller ranges for fsolve to check over.

Remember, that the smallest positive value of z gives you the maximum critical buckling load. So we will first plot to see what range we should solve our root for.

Using the plot commands, scaling=constrained and discont=true, allows you to have a better view of the plot.

Note that there are many solutions to f(z)=0. By looking at the graph you can choose your range for fsolve. We can't tell if a root is on zero or a little bit greater than zero so we'll first try

and Maple will return the answer of 0, which is a solution to the equation, but is not positive so it is not be the solution to the problem that we want. Now we'll try

and Maple will now return 3.828861865, which is the value to be used in computing the critical load.

fsolve does not require that an interval be specified, but without specifying the interval the user has no control over which root Maple will find.

The Matlab routine fzero can also be used to find roots of equations. To use fzero you should define the function in an M-file. Let's find the roots of the following equation:

The general form for using the fzero routine is x = fzero('function_name', x_value_guess or x range). To find our approximate x values, let's plot the function first. Create an M-file called "f1plot.m" and enter the following:

Save and run f1plot.m. From the graph, you can see that there are two roots of the equation. They are near x = -1 and x = 3. To find the exact values, use the fzero routine:

x1 = fzero('f1',-1) or x1 = fzero('f1', [-1,0]) x2 = fzero('f1',3) or x2 = fzero('f1', [2, 3])

Now we'll use Matlab to solve for the cantilevered beam problem. We need to find the root of

that lies between 0 and 1. (Why are we not interested in roots outside of the interval between 0 and 1? ) First, we'll define the function:

Note that Maple can obtain a solution to the cantilevered beam problem using symbolic methods. However, since we are interested in the numerical value between 0 and 1, it is as quick and convenient to use fsolve as it is to use solve followed by evalf.

So what's the difference between fsolve and fzero?

In many cases, fsolve and fzero can be used interchangeably (within their respective programs, i.e. Maple or Matlab). They will both find real roots of an equation. For both routines it is important to have some estimate of the root. In both Maple and Matlab, it is best to see a plot of the function before you execute the fsolve or fzero command. fsolve will find all the real roots of a polynomial. Matlab has other commands to find all the roots of a polynomial.

Another method for solving equations is the method of successive substitution. This is also called the fixed point method. It is especially useful since it can be quickly programmed on hand calculators and spreadsheets.

The process of successive substitution starts with an initial approximate value of x0 and generates additional values x1, x2, x3, ... using the formula: xk+1 = g(xk).

Starting with α = 0.5, we compute successive values for α: 0.705882, 0.578272, 0.645047, 0.606224, etc. After about 25 iterations the values converge to 0.619775

Other formulas for the cantilevered beam problem could also be used. The formula

will converge even faster, in about 10 iterations. Starting with α = 0.5, we obtain 0.594089, 0.614173, 0.618549, 0.619506, etc. converging to the same value, of course. Not all formulas will give convergent methods, some diverge, that is, get farther away from the root. With some problems you may have to try several formulas before you find one that converges to the root you want.

Successive substitution is a good method to use on simple equations that can be easily programmed on calculators, spreadsheets, and similar computational aids, but it is unreliable for general purpose computations.

In the section on symbolic computation, you saw several examples of how Maple can evaluate integrals symbolically. You also saw how Maple uses special functions to obtain expressions for some integrals. However, there are many integrals that can not be solved by any of the rules known to Maple. In addition, there are times when we do not need the symbolic formula for an integral, we only need the numerical value. In these cases, we can obtain values for integrals using numerical integration.

Before discussing numerical methods for integration, we present several examples that we will use throughout our discussion of numerical integration.

You are working with a tank that has a radius that varies by height according to the following relationship (all measurements in feet):

From calculus, we know that we can find the volume of a container with varying radius by using the following expression:

Although this is an integral that we could actually compute using paper and pencil, we probably don't want to evaluate it that way. This is a complex problem. We will return to it later and show how to evaluate it using both Maple and Matlab.

In a factory it has been found that the failure of a critical electric motor has the probability distribution

where m is the number of months to failure. What this means is that the probability that a motor (installed at time 0) will fail between months a and b is the integral of f(m) between a and b. That is,

We do not know the distribution beyond 12 months because these motors are always replaced after a year of use.

First, we ask "What is the probability that a motor survives a year of use?" To find this we compute the probability that a motor will fail. This is

For integrals such as P₁₂ Maple can evaluate it, but the answer is not too informative. All we really need is the numerical value. Using evalf and Int on this integral will produce the value of 0.736120123 Maple determined this value by approximating the integral using a numerical, rather than a symbolic, method.

Matlab does not work symbolically, it works numerically, and it has a procedure, called quad, that can be used to evaluate integrals such as from the volume problem or the motor problem.

We now look at a method used to numerically evaluate integrals, called the trapezoid method. By understanding how the trapezoid method works, you can get an appreciation for how other numerical methods for evaluating integrals work.

We divide the interval [A,B] into N subintervals, each of length h. For notational convenience we will use [A:B] to denote an interval that is divided into several subintervals of length h and we will use [a:b] to denote an interval of length h.

The idea of the trapezoid rule is displayed in Figure 3.1. The area under the graph of f(x) is approximated by the areas of the trapezoids formed by connecting the points (x_j, f(x_j)) by straight-line segments.

Recall that a trapezoid is a quadrilateral with one pair of opposite sides being parallel. The area of a trapezoid is equal to the average of the lengths of the two parallel sides times the distance between the lines. In Figure 3.1, the areas above each length of size h are trapezoids with the vertical sides being parallel.

Notice that the values of the function at the endpoints of each increment are summed with the first and last points of the entire interval having a weight of one-half.

As the value of h is made smaller, the trapezoid formula gives better approximations to the value of the integral. By choosing a small enough value of h, we can obtain the value of the integral to a high degree of accuracy.

Both Maple and Matlab have methods for evaluating integrals numerically. The procedures are similar to the trapezoid method. However, they are more sophisticated in that they use different sizes of subintervals depending on the amount of variability of the function. In this way they can adjust the size of the subinterval so that significant accuracy can be obtained in an efficient way. Where the function has great variation, the values of h are smaller, where the function changes less rapidly, larger values of h are used.

In Maple using the command evalf on an integral invokes the numerical evaluation of the integral. If you know that Maple cannot evaluate the integral or you are not interested in the analytical expression for the integral, you can use the Int command, so that Maple does not try to evaluate the integral, and then follow this with the evalf command.

7.1 Matlab

Matlab does not use symbolic methods so the only way to evaluate integrals in Matlab is numerically. The routine to evaluate integrals in Matlab is called quad (from quadrature, which is another name for numerical integration). The command q = quad(fun,a,b) approximates the integral of function fun from a to b using a method called recursive adaptive Simpson quadrature.

To evaluate the integral for the tank volume we must first define the function in an M-file, tank.m.

Notice the periods, this is because Matlab deals with everything as a matrix and when we place the periods there, Matlab will know to multiple element-wise rather than executing a matrix multiplication.

Matlab will return with ans = 2.9652e+003 In other words, our tank volume is 2965.2 cubic feet. Remember, if you plan to use this value again, name it (such as TVol = quad('tank',0,5). That way, TVol = 2965.2)

Note: Although it is good programming practice to define the function in a separate M-file, the quad routine can be defined with a single line:
TVol = quad('pi*(h.^2+(1/3)*h.^sqrt(3)).^2',0,5)

We now consider evaluating the integral from the motor problem using the quad command.

Matlab should return the answer Q = 0.7361 So 73.61% of the motors fail within a year. Therefore, 26.39% last for a year.

To continue the motor example, suppose that motors that are replaced are rebuilt. We ask "What is the average age of motors that are sent to the shop to be rebuilt?"

The average age of the motors that fail between months a and b is given by the integral

This is an example of an 'expected value,' in this case it is the expected value of m.

This formula combines both the average age of those that fail and those that survive for a year. (Recall that motors are always replaced after a year of use.) The final value is 6.8250 months.

Matlab has a routine, dblquad, that can be used to evaluate two-dimensional integrals. To numerically evaluate double integrals, you can use the dblquad command, which works just like the quad command. This command in Matlab q = dblquad(fun,xmin,xmax,ymin,ymax) calls the quad function to evaluate the double integral fun(x,y) over the rectangle xmin <= x <= xmax, ymin <= y <= ymax. fun(x,y) must accept a vector x and a scalar y and return a vector of values of the integrand.

Let's say you have a sphere with the center at the origin given by the equation:

Now you decided to pull a 4 × 4 square plug, also centered at the origin, out of the sphere as shown below.

How would you calculate the volume of that square plug? Since z is the height of your sphere and varies as we move along the plug, we will now set the equation equal to z.

To solve for the volume, we will have to take the double integral of the height since that is what is varying.

Notice that you need twice the value of z because the plug extends equal values above and below the origin (from �z to +z). Now if we substitute in the equation and the range, we get this equation:

This integral is rather difficult to do by hand so let's practice some more with Matlab. We'll start off by creating an M-file called "plug.m"

To numerically evaluate double integrals, you can use the dblquad command, which works just like the quad command. q = dblquad(fun,xmin,xmax,ymin,ymax) calls the quad function to evaluate the double integral fun(x,y) over the rectangle xmin <= x <= xmax, ymin <= y <= ymax. fun(x,y) must accept a vector x and a scalar y and return a vector of values of the integrand.

You can also do this problem using Maple. You will have to do it using the int command twice. To do this you would type this command into Maple

Maple will compute a symbolic answer. Remember our discussion of floating numbers; we can use the command evalf to get a numerical solution. So our next command will be

Notes

Convergence

CS 310 and NEEP 271

1 Introduction

2 The Two Member Frame

3 The Cantilevered Beam

4 Solving for Roots of Nonlinear Equations

4.1 The Bisection Method

Example 1:

4.2 Root Finding with Mathematical Tools

4.2.1 Maple

4.2.2 Matlab

4.3 Successive Substitution

Example 2:

5 Numerical Evaluation of Integrals

Example 3: The Volume of a Tank.

Example 4: Electric Motor Failure

6 Trapezoid Method

7 Evaluating Integrals with Mathematical Tools

7.1 Matlab

Example 5:

7.2 Two Dimensional Integrals

Example 6: