Then replace the GUID above with the vault GUID, the '0' with the object type ID, the '438' with the object id and the '476' with the file id. Check the documentation for more details. Given the parameter values V = 1 × 106m 3, Q = 1 × 105 m 3 /yr, W = 1 × 106 g/yr, and k = 0.25 m 0.5 /g 0.5 /yr, use the modified secant method to solve for the steady-state concentration. Employ an initial guess of c = 4 g/m 3 and δ = 0.5.

This code calculates roots of continuous functions within a given interval and uses the Bisection method. The program assumes that the provided points produce a change of sign on the function under study. If a change of sign is found, then the root is calculated using the Bisection algorithm (also known as the Half-interval Search). If there is no change of sign, an error is displayed. You could try with other low or high values, or you could improve the code to find two values with a different sign before going on.

In general, when we work with numerical methods we must be aware that errors may result for a number of reasons. First, a root may be calculated when it should not be. It could happen if a point is so close to zero that a root is found due to round-off error. Second, two roots may be so close together that the program never finds the opposite signs between them...

You will provide the function, the interval (both low and high values) and a tolerance. This suggested version of the method will list the evaluated intervals and the final approximation found using your tolerance.

functionm = bisection(f, low, high, tol)
disp('Bisection Method');

% Evaluate both ends of the interval
y1 = feval(f, low);
y2 = feval(f, high);
i = 0;

% Display error and finish if signs are not different
if y1 * y2 > 0
disp('Have not found a change in sign. Will not continue...');
m = 'Error'
return
end

% Work with the limits modifying them until you find
% a function close enough to zero.
disp('Iterlowhighx0');
while (abs(high - low) >= tol)
i = i + 1;
% Find a new value to be tested as a root
m = (high + low)/2;
y3 = feval(f, m);
if y3 0
fprintf('Root at x = %f nn', m);
return
end
fprintf('%2i t %f t %f t %f n', i-1, low, high, m);

% Update the limits
if y1 * y3 > 0
low = m;
y1 = y3;
else
high = m;
end
end

% Show the last approximation considering the tolerance
w = feval(f, m);
fprintf('n x = %f produces f(x) = %f n %i iterationsn', m, y3, i-1);
fprintf(' Approximation with tolerance = %f n', tol);

Let’s say that we want to find a root of the function
y = 5x⁴ - 2.7x² - 2x + 0.5

and want to explore the interval [0.1, 0.5], we could call the function from another m-file, like this:

my_fun = @(x) 5*x^4 - 2.7*x^2 - 2*x + .5;
low = .1;
high = 0.5;
tolerance = .00001;
x = bisection(my_fun, low, high, tolerance);

the result is:

Bisection Method
Iter low high x0
0 0.100000 0.500000 0.300000
Root at x = 0.200000

If we plot the function, we get a visual way of finding roots. In this case, this is the function

Now, another example and let’s say that we want to find the root of another function
y = 2.5x² - 3x + 0.5

using another interval, like [0, 0.5], we can use this code to call the half-interval search

my_fun = @(x) 2.5*x^2 - 3*x + .5;
low = 0;
high = 0.5;
tolerance = .00001;
x = bisection(my_fun, low, high, tolerance);

and we get this information from the proposed code above:

Bisection Method
Iterlowhighx0
0 0.000000 0.500000 0.250000
1 0.000000 0.250000 0.125000
2 0.125000 0.250000 0.187500
3 0.187500 0.250000 0.218750
4 0.187500 0.218750 0.203125
5 0.187500 0.203125 0.195313
6 0.195313 0.203125 0.199219
7 0.199219 0.203125 0.201172
8 0.199219 0.201172 0.200195
9 0.199219 0.200195 0.199707
10 0.199707 0.200195 0.199951
11 0.199951 0.200195 0.200073
12 0.199951 0.200073 0.200012
13 0.199951 0.200012 0.199982
14 0.199982 0.200012 0.199997
15 0.199997 0.200012 0.200005

x = 0.200005 produces f(x) = -0.000009
15 iterations
Approximation with tolerance = 0.000010

