The main report consists of the following sections:

• Verification Section - My secant method and the driving routines that call it are tested on 4 different functions, and compared to the Web-Based Secant Solver
• Analysis Section - The behavior of my secant.m code is tested with oscillatory and potential divide-by-zero non-convergent cases, and the convergence rate of the secant method is explored and contrasted with the bisection method
• Further Explorations -A quick comparison is made between the 3 convergent criteria for the secant method that were mentioned in class
• Reference Section - Links are given to all of the significant files used in the project, along with short functional descriptions where appropriate

 

 
 
 
 
 
 
 

Verification Section
In order to convince myself that  my secant.m file was behaving properly, I executed the basic run_me.m function to call differing functions, and then compared the results to the expected answers. I chose the most 'interesting of the examples to create my movie from, which turned out to be the 3rd example in the list. The examples I used are:

1. Problem from textbook page 125
• The original results
• Results after modifying secant  function to match output from web-based solver
  
2. Problem from previous Test
3. Problem from textbook, page 105
4. Problem from textbook page 101
5. Movie that was produced during verification (ver_pg105_movie.avi)

Plot of the function:

secant.m output:

=====================================================================================

Math4650 - Project 1 - Secant Method Matlab Implementation - Mike Lewis

Input Parameters
-------------------------------------------------------------------------------------
           file = f_pg125.m
       function = f(x) = x^5 + x^3 + 3
            x_0 = +1.000000e+000
            x_1 = -1.000000e+000
 Max Iterations =  100
      Tolerance =  1.000000e-008

Start Values (adjusted for polarity if required):
-------------------------------------------------------------------------------------
   x_0 = -1.000000e+000
   x_1 = +1.000000e+000
f(x_0) = +1.000000e+000
f(x_1) = +5.000000e+000
 

------------------------------ Function Output --------------------------------------
  Initial x0    x = -1.000000e+000    f(x_0) = +1.000000e+000
  Initial x1    x = +1.000000e+000    f(x_1) = +5.000000e+000
Iteration  1    x = -1.500000e+000    f(x_0) = -7.968750e+000
Iteration  2    x = -1.055749e+000    f(x_0) = +5.116487e-001
Iteration  3    x = -1.114158e+000    f(x_0) = -9.991207e-002
Iteration  4    x = -1.104615e+000    f(x_0) = +7.593217e-003
Iteration  5    x = -1.105289e+000    f(x_0) = +1.011787e-004
Iteration  6    x = -1.105299e+000    f(x_0) = -1.045169e-007
Iteration  7    x = -1.105299e+000    f(x_0) = +1.436185e-012

---------------------- Halt Due To Convergence Limit! -------------------------------

    root = -1.105299e+000    f(x_0) = 1.436185e-012
======================================================================================

Secant Solver  on the Web results: Run it yourself: Secant Solver on the Web

x*x*x*x*x+x*x*x+3 = 0
n	x
0	-1
1	1
2	-1.5
3	0.03614457831325302
4	-0.38400233921340904
5	-19.349266937033033
6	-0.38402280789503607
7	-0.38404327647589853
8	-5.70881599561185
9	-0.3865439592063016
10	-0.3890422849830717
11	-5.585863762891204
12	-0.3917573808811774
13	-0.3944696182413761
14	-5.418342008099811
15	-0.3975167811771767
16	-0.4005602171324961
17	-5.240163796091338
18	-0.40401759378038715
19	-0.4074699851557766
20	-5.048551119912421
21	-0.41144767124239934
22	-0.41541846858630543
23	-4.841074846865122
24	-0.4200734034809859
25	-0.42471843078936866
26	-4.614549263944573
27	-0.4302834369193853
28	-0.4358334581564652
29	-4.364745070627189
30	-0.442671094020024
31	-0.44948454611068894
32	-4.085912945668211
33	-0.45819474922210546
34	-0.4668624082500378
35	-3.770051160791487
36	-0.47851796730561924
37	-0.49008924610174964
38	-3.405827675330812
39	-0.5068118943443367
40	-0.5233365119160505
41	-2.9773616228549145
42	-0.5499104665610384
43	-0.5758828978953849
44	-2.4652310059570004
45	-0.6249239272893554
46	-0.6712498483006798
47	-1.8653672121919507
48	-0.7780506124120914
49	-0.8642026116794079
50	-1.298648897771246
51	-1.0352400683077987
52	-1.0867749986207573
53	-1.107392174190506
54	-1.1052392281314827
55	-1.1052983585732785
56	-1.1052985460229752
57	-1.1052985460061695
58	-1.1052985460061695

Summary of Results:

• Book result is given as -1.10530, which is within 1 place in 10^-6, difference due to the different number of digits displayed in the answer.
• Amazingly, the secant solver on the web got essentially the same answer, but took a lot more iterations to do it! This is a very unexpected result! I decided to try to figure out why. The professor hinted that the problem had to be in the code that swaps the order of the input variables, based on which is greater. Possibly this code did not exist in the web solver. So, I edited a copy of my secant.m file, calling the new file secant_modified.m. I also needed to create a new calling function, run_me_modified.m, the only change required being to call the new instead of the old secant function. Then I re-ran using the same input arguments:

Plot of the modified function:

Summary of Results:

• Well, that was the difference. Now the results agree, allowing for a slight difference in numbering of samples (my numbering is one higher since I list both starting values), and allowing for a different ordering of the values at the very start. The discrpancies are  due to differing precision of the display used between the two programs, except for the few aded entries shown at the end of the Web Solver output. It appears that the Web Solver is using a different criteria for determining convergence, which explains the extra samples at the end. Since my version seemed to get the correct value with less iterations, I have decided to use the unmodified versions of secant.m and run_me.m for the remainder of the report, unless I find a reason later to change back...

back to top of page

Plot of the Function

secant.m output:

=====================================================================================

Math4650 - Project 1 - Secant Method Matlab Implementation - Mike Lewis

Input Parameters
-------------------------------------------------------------------------------------
           file = f_prevtest.m
       function = f(x) = (x-1)*(x-3)*(x-5)*(x-10)*(x-12)
            x_0 = +4.000000e+000
            x_1 = +8.000000e+000
 Max Iterations =  100
      Tolerance =  1.000000e-008

Start Values (adjusted for polarity if required):
-------------------------------------------------------------------------------------
   x_0 = +4.000000e+000
   x_1 = +8.000000e+000
f(x_0) = -1.440000e+002
f(x_1) = +8.400000e+002
 

------------------------------ Function Output --------------------------------------
  Initial x0    x = +4.000000e+000    f(x_0) = -1.440000e+002
  Initial x1    x = +8.000000e+000    f(x_1) = +8.400000e+002
Iteration  1    x = +4.585366e+000    f(x_0) = -9.462086e+001
Iteration  2    x = +5.707050e+000    f(x_0) = +2.433916e+002
Iteration  3    x = +4.899362e+000    f(x_0) = -2.699501e+001
Iteration  4    x = +5.024704e+000    f(x_0) = +6.986288e+000
Iteration  5    x = +4.998935e+000    f(x_0) = -2.981128e-001
Iteration  6    x = +4.999989e+000    f(x_0) = -2.953363e-003
Iteration  7    x = +5.000000e+000    f(x_0) = +1.281698e-006
Iteration  8    x = +5.000000e+000    f(x_0) = -5.471179e-012

---------------------- Halt Due To Convergence Limit! -------------------------------

    root = 5.000000e+000    f(x_0) = -5.471179e-012
======================================================================================

Secant Solver  on the Web results: Run it yourself: Secant Solver on the Web

(x-1)*(x-3)*(x-5)*(x-10)*(x-12) = 0
n	x
0	4
1	8
2	4.585365853658536
3	4.931062842485611
4	5.016502339782172
5	4.999520311029412
6	4.999996810958355
7	5.000000000623025
8	4.999999999999999
9	5
10	5

Summary of Results:

• This was the problem from a previous test. Due to the formulation of the function, it is easy to see that the root found is a correct root.
• The web solver result is in total agreement here for final result, and iterations (accounting for my differing counting scheme and counting the output line as the final iteration), but some subtle differences in the values for x along way which are the result of restoring the sign checking into the secant function, as shown in the last test before this one.

back to top of page

Results For Problem from page 105, for input function x³- 2x² + x - 3:

Plot of Function:

secant.m output :

=====================================================================================

Math4650 - Project 1 - Secant Method Matlab Implementation - Mike Lewis

Input Parameters
-------------------------------------------------------------------------------------
           file = f_pg105.m
       function = f(x) =  x^3 - (2 * x^2) + x - 3
            x_0 = +4.000000e+000
            x_1 = +1.000000e+000
 Max Iterations =  100
      Tolerance =  1.000000e-008

Start Values (adjusted for polarity if required):
-------------------------------------------------------------------------------------
   x_0 = +1.000000e+000
   x_1 = +4.000000e+000
f(x_0) = -3.000000e+000
f(x_1) = +3.300000e+001
 

------------------------------ Function Output --------------------------------------
  Initial x0    x = +1.000000e+000    f(x_0) = -3.000000e+000
  Initial x1    x = +4.000000e+000    f(x_1) = +3.300000e+001
Iteration  1    x = +1.250000e+000    f(x_0) = -2.921875e+000
Iteration  2    x = +1.060000e+001    f(x_0) = +9.738960e+002
Iteration  3    x = +1.277968e+000    f(x_0) = -2.901256e+000
Iteration  4    x = +5.213337e+000    f(x_0) = +8.954824e+001
Iteration  5    x = +1.401468e+000    f(x_0) = -2.774116e+000
Iteration  6    x = +4.096162e+000    f(x_0) = +3.626670e+001
Iteration  7    x = +1.592944e+000    f(x_0) = -2.439948e+000
Iteration  8    x = +2.991019e+000    f(x_0) = +8.856872e+000
Iteration  9    x = +1.894908e+000    f(x_0) = -1.482443e+000
Iteration 10    x = +2.362419e+000    f(x_0) = +1.385089e+000
Iteration 11    x = +2.136600e+000    f(x_0) = -2.398147e-001
Iteration 12    x = +2.169928e+000    f(x_0) = -2.995296e-002
Iteration 13    x = +2.174685e+000    f(x_0) = +8.118269e-004
Iteration 14    x = +2.174559e+000    f(x_0) = -2.627157e-006
Iteration 15    x = +2.174559e+000    f(x_0) = -2.291780e-010

---------------------- Halt Due To Convergence Limit! -------------------------------

    root = 2.174559e+000    f(x_0) = -2.291780e-010

======================================================================================
Secant Solver  on the Web results: Run it yourself: Secant Solver on the Web

x*x*x- 2*x*x+ x - 3 = 0
n	x
0	1
1	4
2	1.25
3	1.4736842105263157
4	3.838061646986655
5	1.6800488695274716
6	1.8392327759609992
7	2.3626354666543623
8	2.1278301005834734
9	2.1688300953974173
10	2.174750846910788
11	2.174558643528284
12	2.1745594101906423
13	2.17455941029298
14	2.17455941029298

Summary of Results:

• The result for the root agrees with the book to the number of decimal places provided by secant.m.
• This run also seems to illustrate the differences in convergence speed between the newton method, which is used with this data in the book, and the secant method, used here. My algorithm took 15 iterations after the initial points to determine the answer, while the newton example indicates just 6 iterations after the initial guess. There is probably some effect of the second point value required by the secant method also affecting the convergence rate. Since this example did have 15 iterations, I thought it ws the most interesting movie generated by all of the 4 verification runs, so I  used it for my movie case.

• If you did not load it from the top of the page, you can also load it here: ver_pg105_movie.avi
• With regard to the web solver, we again reach the same answer (allowing for my lower precision display), but we get there differently, again due to the use of the sign check on the inputs after each iteration that is used in my implementation. In this case, the web solver uses less iterations then my approach. The difference is not nearly as striking as the original difference I discovered, however, so I don't think it warrants discarding my approach.

Plot of Function:

secant.m output:

=====================================================================================

Math4650 - Project 1 - Secant Method Matlab Implementation - Mike Lewis

Input Parameters
-------------------------------------------------------------------------------------
           file = f_pg101.m
       function = f(x) = x^3 + (2 * x^2) + (10 * x) - 20
            x_0 = +1.000000e+000
            x_1 = +2.000000e+000
 Max Iterations =  100
      Tolerance =  1.000000e-008

Start Values (adjusted for polarity if required):
-------------------------------------------------------------------------------------
   x_0 = +1.000000e+000
   x_1 = +2.000000e+000