Again, plotting the function is a good idea to better know what we’re doing...

Maybe you're also interested in:

Other ideas about root finding

From 'Bisection Method' to home
From 'Bisection Method' to Generic Programming with Matlab

Top

Introduction	Exercise 1
A Sample Problem	Exercise 2
The Bisection Idea	Exercise 3
Programming Bisection	Exercise 4
Variable Function Names	Exercise 5
Exercise 6
Convergence criteria	Exercise 7
The secant method	Exercise 8
Regula Falsi	Exercise 9
Muller's method (Extra)	Extra Credit

The root finding problem is the task of finding one or morevalues for a variable so that an equation

is satisfied. Denote the desired solution as .As a computational problem, we are interested ineffective and efficient ways of finding a value that is ``closeenough' to . There are two common ways to measuring ``closeenough.'

1. The ``residual error,' , is small, or
2. The ``approximation error' or ``true error,' is small.

Of course, the true error is not known, so it must be approximatedin some way. In this lab we will look at some nonlinear equationsand some simple methods of finding an approximate solution. Wewill use an estimate of the true error when it is readily availableand will use the residual error the rest of the time.

You will see in this lab that error estimates can be misleading and thereare many ways to adjust and improve them. We will not examine moresophisticated error estimation in this lab.

This lab will take two sessions.If you print this lab, you may find the pdf versionmore appropriate.