f(x_0) = -7.000000e+000
f(x_1) = +1.600000e+001
 

------------------------------ Function Output --------------------------------------
  Initial x0    x = +1.000000e+000    f(x_0) = -7.000000e+000
  Initial x1    x = +2.000000e+000    f(x_1) = +1.600000e+001
Iteration  1    x = +1.304348e+000    f(x_0) = -1.334758e+000
Iteration  2    x = +1.376054e+000    f(x_0) = +1.531733e-001
Iteration  3    x = +1.368672e+000    f(x_0) = -2.872217e-003
Iteration  4    x = +1.368808e+000    f(x_0) = -6.018647e-006
Iteration  5    x = +1.368808e+000    f(x_0) = +2.372076e-010

---------------------- Halt Due To Convergence Limit! -------------------------------

    root = 1.368808e+000    f(x_0) = 2.372076e-010
======================================================================================
Secant Solver  on the Web results: Run it yourself: Secant Solver on the Web

(x*x*x)+(2*x*x)+(10*x)-20 = 0
n	x
0	1
1	2
2	1.3043478260869565
3	1.3579123046578667
4	1.3690133259925659
5	1.3688074597219246
6	1.3688081077828753
7	1.3688081078213725
8	1.3688081078213725

Summary of Results:

• The result again agrees with the book, although I have no intermediate terms to compare to this time, having just lifted the problem and the answer from the book on page 101 (see problem 8).. The web solver agrees with the result as well.

back to top of page