You can find Matlab code on the internet and in books (Quarteroni,Sacco and Saleri's book is an example).This code appears in a variety of styles. The coding style used in examples in these labs has an underlying consistency and is, in myopinion, one of the easiest to read, understand and debug.Some of the style rules followed in these labs includes:

• One statement per line, with minor exceptions.
• Significant variables are named with long names starting withnouns and followed by modifiers, and with beginnings of modifierscapitalized.
• Loop counters and other insignificant variables are short, oftena single letter such as k, m or n. iand j are avoided as variable names.
• The interior blocks of loops and if-tests are indented.
• Function files always have comments following the functionstatement so that they are available using the help command.Function usage is indicated by repeating the signature line amongthe comments.

These rules of style improve clarity of the code so that people andnot just computers can make sense of it.

Clarity is important for debugging because if code is harder to read thenit is harder to debug. Debugging is one of the most difficult and least rewardingactivities you will engage in, so anything that simplifies debugging pays off.Clarity is also important when others read your code. Since I need to readyour code in order to give you a grade, please format your code similarlyto that in the labs.

Suppose we want to know if there is a solution to

whether the solution is unique, and what its value is. This is notan algebraic equation, and there is little hope of forming anexplicit expression for the solution.

Since we can't solve this equation exactly, it's worth knowing whetherthere is anything we can say about it. In fact, if there is a solution to the equation, we can probably find a computerrepresentable number x which approximately satisfies theequation (has low residual error) and which is close to (has alow approximation error). Each of the following two exercises uses aplot to illustrate the existance of a solution.

Note: Matlab has the capability of plotting somefunctions by name. We won't be using this capability in thiscourse because it is too specialized. Instead, we will constructplots the way you probably learned when you first learned aboutplotting: by constructing a pair of vectors of -values (abscissæ) andcorresponding -values (ordinates) and then plotting the pointswith lines connecting them.

Exercise 1: The following steps show how to use Matlab to plot the functions and together in the same graph from to .

1. Define a Matlab (vector) variable xplot to take on 100evenly-spaced values between -pi and pi. You can uselinspace to do this.
2. Define the (vector) variable y1plot by y1plot=cos(xplot).This is for the curve .
3. Define the (vector) variable y2plot by y2plot=xplot.This is for the line .
4. Plot both y1plot and y2plot by plotting the first,then hold on, plot the second, and hold off. You saw this donebefore in Lab 2.Note: If you read the help files for the ``plot' function,you will find that an alternative single command isplot(xplot,y1plot,xplot,y2plot).
5. The two lines intersect, clearly showing that a solution to the equation exists. Read an approximate value of at the intersection pointfrom the plot.
6. Include the Matlab commands you used in your summary file. Send me the plotin the form of a jpg file. To create such a file, use the commandor use the ``File' menu on the plot window, ``save' and choose JPEG image.

Exercise 2:

1. Write an m-file namedcosmx.m (``COS Minus X') that defines the function
   
   Recall that a function m-file is one that starts with the wordfunction. In this case it should start off as(This is very simple-the object of this exercise is toget going.) Include a copy of this file with your summary.
2. Now use your cosmx function to compute cosmx(0.5).Your result should be about 0.37758 and should not result in extraneousprinted or plotted information.
3. Now use your cosmx function in the same plot sequenceas above to plot cosmx and a line representing the horizontal axis on the interval. Send me the plot file as exer2.jpg.

   Hint: The horizontal axis is the line . To plot this asa function, you need a vector of at least two -values between and and a corresponding vector of -values all of whichare zero. Think of a way to generate an appropriate set of -valuesand -values when all the -values are zero. You could use theMatlab function zeros.

4. Does the value of at the point where the curve crosses the-axis agree with the previous exercise?

The idea of the bisection method is very simple. We assume that weare given two values and (), and that the function is positive at one of these values (say at ) andnegative at the other. (When this occurs, we say that isa ``change-of-sign interval' for the function.) Assuming is continuous, we know that there must be at leastone root in the interval.

Intuitively speaking, if you divide a change-of-sign interval in half, oneor the other half must end up being a change-of-sign interval, and it isonly half as long. Keep dividing in half until the change-of-sign intervalis so small that any point in it is within the specified tolerance.

More precisely, consider the point .If we are done. (This is pretty unlikely.)Otherwise, depending on the sign of , we know that theroot lies in or . In any case, ourchange-of-sign interval is now half as large as before. Repeat thisprocess with the new change of sign interval until the interval issufficiently small and declare victory.

We are guaranteed to converge. We can even compute themaximum number of steps this could take, because if the originalchange-of-sign interval has length then after one bisection the current change-of-sign interval is of length , etc. We know in advance how wellwe will approximate the root . These are very powerfulfacts, which make bisection a robust algorithm--that is, itis very hard to defeat it.

Exercise 3:

1. If we know the start points and and the interval size tolerance , we can predict beforehandthe number of steps required to reach the specified accuracy. Thebisection method will always find the root in that number or fewer steps.What is the formula for that number?
2. Give an example of a continuous function that has only one root in theinterior of the interval [-2, 1], but for which bisection could not be used.

Before we look at a sample bisection program, let's discuss someprogramming issues. If you haven't done much programming before, thisis a good time to try to understand the logic behind how we choosevariables to set up, what names to give them, and how to control thelogic.

First of all, we should write the bisection algorithm as a Matlabfunction. There are two reasons for writing it as a functioninstead of a script m-file (without the function header) or by typingeverything at the command line. The first reason is a practical one:we will be executing the algorithm several times, with differing endpoints, and functions allow us to use temporary variables withoutworrying about whether or not they have already been used before. Thesecond reason is pedantic: you should learn how to do this nowbecause it is an important tool in your toolbox.

A precise description of the bisection algorithm is presentedby Quarteroni, Sacco, and Saleri in Section 6.2.1, on pages 250-251. The algorithmcan be expressed in the following way. (Note: the product is negative if and only if and are ofopposite sign and neither is zero.)

Given a function and so that , thensequences and for with for each can be constructed bystarting out with , then, for each ,

1. Set .
2. If is small enough, or if exit the algorithm.
3. If , then set and .
4. If , then set and .
5. Return to step 1.

The bisection algorithm can be described in a manner more appropriate forcomputer implementation in the following way. Any algorithm intended fora computer program must address the issue of when to stop theiteration. In the following algorithm,

Algorithm Bisect(,)

1. Set :=.
2. If , or if happens to be exactly zero, then accept as the approximate root, and exit.
3. If sign() * sign() < 0, then :=; otherwise, :=.
4. Return to step 1.

Remarks:

• The latter form of the algorithm is expressed in a form that iseasily translated into a loop, with an explicit exit condition.
• The latter form of the algorithm employs a signum function ratherthan the value of in order to determine sign. This is a better strategybecause of roundoff errors. If and are each less than in Matlab, then their product is zero, not positive.

The language of this description, an example of ``pseudocode,' isbased on a computer language called Algol but is intended merely to bea clear means of specifying algorithms in books and papers. (The term``pseudocode' is used for any combination of a computer codinglanguage with natural language.) One objective of this lab is to showhow to rewrite algorithms like it in the Matlab language. As astarting point, let's fix to be the function cosmxthat you just wrote. Later you will learn how to add it to thecalling sequence. Let's also fix .

Here is a Matlab function that carries out the bisection algorithmfor our cosmx function. In Matlab, returned values are separated from arguments, sothe variables x and itCount are written on the left of an equal sign.

Note: I have violated my own rule abouthaving functions print things, but the disp statementmerely prints progress of the function and will be discarded once youare confident the code is correct.

Remarks about the code: Compare the code and thealgorithm. The code is similar to the algorithm, but there are somesignificant differences. For example, the algorithm generatessequences of endpoints and andmidpoints . The code keeps track of only the most recent endpoint and midpoint values.This is a typical strategy in writing code--use variables to representonly the most recent values of a sequence when the full sequence doesnot need to be retained.

A second difference is in the looping. There is an error exit in thecase that the convergence criterion is never satisfied. If thefunction ever exits with this error, it is because either yourestimate of the maximum number of iterations is too small or becausethere is a bug in your code.

Note the symbol to represent equality in a test. Use ofa single = sign causes a syntax error in Matlab. In languages such as C no syntax error is produced but the statement iswrong nonetheless and can cause serious errors that are very hard to find.

Exercise 4: Create the file bisect_cosmx.mby typing it in yourself, or by using copy-and-paste to copy it from theweb browser window.

1. Replace the ``???' with your result from Exercise 3.This value should be an integer, so be sure to increase the value fromExercise 3 to the next larger integer.(If you are not confident of your result, use 10000. If the codeis correctly written, it will work correctly for any value larger thanyour result from Exercise 3.)
   Hint: Recall that .
2. Why is the name EPSILON all capitalized?
3. The variable EPSILON and the reserved Matlabvariable eps have similar-sounding names. As EPSILONis used in this function, is it related to eps?If so, how?
4. In your own words, what does the sign(x) function do? What if x is 0? The sign function can avoid representation difficulties whenfa and fm are near the largest or smallest numbers that can be represented on the computer.
5. The disp command is used only to monitor progress. Notethe use of ... as the continuation character.
6. What is the result if the Matlab error function is called?
7. Based on your understanding of the code, what would the final value of itCount be if aand b are already closer together than the tolerance when the functionstarts?
8. Try the command: Is the value z close to a root of the equation cos(z)=z?
   Warning: You should not see the error message! If you do,your value for the maximum number of iterations is too small.
9. How many iterations were required? Is this value smaller than thevalue from your formula in Exercise 3?
10. Type the command help bisect_cosmx at the Matlabcommand line. You should see your comments reproduced as a helpmessage. This is one reason for putting comments at the top.
11. For what input numbers will this program producean incorrect answer (i.e., return a value that is not closeto a root)? Add a check before the loop so that the function willcall the Matlab error function instead of returning an incorrect answer.
12. In your opinion, why did I use abs in the convergence test?

Now we have written an m-file called bisect_cosmx.m that canfind roots of the function cosmx on any interval, but we reallywant to use it on an arbitrary function. To do this, we couldcreate an m-file for the other function, but we would have to callthat m-file cosmx.m even if there were no cosines in it, becausethat is what the bisect_cosmx function expects. Suppose that wehave three files, called f1.m through f3.m, each of whichevaluates a function that we want to use bisection on. We couldrename f1.m to cosmx.m, and change its name inside thefile as well, and repeat this for the other files. But who wants to dothat?

It is more convenient to introduce a dummy variable for the functionname, just as we used the dummy variables a and b fornumerical values. In Matlab, function names can be passed as ``function handles.'

In the bisection function given below, the name of the function is specifiedin the function call using the @ character as:

Warning: If you forget the @ in front of thefunction name, Matlab will produceerrors that seem to make no sense at all (even less than usual). Keepin mind that function names in calling sequences must havethe @ character in front of them.Check for this first when you get errors in a function.

Remark: Quarteroni, Sacco, and Saleri have an m-filenamed bisect.m on p. 252. Their function includes a variablenamed fun that contains a string representing the function whoseroot is to be found. They find function values using eval. Incontrast, bisect.m presented above uses a variable named functo contain a function handle to evaluate the function.Using function handles is moreflexible, allowing much more complicated functions, and is similar to theway that other programming languages (Fortran, C, C++, Java, etc.)work.

Exercise 5: Copy the file bisect_cosmx.mto a new file named bisect.m, or use the ``Save as' menu choiceto make a copy with the new name bisect.m.

1. In the declaration line for bisect,include a dummy function name. We might as well call it func, sothe line becomes
2. Change the comments describing the purpose of the function,and explain the variable func.
3. Where we evaluated the function by writing cosmx(x),you now write func(x). For instance, the line must be rewritten as: Make this change for fa, and a similar change for fb and fx.
4. When we use bisect, we must include the function namepreceeded by @Execute this command and make sure that you get the same resultas before. (This is called ``regression testing.')

   Warning: If you forget the @ before cosmx anduse the commandyou will get the following mysterious error message:

5. It is always a good idea to check that your answers are correct. Consider the simple problem . Write a simple function m-file named f0.m (based on cosmx.m but defining the function ) and show that you get the correct solution (x=1 to within EPSILON) starting from the change-of-sign interval [0,3].

Secant Method Vs Newton Method

Remark for advanced users: Matlab allows a simple function to be defined withouta name or an m-file. You can then use it anywhere a function handle could be used. For example, for the function , you could writeor even asBe very careful when using this last form because there is no@ appearing in the bisect call but it is a functionhandle nonetheless.

Consider the function f0=x-1 on the interval [0,3].Does your new function bisect0 find the correct answer? (It should.)

The amount the residual differs from zero is called the residualerror and the difference between the solution you found and the truesolution (found by hand or some other way) is called the trueerror. What are the residual and true errors when you used bisect0 to find the root of the function f0 above?How many iterations did it take?

What are the residual and true errors when usingbisect0 to find the root of the function f5=(x-1)^5?How many iterations did it take?

To summarize your comparison, fill in the following table

The bisect0 function has the value EPSILON = 1.0e-10;in it. Does the table show that the residual error is always smallerthan EPSILON? What about the true error?

The secant method is described by Quarteroni, Sacco, and Saleri in Section 6.2.2. Instead of dividing the interval in half, as is donein the bisection method, it regards the function as approximately linear,passing through the two points and and then finds theroot of this linear function. The interval need no longer bea change-of-sign interval, and the next value of need no longerlie inside the interval . In consequence, residual convergencemust be used in the algorithm instead of true convergence.

It is easy to see that if a straight line passes through the two points and , then it crosses the -axis () when.

The secant method can be described in the following manner.

Algorithm Secant(,)

1. Set
   
   
2. If , then accept as the approximate root, and exit.
3. Replace with and with .
4. Return to step 1.

The secant method is sometimes much faster than bisection, but sinceit does not maintain an interval inside which the solution must lie, thesecant method can fail to converge at all. Unlike bisection, thesecant method can be generalized to two or more dimensions, andthe generalization is usually called Broyden's method.

Exercise 8:

1. Write an m-file named secant.m, based partly on bisect.m, and beginning with the lines
2. Because residual error is used for convergence, and becausethe number of iterations is unknown, choose the values
3. Comment any ``disp' statements so they do not provide a distraction.
4. Complete secant.m according to the secant algorithm given above.
5. Test secant.m with the linear function f0(x)=x-1. Youshould observe convergence to the exact root in a single iteration. If youdo not, go back and correct your code. Explain why it should only takea single iteration.
6. Repeat the experiments done with bisection above, usingthe secant method, and fill in the following table.
7. In the above table, you can see that the secant method can be eitherfaster or slower than bisection. You may also observe convergence failures:either convergence to a value that is not near a root or convergence toa value substantially less accurate than expected. Regarding thebisection roots as accurate, are there any examplesof convergence to a value that is not near a root in the table? If so, which?Are there any examples of inaccurate roots? If so, which?

The regula falsi algorithm can be described in the following manner, reminiscentof the bisection and secant algorithms given above. Since the interval is supposed to begin as a change-of-sign interval, convergence isguaranteed.

Algorithm Regula(,)

1. Set
   
   
2. If , then accept as the approximate root, and exit.
3. If sign() * sign() , then set , to keep the change-of-sign interval.
4. Set .
5. Return to step 1.

Remark:It is critical to observe that Step 3 maintains theinterval as a change-of-sign interval that is a subinterval ofthe initial change-of-sign interval. Thus, regula falsi, unlike the secant method,must converge, although convergence might take a long time.

Exercise 9:

1. Starting from your bisect0.m file, write an m-file namedregula.m to carry out the regula falsi algorithm.
2. Test your work by finding the root of f0(x)=x-1 onthe interval [-1,2]. You should get the exact answer ina single iteration.
3. Repeat the experiments from the previous exercise but using regula.m and fill in the following table.You should observe both faster and slower convergence, comparedwith bisection and secant. You should not observe lack of convergenceor convergence to an incorrect solution.
4. Function f5, (x-1)^5 turns out to be verydifficult for regula falsi! Loosen your convergence criterion to a tolerance of and increase the maximum allowable number ofiterations to 500,000 and fill in the following line of the table.(Be sure that any ``disp' statements are commented out.)You should observe convergence, but it takes a very large number of iterations.
5. regula.m and bisect.m both keep the current iterate,x in a change-of-sign interval. Why would it be wrong to usethe same convergence criterion in regula.m as was used inbisect.m?

In the secant method, the function is approximated by its secant (thelinear function passing through the two endpoints) and then thenext iterate is taken to be the root of this linear function. Muller'smethod goes one better and approximates the function by the quadraticfunction passing through the two endpoints and the current iterate. Thenext iterate is then taken to be the root of this quadratic that isclosest to the current iterate.

In more detail, suppose I want a function with signature similar toour functions above:

1. Choose a convergence criterion and a maximumnumber of iterations ITMAX = 100
2. Choose a third point, not quite at the center of the change-of-signinterval [a,b]
3. Evaluate y0, y1, and y2 as values offunc at x0, x1, and x2.
4. Determine coefficients A, B, andC of the polynomial passing through the three points(x0,y0), (x1,y1) and (x2,y2), :
5. If the polynomial has real roots, find them and choose theone closer to x2, otherwise choose something reasonable:
6. Discard the point (x0, y0) and add the new point
7. Exit the loop and return result = x2 if the residual (y2) is smaller than .
8. Go back to 4 unless ITMAXis exceeded, in which case print an error message.

Exercise 10:

1. Following the outline above, create a file muller.m thatcarries out Muller's method.
2. Test muller.m on the simple linear function starting with thechange-of-sign interval [0,2]. You should observe convergencein a single iteration. If you do not observe convergence in a singleinteration, you probably have a bug. Fix it before continuing. Explain why (one sentence, please) it requires only a single iteration.
3. Test muller.m on the function f1 on the change-of-sign interval [0, 5]. Again, you should observe convergence in a single iteration. If you do not, you probably have abug. Fix it before continuing.Explain why (one sentence, please) it requires only a single iteration.
4. Repeat the experiments from the previous exercise but using muller.m and fill in the following table.

   You should observe that, since the interval at each iteration need notbe a change-of-sign interval, it is possible for the method to fail, as with secant. It is possible to combine several methods insuch a way that the speed of secant or muller ispartially retained but failure is avoided by maintaining achange-of-sign interval at each iteration. One such method is calledBrent's method. See, for example,http://en.wikipedia.org/wiki/Brent%27s_method

Back to MATH2070 page.

Secant Method Algorithm

Secant Method Calculator

Mike Sussman2016-08-26