1. Behavior of the secant.m function with data that is not well-behaved to see what happens to the method (and my program) when the algorithm does not converge
2.  What happens when the denominator term, (f(x_n) - f(x_n-1), approaches zero
3. What can we say about the rate of convergence of the secant method?

• Convergence Example 1, Another Look at The Page 125 Function
• Convergence Example 2, Another Look at The Page 105 Function
• A Comparison of The Convergence  Rate for the Secant and Bisection Methods, Using the Page 125 Function

Plot of function:

secant.m output:
=====================================================================================

Math4650 - Project 1 - Secant Method Matlab Implementation - Mike Lewis

Input Parameters
-------------------------------------------------------------------------------------
           file = f_pg105.m
       function = f(x) =  x^3 - (2 * x^2) + x - 3
            x_0 = +0.000000e+000
            x_1 = +8.000000e+000
 Max Iterations =  100
      Tolerance =  1.000000e-008

Start Values (adjusted for polarity if required):
-------------------------------------------------------------------------------------
   x_0 = +0.000000e+000
   x_1 = +8.000000e+000
f(x_0) = -3.000000e+000
f(x_1) = +3.890000e+002
 

------------------------------ Function Output --------------------------------------
  Initial x0    x = +0.000000e+000    f(x_0) = -3.000000e+000
  Initial x1    x = +8.000000e+000    f(x_1) = +3.890000e+002
Iteration  1    x = +6.122449e-002    f(x_0) = -2.946043e+000
Iteration  2    x = +3.404064e+000    f(x_0) = +1.667387e+001
Iteration  3    x = +5.631710e-001    f(x_0) = -2.892536e+000
Iteration  4    x = +2.769794e+001    f(x_0) = +1.973954e+004
Iteration  5    x = +5.671466e-001    f(x_0) = -2.893738e+000
Iteration  6    x = -9.001529e+000    f(x_0) = -9.034281e+002
Iteration  7    x = +5.938930e-001    f(x_0) = -2.902053e+000
Iteration  8    x = -8.773788e+000    f(x_0) = -8.411331e+002
Iteration  9    x = +5.953903e-001    f(x_0) = -2.902529e+000
Iteration 10    x = -8.762631e+000    f(x_0) = -8.381572e+002
Iteration 11    x = +5.954664e-001    f(x_0) = -2.902553e+000
Iteration 12    x = -8.762067e+000    f(x_0) = -8.380071e+002
Iteration 13    x = +5.954703e-001    f(x_0) = -2.902555e+000
Iteration 14    x = -8.762039e+000    f(x_0) = -8.379996e+002
Iteration 15    x = +5.954705e-001    f(x_0) = -2.902555e+000
Iteration 16    x = -8.762037e+000    f(x_0) = -8.379992e+002                                               *
                                               *
                                               *(skip redundant data, results oscillating between identical points)                                               *
                                               *
                                               *Iteration 95    x = +5.954705e-001    f(x_0) = -2.902555e+000
Iteration 96    x = -8.762037e+000    f(x_0) = -8.379992e+002
Iteration 97    x = +5.954705e-001    f(x_0) = -2.902555e+000
Iteration 98    x = -8.762037e+000    f(x_0) = -8.379992e+002

----------------------- Halt Due To Iteration Limit! --------------------------------
Final attempt at root = -8.762037e+000
         Final f(x_0) = -8.379992e+002
======================================================================================
Summary of Results:
So, just changing the input parameters from (1,4) to (0,8), caused this to not converge! Including 0 on the low side may be the problematic part. It is clear from the plot that the function has what could be a 'flat-spot' in the area of interest - i.e. a place where the first derivitive is 0. This is most likely the cause of the trouble. One thing is for sure, we are not going to find a root this way...
You can refer back to the earlier run from here for comparison:  (pg105 using more well behaved starting points)

back to top of page

Plot of function:

secant.m output:
=====================================================================================

Math4650 - Project 1 - Secant Method Matlab Implementation - Mike Lewis

Input Parameters
-------------------------------------------------------------------------------------
           file = f_pg105.m
       function = f(x) =  x^3 - (2 * x^2) + x - 3
            x_0 = +0.000000e+000
            x_1 = +4.000000e+000
 Max Iterations =  100
      Tolerance =  1.000000e-008

Start Values (adjusted for polarity if required):
-------------------------------------------------------------------------------------
   x_0 = +0.000000e+000
   x_1 = +4.000000e+000
f(x_0) = -3.000000e+000
f(x_1) = +3.300000e+001
 

------------------------------ Function Output --------------------------------------
  Initial x0    x = +0.000000e+000    f(x_0) = -3.000000e+000
  Initial x1    x = +4.000000e+000    f(x_1) = +3.300000e+001
Iteration  1    x = +3.333333e-001    f(x_0) = -2.851852e+000
Iteration  2    x = +6.750000e+000    f(x_0) = +2.201719e+002
Iteration  3    x = +4.153846e-001    f(x_0) = -2.858032e+000
Iteration  4    x = -3.753038e+001    f(x_0) = -5.572022e+004
Iteration  5    x = +3.352714e-001    f(x_0) = -2.851856e+000
Iteration  6    x = -1.474048e+003    f(x_0) = -3.207187e+009
Iteration  7    x = +3.333346e-001    f(x_0) = -2.851852e+000

----------------------- Halt Due To Denominator < sqr(eps)! --------------------------------
Final attempt at root = +3.333333e-001
         Final f(x_0) = -2.851852e+000
      f(x_1) - f(x_0) = -1.718625e-012
======================================================================================

Summary of Results:
Now we have changed the original input parameters from (1,4) to (0,4), and we have what would potentially have resulted in a division by zero condition! Clearly, this function is a bit problematic for the secant method (or at least my implementation of it) around this region. The denominator is at  -1.718625e-012, getting pretty small and our root is not close to the other root found. The square root of machine epsilon ((eps)^0.5 = 1.4901e-008) is used as a test criteria for divide-by-zero, since, by definition CPU calculations become unpredictable near machine epsilon itself (eps for this system according to Matlab = 2.2204e-016, a very familiar result from class).

You can refer back to the earlier run from here for comparison:  (pg105 using more well behaved starting points)

back to top of page

Convergence Plot For Problem from page 125, for input function x⁵ + x³ + 3,
Plot of Function:

Plot of Convergence:

secant.m output:
=====================================================================================

Math4650 - Project 1 - Secant Method Matlab Implementation - Mike Lewis

Input Parameters
-------------------------------------------------------------------------------------
           file = f_pg125.m
       function = f(x) = x^5 + x^3 + 3
            x_0 = +0.000000e+000
            x_1 = +1.000000e+000
 Max Iterations =  100
      Tolerance =  1.000000e-014

Start Values (adjusted for polarity if required):
-------------------------------------------------------------------------------------
   x_0 = +0.000000e+000
   x_1 = +1.000000e+000
f(x_0) = +3.000000e+000
f(x_1) = +5.000000e+000
 

------------------------------ Function Output --------------------------------------
  Initial x0    x = +0.000000e+000    f(x_0) = +3.000000e+000
  Initial x1    x = +1.000000e+000    f(x_1) = +5.000000e+000
Iteration  1    x = -1.500000e+000    f(x_0) = -7.968750e+000
Iteration  2    x = -4.102564e-001    f(x_0) = +2.919328e+000
Iteration  3    x = -1.525641e+001    f(x_0) = -8.300834e+005
Iteration  4    x = -4.103086e-001    f(x_0) = +2.919294e+000
Iteration  5    x = -4.924636e+000    f(x_0) = -3.012913e+003
Iteration  6    x = -4.146784e-001    f(x_0) = +2.916431e+000
Iteration  7    x = -4.865734e+000    f(x_0) = -2.839553e+003
Iteration  8    x = -4.192453e-001    f(x_0) = +2.913359e+000
Iteration  9    x = -4.749976e+000    f(x_0) = -2.522175e+003
Iteration 10    x = -4.242420e-001    f(x_0) = +2.909902e+000
Iteration 11    x = -4.630488e+000    f(x_0) = -2.225078e+003
Iteration 12    x = -4.297356e-001    f(x_0) = +2.905984e+000
Iteration 13    x = -4.504413e+000    f(x_0) = -1.942739e+003
Iteration 14    x = -4.358215e-001    f(x_0) = +2.901497e+000
Iteration 15    x = -4.371016e+000    f(x_0) = -1.676068e+003
Iteration 16    x = -4.426220e-001    f(x_0) = +2.896295e+000
Iteration 17    x = -4.229304e+000    f(x_0) = -1.425795e+003
Iteration 18    x = -4.502985e-001    f(x_0) = +2.890179e+000
Iteration 19    x = -4.078077e+000    f(x_0) = -1.192738e+003
Iteration 20    x = -4.590679e-001    f(x_0) = +2.882866e+000
Iteration 21    x = -3.915866e+000    f(x_0) = -9.777898e+002
Iteration 22    x = -4.692298e-001    f(x_0) = +2.873939e+000
Iteration 23    x = -3.740857e+000    f(x_0) = -7.819304e+002
Iteration 24    x = -4.812105e-001    f(x_0) = +2.862766e+000
Iteration 25    x = -3.550790e+000    f(x_0) = -6.062178e+002
Iteration 26    x = -4.956379e-001    f(x_0) = +2.848333e+000
Iteration 27    x = -3.342829e+000    f(x_0) = -4.517724e+002
Iteration 28    x = -5.134764e-001    f(x_0) = +2.828923e+000
Iteration 29    x = -3.113433e+000    f(x_0) = -3.197283e+002
Iteration 30    x = -5.362788e-001    f(x_0) = +2.801413e+000
Iteration 31    x = -2.858271e+000    f(x_0) = -2.111245e+002
Iteration 32    x = -5.666859e-001    f(x_0) = +2.759578e+000
Iteration 33    x = -2.572461e+000    f(x_0) = -1.266768e+002
Iteration 34    x = -6.094489e-001    f(x_0) = +2.689555e+000
Iteration 35    x = -2.251948e+000    f(x_0) = -6.633537e+001
Iteration 36    x = -6.734489e-001    f(x_0) = +2.556045e+000
Iteration 37    x = -1.898728e+000    f(x_0) = -2.852347e+001
Iteration 38    x = -7.742185e-001    f(x_0) = +2.257748e+000
Iteration 39    x = -1.536922e+000    f(x_0) = -9.205888e+000
Iteration 40    x = -9.244318e-001    f(x_0) = +1.534895e+000
Iteration 41    x = -1.243392e+000    f(x_0) = -1.894262e+000
Iteration 42    x = -1.067199e+000    f(x_0) = +4.002630e-001
Iteration 43    x = -1.117563e+000    f(x_0) = -1.390262e-001
Iteration 44    x = -1.104579e+000    f(x_0) = +7.995566e-003
Iteration 45    x = -1.105285e+000    f(x_0) = +1.471338e-004
Iteration 46    x = -1.105299e+000    f(x_0) = -1.600563e-007
Iteration 47    x = -1.105299e+000    f(x_0) = +3.200107e-012
Iteration 48    x = -1.105299e+000    f(x_0) = +4.440892e-016

---------------------- Halt Due To Convergence Limit! -------------------------------

back to top of page

Plot of Function:

Plot of Convergence:

=====================================================================================

Math4650 - Project 1 - Secant Method Matlab Implementation - Mike Lewis

Input Parameters
-------------------------------------------------------------------------------------
           file = f_pg105.m
       function = f(x) =  x^3 - (2 * x^2) + x - 3
            x_0 = +9.000000e+001
            x_1 = +8.000000e+001
 Max Iterations =  100
      Tolerance =  1.000000e-014

Start Values (adjusted for polarity if required):
-------------------------------------------------------------------------------------
   x_0 = +8.000000e+001
   x_1 = +9.000000e+001
f(x_0) = +4.992770e+005
f(x_1) = +7.128870e+005
 

------------------------------ Function Output --------------------------------------
  Initial x0    x = +8.000000e+001    f(x_0) = +4.992770e+005
  Initial x1    x = +9.000000e+001    f(x_1) = +7.128870e+005
Iteration  1    x = +5.662670e+001    f(x_0) = +1.752187e+005
Iteration  2    x = +4.398873e+001    f(x_0) = +8.128955e+004
Iteration  3    x = +3.305140e+001    f(x_0) = +3.395043e+004
Iteration  4    x = +2.520741e+001    f(x_0) = +1.476851e+004
Iteration  5    x = +1.916818e+001    f(x_0) = +6.324090e+003
Iteration  6    x = +1.464536e+001    f(x_0) = +2.723906e+003
Iteration  7    x = +1.122339e+001    f(x_0) = +1.170041e+003
Iteration  8    x = +8.646679e+000    f(x_0) = +5.025860e+002
Iteration  9    x = +6.706449e+000    f(x_0) = +2.153859e+002
Iteration 10    x = +5.251373e+000    f(x_0) = +9.191422e+001
Iteration 11    x = +4.168192e+000    f(x_0) = +3.883797e+001
Iteration 12    x = +3.375586e+000    f(x_0) = +1.604980e+001
Iteration 13    x = +2.817350e+000    f(x_0) = +6.305030e+000
Iteration 14    x = +2.456162e+000    f(x_0) = +2.208063e+000
Iteration 15    x = +2.261499e+000    f(x_0) = +5.989060e-001
Iteration 16    x = +2.189049e+000    f(x_0) = +9.495726e-002
Iteration 17    x = +2.175397e+000    f(x_0) = +5.437503e-003
Iteration 18    x = +2.174568e+000    f(x_0) = +5.450293e-005
Iteration 19    x = +2.174559e+000    f(x_0) = +3.181802e-008
Iteration 20    x = +2.174559e+000    f(x_0) = +1.865175e-013
Iteration 21    x = +2.174559e+000    f(x_0) = -4.440892e-016

---------------------- Halt Due To Convergence Limit! -------------------------------

    root = 2.174559e+000    f(x_0) = -4.440892e-016
======================================================================================
 
 

Comparison of Convergence of the Secant and Bisection Methods
Plot of secant function

Plot of Bisection function

 

Comparison of Convergence Rate for Both the Secant and Bisection Methods:

secant.m output:
=====================================================================================

Math4650 - Project 1 - Secant Method Matlab Implementation - Mike Lewis

Input Parameters
-------------------------------------------------------------------------------------
           file = f_pg125.m
       function = f(x) = x^5 + x^3 + 3
            x_0 = -8.000000e+000
            x_1 = +1.000000e+000
 Max Iterations =  100
      Tolerance =  1.000000e-014

Start Values (adjusted for polarity if required):
-------------------------------------------------------------------------------------
   x_0 = +1.000000e+000
   x_1 = -8.000000e+000
f(x_0) = +5.000000e+000
f(x_1) = -3.327700e+004
 

------------------------------ Function Output --------------------------------------
  Initial x0    x = +1.000000e+000    f(x_0) = +5.000000e+000
  Initial x1    x = -8.000000e+000    f(x_1) = -3.327700e+004
Iteration  1    x = +9.986479e-001    f(x_0) = +4.989207e+000
Iteration  2    x = +3.736253e-001    f(x_0) = +3.059437e+000
Iteration  3    x = -6.172791e-001    f(x_0) = +2.675175e+000
Iteration  4    x = -7.515808e+000    f(x_0) = -2.440316e+004
Iteration  5    x = -6.180353e-001    f(x_0) = +2.673760e+000
Iteration  6    x = -2.046162e+000    f(x_0) = -4.143421e+001
Iteration  7    x = -7.046062e-001    f(x_0) = +2.476511e+000
Iteration  8    x = -1.791529e+000    f(x_0) = -2.120528e+001
Iteration  9    x = -8.182706e-001    f(x_0) = +2.085266e+000
Iteration 10    x = -1.424082e+000    f(x_0) = -5.745044e+000
Iteration 11    x = -9.796024e-001    f(x_0) = +1.157865e+000
Iteration 12    x = -1.181026e+000    f(x_0) = -9.450385e-001
Iteration 13    x = -1.090507e+000    f(x_0) = +1.609611e-001
Iteration 14    x = -1.103680e+000    f(x_0) = +1.796325e-002
Iteration 15    x = -1.105335e+000    f(x_0) = -4.074499e-004
Iteration 16    x = -1.105298e+000    f(x_0) = +9.976886e-007
Iteration 17    x = -1.105299e+000    f(x_0) = +5.521095e-011
Iteration 18    x = -1.105299e+000    f(x_0) = +4.440892e-016

---------------------- Halt Due To Convergence Limit! -------------------------------

    root = -1.105299e+000    f(x_0) = 4.440892e-016
======================================================================================

bisection.m output:
=====================================================================================

Math4650 - The A. Knyazev Bisection Method Matlab Implementation - Modified To Use Snapshot
Used By Permission
Initial values: a= -8.000000e+000 , b= 1.000000e+000 , f(a)= -3.327700e+004 ,f(b)= 5.000000e+000

Iteration n= 1 , c= -3.500000e+000 , f(c)= -5.650938e+002 ,bisection_error= 4.500000e+000
Iteration n= 2 , c= -1.250000e+000 , f(c)= -2.004883e+000 ,bisection_error= 2.250000e+000
Iteration n= 3 , c= -1.250000e-001 , f(c)= 2.998016e+000 ,bisection_error= 1.125000e+000
Iteration n= 4 , c= -6.875000e-001 , f(c)= 2.521459e+000 ,bisection_error= 5.625000e-001
Iteration n= 5 , c= -9.687500e-001 , f(c)= 1.237636e+000 ,bisection_error= 2.812500e-001
Iteration n= 6 , c= -1.109375e+000 , f(c)= -4.564164e-002 ,bisection_error= 1.406250e-001
Iteration n= 7 , c= -1.039063e+000 , f(c)= 6.669962e-001 ,bisection_error= 7.031250e-002
Iteration n= 8 , c= -1.074219e+000 , f(c)= 3.299894e-001 ,bisection_error= 3.515625e-002
Iteration n= 9 , c= -1.091797e+000 , f(c)= 1.472078e-001 ,bisection_error= 1.757813e-002
Iteration n= 10 , c= -1.100586e+000 , f(c)= 5.206799e-002 ,bisection_error= 8.789063e-003
Iteration n= 11 , c= -1.104980e+000 , f(c)= 3.537743e-003 ,bisection_error= 4.394531e-003
Iteration n= 12 , c= -1.107178e+000 , f(c)= -2.097038e-002 ,bisection_error= 2.197266e-003
Iteration n= 13 , c= -1.106079e+000 , f(c)= -8.695983e-003 ,bisection_error= 1.098633e-003
Iteration n= 14 , c= -1.105530e+000 , f(c)= -2.574042e-003 ,bisection_error= 5.493164e-004
Iteration n= 15 , c= -1.105255e+000 , f(c)= 4.831192e-004 ,bisection_error= 2.746582e-004
Iteration n= 16 , c= -1.105392e+000 , f(c)= -1.045144e-003 ,bisection_error= 1.373291e-004
Iteration n= 17 , c= -1.105324e+000 , f(c)= -2.809331e-004 ,bisection_error= 6.866455e-005
Iteration n= 18 , c= -1.105289e+000 , f(c)= 1.011129e-004 ,bisection_error= 3.433228e-005
Iteration n= 19 , c= -1.105307e+000 , f(c)= -8.990518e-005 ,bisection_error= 1.716614e-005
Iteration n= 20 , c= -1.105298e+000 , f(c)= 5.605078e-006 ,bisection_error= 8.583069e-006
Iteration n= 21 , c= -1.105302e+000 , f(c)= -4.214974e-005 ,bisection_error= 4.291534e-006
Iteration n= 22 , c= -1.105300e+000 , f(c)= -1.827225e-005 ,bisection_error= 2.145767e-006
Iteration n= 23 , c= -1.105299e+000 , f(c)= -6.333569e-006 ,bisection_error= 1.072884e-006
Iteration n= 24 , c= -1.105299e+000 , f(c)= -3.642405e-007 ,bisection_error= 5.364418e-007
Iteration n= 25 , c= -1.105298e+000 , f(c)= 2.620420e-006 ,bisection_error= 2.682209e-007
Iteration n= 26 , c= -1.105298e+000 , f(c)= 1.128090e-006 ,bisection_error= 1.341105e-007
Iteration n= 27 , c= -1.105299e+000 , f(c)= 3.819248e-007 ,bisection_error= 6.705523e-008
Iteration n= 28 , c= -1.105299e+000 , f(c)= 8.842187e-009 ,bisection_error= 3.352761e-008
Iteration n= 29 , c= -1.105299e+000 , f(c)= -1.776992e-007 ,bisection_error= 1.676381e-008
Iteration n= 30 , c= -1.105299e+000 , f(c)= -8.442848e-008 ,bisection_error= 8.381903e-009
Iteration n= 31 , c= -1.105299e+000 , f(c)= -3.779315e-008 ,bisection_error= 4.190952e-009
Iteration n= 32 , c= -1.105299e+000 , f(c)= -1.447548e-008 ,bisection_error= 2.095476e-009
Iteration n= 33 , c= -1.105299e+000 , f(c)= -2.816646e-009 ,bisection_error= 1.047738e-009
Iteration n= 34 , c= -1.105299e+000 , f(c)= 3.012771e-009 ,bisection_error= 5.238689e-010
Iteration n= 35 , c= -1.105299e+000 , f(c)= 9.806245e-011 ,bisection_error= 2.619345e-010
Iteration n= 36 , c= -1.105299e+000 , f(c)= -1.359291e-009 ,bisection_error= 1.309672e-010
Iteration n= 37 , c= -1.105299e+000 , f(c)= -6.306147e-010 ,bisection_error= 6.548362e-011
Iteration n= 38 , c= -1.105299e+000 , f(c)= -2.662759e-010 ,bisection_error= 3.274181e-011
Iteration n= 39 , c= -1.105299e+000 , f(c)= -8.410694e-011 ,bisection_error= 1.637090e-011
Iteration n= 40 , c= -1.105299e+000 , f(c)= 6.977974e-012 ,bisection_error= 8.185452e-012
Iteration n= 41 , c= -1.105299e+000 , f(c)= -3.856471e-011 ,bisection_error= 4.092726e-012
Iteration n= 42 , c= -1.105299e+000 , f(c)= -1.579359e-011 ,bisection_error= 2.046363e-012
Iteration n= 43 , c= -1.105299e+000 , f(c)= -4.408029e-012 ,bisection_error= 1.023182e-012
Iteration n= 44 , c= -1.105299e+000 , f(c)= 1.285194e-012 ,bisection_error= 5.115908e-013
Iteration n= 45 , c= -1.105299e+000 , f(c)= -1.561418e-012 ,bisection_error= 2.557954e-013
Iteration n= 46 , c= -1.105299e+000 , f(c)= -1.381117e-013 ,bisection_error= 1.278977e-013
Iteration n= 47 , c= -1.105299e+000 , f(c)= 5.737633e-013 ,bisection_error= 6.394885e-014
Iteration n= 48 , c= -1.105299e+000 , f(c)= 2.176037e-013 ,bisection_error= 3.197442e-014
Iteration n= 49 , c= -1.105299e+000 , f(c)= 3.996803e-014 ,bisection_error= 1.598721e-014
Iteration n= 50 , c= -1.105299e+000 , f(c)= -4.929390e-014 ,bisection_error= 7.993606e-015
Converged within the tolerance 1.000000e-014 . The zero found is -1.105299e+000
======================================================================================

Summary of Results:

• From class we learned that the convergence of the secant method should be super-linear,  meaning |e_n+1| <= C|e_n|^alpha, where alpha is between 1 and 2 (in the case of the secant method, the value for alpha comes to ~ 1.62). I used the analysis technique shown on the project assignment webpage to produce the error graph. The downward rolloff clearly indicates that the secant method convergence is at least super-linear. 3 examples were analyzed, all showing the same trending of the convergence rate. The final example contrasts the secant method convergence rate with that of the bisection method. This shows that the secant method converegnce is faster than that of the bisection method, as is expected from theory. In fact, the bisection plot shows an approximate straight line shape, indicating linear convergence, in perfect agreement with what we learned in class.

• Another interesting point that came out of this study was how hard it was to find proper examples for this study. It was relatively easy to find oscillatory cases that did not converge, and convergent cases that converged within 15 or so samples worst case, but finding convergent cases with more than 15 samples proved to be difficult, at least using any of the 4 functions I chose as the basis of this project.

•  A final interesting trait of the dual secant/bisection co-convergence rate plot is the apparent indication that the bisection method seems to converge more quickly than the secant method during the first part of the graph. The implication is that in special cases where the starting points for the bisection are chosen judiciusly, the graph seems to indicate that it is possible for the bisection method to converge more quickly than the secant method using the same starting points. This result was borne out in a few cases while I was playing with runs of the 2 methods attempting to find good candidate function/starting point conditions to use.

I was intrigued by the discussion of the 'best' method to use to determine when the secant method has converged. I decided to explore this more fully even though this was strictly speaking not part of the actual assignment. This section contrasts the 2 other secant.m methods against the one that was used in the main portion of the report, to see the effect of changing the convergence criteria between the 3 choices mentioned in class.

Example using text criteria on top of page 124:

secant_text_conv('f_pg125',-1,1,100,1e-8)
Input Parameters
-------------------------------------------------------------------------------------
           file = f_pg125.m
       function = f(x) = x^5 + x^3 + 3
            x_0 = -1.000000e+000
            x_1 = +1.000000e+000
 Max Iterations =  100
      Tolerance =  1.000000e-008

Start Values (adjusted for polarity if required):
-------------------------------------------------------------------------------------
   x_0 = -1.000000e+000
   x_1 = +1.000000e+000
f(x_0) = +1.000000e+000
f(x_1) = +5.000000e+000
 

------------------------------ Function Output --------------------------------------
  Initial x0    x = -1.000000e+000    f(x_0) = +1.000000e+000
  Initial x1    x = +1.000000e+000    f(x_1) = +5.000000e+000
Iteration  1    x = -1.500000e+000    f(x_0) = -7.968750e+000
Iteration  2    x = -1.055749e+000    f(x_0) = +5.116487e-001
Iteration  3    x = -1.114158e+000    f(x_0) = -9.991207e-002
Iteration  4    x = -1.104615e+000    f(x_0) = +7.593217e-003
Iteration  5    x = -1.105289e+000    f(x_0) = +1.011787e-004
Iteration  6    x = -1.105299e+000    f(x_0) = -1.045169e-007
Iteration  7    x = -1.105299e+000    f(x_0) = +1.436185e-012
Iteration  8    x = -1.105299e+000    f(x_0) = +4.440892e-016

----------------------- Halt Due To Denominator < sqr(eps)! --------------------------------
Final attempt at root = -1.105299e+000
         Final f(x_0) = +4.440892e-016
      f(x_1) - f(x_0) = +1.435740e-012
 

Re-run example using text criteria on top of page 124, divide-by-zero temporarily commented out

EDU>> secant_text_conv('f_pg125',-1,1,100,1e-8)
Input Parameters
-------------------------------------------------------------------------------------
           file = f_pg125.m
       function = f(x) = x^5 + x^3 + 3
            x_0 = -1.000000e+000
            x_1 = +1.000000e+000
 Max Iterations =  100
      Tolerance =  1.000000e-008

Start Values (adjusted for polarity if required):
-------------------------------------------------------------------------------------
   x_0 = -1.000000e+000
   x_1 = +1.000000e+000
f(x_0) = +1.000000e+000
f(x_1) = +5.000000e+000
 

------------------------------ Function Output --------------------------------------
  Initial x0    x = -1.000000e+000    f(x_0) = +1.000000e+000
  Initial x1    x = +1.000000e+000    f(x_1) = +5.000000e+000
Iteration  1    x = -1.500000e+000    f(x_0) = -7.968750e+000
Iteration  2    x = -1.055749e+000    f(x_0) = +5.116487e-001
Iteration  3    x = -1.114158e+000    f(x_0) = -9.991207e-002
Iteration  4    x = -1.104615e+000    f(x_0) = +7.593217e-003
Iteration  5    x = -1.105289e+000    f(x_0) = +1.011787e-004
Iteration  6    x = -1.105299e+000    f(x_0) = -1.045169e-007
Iteration  7    x = -1.105299e+000    f(x_0) = +1.436185e-012
Iteration  8    x = -1.105299e+000    f(x_0) = +4.440892e-016
Iteration  9    x = -1.105299e+000    f(x_0) = +4.440892e-016

---------------------- Halt Due To Convergence Limit! -------------------------------

    root = -1.105299e+000    f(x_0) = 4.440892e-016
 

Example using algorithm criteria on bottom of page 124:

secant_algo_conv('f_pg125',-1,1,100,1e-8)
Input Parameters
-------------------------------------------------------------------------------------
           file = f_pg125.m
       function = f(x) = x^5 + x^3 + 3
            x_0 = -1.000000e+000
            x_1 = +1.000000e+000
 Max Iterations =  100
      Tolerance =  1.000000e-008

Start Values (adjusted for polarity if required):
-------------------------------------------------------------------------------------
   x_0 = -1.000000e+000
   x_1 = +1.000000e+000
f(x_0) = +1.000000e+000
f(x_1) = +5.000000e+000
 

------------------------------ Function Output --------------------------------------
  Initial x0    x = -1.000000e+000    f(x_0) = +1.000000e+000
  Initial x1    x = +1.000000e+000    f(x_1) = +5.000000e+000
Iteration  1    x = -1.500000e+000    f(x_0) = -7.968750e+000
Iteration  2    x = -1.055749e+000    f(x_0) = +5.116487e-001
Iteration  3    x = -1.114158e+000    f(x_0) = -9.991207e-002
Iteration  4    x = -1.104615e+000    f(x_0) = +7.593217e-003
Iteration  5    x = -1.105289e+000    f(x_0) = +1.011787e-004
Iteration  6    x = -1.105299e+000    f(x_0) = -1.045169e-007

---------------------- Halt Due To Convergence Limit! -------------------------------

    root = -1.105299e+000    f(x_0) = -1.045169e-007

Summary of Results:
The results from these two runs are interesting.  I ended up using the citeria of | f(x_n)| < tolerance, even though the algorithm method achieved results that appear to be at least as good, because this method makes more intuitive sense to me than the algorithm method. Review the test run using my selected method here.

• Text-based method (top of page 124)  |f(x_n - f(x_n-1)| < (tolerance * f(x_n)) . This method always exited due to either divide-by-zero criteria in a convergent case, or due to iteration. In short, the exit criteria was never met.  Possibly a coding error, or is the divide-by-zero test redundant in this case? The method ends up doing calculations using a value as small as (2  * eps), an unsafe practice in my opinion.  I re-ran this routine with the divide-by-zero commented out, it now exited at convergence but seems to be not a safe method, running that close to machine eps, and I can't guarantee that I would not now trip a divide-by-zero case. I discarded this approach.
• Algorithm-based method (bottom of page 124)  |(x_n - x_n-1)|/|f(x_n - f(x_n-1)| * f(x_n) < tolerance (also referred to as d in the algorithm code). This approach seems to converge the slightly quicker in one or 2 cases (like this one), than the method I employed. Since this value is not exactly equal to the next root (the method shows that the next xn = current x_n - |(x_n - x_n-1)|/|f(x_n - f(x_n-1)| * f(x_n), I am a bit confused as to what this intermediate value really is. I prefer using the simpler and easier to comprehend | f(x_n)| < tolerance criteria instead. In any case, there are no differences in the calculated root, and I only ever perceived a single iteration difference in speed of convergence between the book algorithm method and my own.

Back to top of page

In this section I descibe the files that were used in the project.
 
 

File	File Type	My Modification Level	Description of What the File Is
matlab m files
secant.m	m file	created using template supplied	This is the main file assigned to be produced. It implements the secant method using the algorithm in the textbook with modifications to the convergence criteria, small re-ordering changes, and the addition of a division-by-zero test. This file may also be run standalone, to produce a text output and a snapshot of the final function showing all the points plotted.  The syntax is secant(function_mfile_prefix,x_0,x_1,iterations,tolerance)
secant_modified.m	m file	copied	This is copied from secant.m, and then changed to comment out the sign check and order swap on the inputs to the secant method on each iteration, to discover if I could get my results to more closely match the  Web Solver results.
bisection.m	m file	supplied	This is a supplied file that runs the bisection method. It was used to help contrast the linear versus super-linear convergence rates of the bisection and secant methods. Changes were made to this code to make it ccompatible with the interface used by the  run_me_bisection file.
secant_text_conv.m	m file	copied	This is file uses the convergence criteria shown at the top of page 124. It is used in the exploratory section.
secant_algo_conv.m	m file	copied	This is file uses the convergence criteria shown at the bottom of page 124, used in the algorithm given in the text for the secant method. It is used in the exploratory section.
get_convergence_movie.m	m file	supplied	This is the file that produces the movie output, and was supplied by the professor. Minor changes were implemented in this file to make it compatible with using dynamic memory allocation, instead of the original approach. Additionally, portions of the plotting code were copied in order to be used elsewhere
run_me.m	m file	supplied	This is the top-level driver file, used to call the secant method, and create the movie from the secant.m output. Changs were made to use dynamic memory allocation,  to provide a flag to enable/disable movie creation during the call, to allow the function m file and the other secant method parameters to be specified on the command line, and the call to the bisection code was moved to a simliar routine copied from this one.
run_me_modified.m	m file	copied	This is copied from run_me.m, and then changed to call the secant_modified.m file (see above)  to discover if I could get my results to more closely match the  Web Solver results.
run_me_bisection.m	m file	copied	This is the companion file to run_me.m, used to call the bisection routine instead of the secant routine, using a similar calling structure
Function Files f_pg101.m f_pg105.m f_pg125.m f_prevtest.m	m files	created	These are the function files, called by the secant or bisection methods. The file name prefix is passed into the run_me file, and used to define which function is being evaluated.
snapshot.m	m file	created	Used to produce a plot of the convergence sequence at completion. Corresponds to the final frame that would be produced by the get_convergence_movie function. Some of the supplied code from get_convergence_movie was shamelessly copied to implement this routine...
get_func_def.m	m file	created	This file is used to extract the function definition from the funvtion m file being used so that it can be added to the plot snapshot that is created by the call to the snapshot function.
Image Files
Various	png files	created	These are the files created by running the matlab functions above. In some cases, extra command line processing was used, as in the case of the convergence plots (the semilogy approach shown on the project web page was run on the output vector after the secant and/or bisection methods were run). Files beginning with ver_ correspond to the verification section, and files beginning with anal_ correspond to the analysis section. The rest of the file name attempts to be somewhat self-explanatory, showing the function being modeled or analyzed, and something about what the test is about.
Movie Files
ver_pg105_movie.avi	avi file	created	This was created by running run_me('pg105',4,1,100,1e-8,1) to call the my secant method to evaluate the page 105 function using interval (1,4), 100 samples, a toerance of 1e-8 and passing the movie creation flag in as 'true'. This was the most intersteing movie from the verification section, in my opinion. NOTE: - The link to this file runs the movie on my windows system, but may not work on the linux system if a suitable avi file viewer is not properly installed and configured.
HTML Files
project1_main.html	html file	created	The main project page, ties the project together and has links to the other 2 pages.
project1_report.html	hrml file	created	This file. Contains the main body of the assigned work for project 1.
project1_feedback.html	html file	created	A short informal feedback form covering my likes and dislikes from this asignment.  Not assigned, but thought it might prove useful. and/or interesting...

back to top of page