[100] 005 Chapter 5 Root and Minimum Search, Tutorial by www.msharpmath.com
[100] 005 revised on 2012.10.30 / cemmathThe Simple is the Best
Chapter 5 Root and Minimum Search
5-1 Roots of a Single Nonlinear Equation
5-2 Collection of Roots and Complex Root
5-3 Coupled Nonlinear Equations
5-4 Pseudo-Linear System
5-5 Minimum Search
5-6 Summary
The term ‘roots’ are used as a terminology for a set of solution to a given system of equations, either single or coupled. Unfortunately, there exist no cure-all methods to solve all kinds of nonlinear equations. In this regard, experience and gifted talent of intuition could be the most powerful approach in solving nonlinear equations. Although no guarantee of success is provided, root-finding methods available in the literature are discussed in this chapter. Also, several methods to find minimum of a function (without constraints though) are addressed.
Section 5-1 Roots of a Single Nonlinear Equation
In this section, we introduce the Umbrella ‘solve’ to find roots of a nonlinear equation of a single-variable. The default method of finding roots is based on the secant method which is a modified version of the Newton method. To aid users’ choice, several well-known methods are also provided with appropriate Spokes.
■ Umbrella ‘solve’. Our mission is to find a root x which satisfies the following equation
fx=g(x)
where fx and g(x) are functions of a single-variable x When g(x)≡0, the root is frequently called the zero of the equation fx=0.
Various methods available in the literature are implemented through the Umbrella ‘solve’ and Spokes listed below
solve .x ( <opt>, f(x) = g(x) ) .bisect( a,b ) // bisection
solve .x ( <opt>, f(x) = g(x) ) .secant( a,b ) // secant method
solve .x ( <opt>, f(x) = g(x) ) .muller( a,b,c ) // Muller method
solve .x ( <opt>, f(x) = g(x) ) .march ( a,b,dx ) // marching method
solve .x ( <opt>, f(x) = g(x) ) .span[n=100,g=1](a,b) // bisections
solve .x {a} ( <opt>, f(x) = g(x) ) .newton ( f'(x)-g'(x) ) // Newton
solve .x {a} ( <opt>, f(x) = g(x) ) .multiple( f'(x)-g'(x) ) // multiple
These methods will be discussed briefly together with the performance test. To avoid any confusion in using syntax, both the commands in the following
solve .x { c } ( f(x) = g(x) ) .bisect( a,b ) // bisection method
solve .x ( f(x) = g(x) ) .bisect( a,b ) // bisection method
are exactly the same, since the Hub field { c } is eliminated by Spoke ‘bisect’.
Other Spokes which control the solution procedures are
.peep show all the iteration procedures
.peeptail show the final stage of the iteration
.return( … ) return manipulated data
.maxiter(200) maximum number of iteration
.tol/abstol(1.e-7) absolute tolerance. tol is equivalent to abstol
.reltol(1.e-7) relative tolerance
where constants inside the parenthesis are the default values, and slash means that either Spoke can be alternatively used.
■ Default Method. If not specified explicitly, the default method to find a root is the secant method near x=a. The second point x=b is selected
b=a+h, h=aϵ if h<ϵ, h=ϵ
where ϵ=10-3. Therefore, the syntax for the default method reads
solve .x { a } ( <opt>, f(x) = g(x) ) // a=0.41 by default
which is simple enough for repeated use.
(Example 5-1-1) Let us find a root of the following equation
ex-3=cosx
The motive of expression-look-alike approach leads us to write
#> solve.x { 0 } ( exp(x)-3 = cos(x) );
which results in
ans = 1.2098915 ◀
(Example 5-1-2) However, for an equation ex-3=sinx, we find that
#> solve.x { 0 } ( exp(x)-3 = sin(x) );
yields the following warning
The secant method fails due to near-zero slope
ans = 1.#INF000
This is because the derivative of ex-3-sinx at x=0 vanishes to zero. Therefore, we have to select another starting point, e.g.
#> solve.x { 0.5 } ( exp(x)-3 = sin(x) );
ans = 1.3818344
We find successfully the root x=1.3818344. As in this example, it is crucial to select a starting point for the secant method. ◀
■ Spoke ‘peep’, ‘return’ and ‘tol’. Not only the root but also the procedure of root-finding may be of interest. In this case, Spoke ‘peep’ enables to watch the iteration procedure. Also, Spoke ‘tol’ controls the termination of iteration. For example
#> solve.x { 0.5 } ( exp(x)-3 = sin(x) ) .tol(0.01).peep;
shows the iteration procedure and the final root as follows.
------secant method ------
iter x f(x)
0 0.5 -1.831
0 0.501 -1.83
1 2.8707546 14.38
2 0.76768751 -1.54
3 0.97105291 -1.185
4 1.6498706 1.209
5 1.3069639 -0.2705
6 1.3696323 -0.04593
7 1.3824521 0.002345
ans = 1.3824521
It can be seen that the total number of iterations is 7 until the absolute tolerance of 0.01 is met. Also, an error 1.3824521-1.3818344=0.0006177 has occurred.
An interesting feature of Spoke ‘return’ is that it can manipulate the final return data, even changing from scalar to matrix. If one argument is returned, the results is a scalar, otherwise a matrix. This can be confirmed by the following commands
#> solve.x { 0.5 } ( exp(x)-3 = sin(x) ) .return( x/100 );
#> solve.x { 0.5 } ( exp(x)-3 = sin(x) ) .return( x+1, x, sin(x) );
ans = 0.0138183
ans =
[ 2.3818 ]
[ 1.3818 ]
[ 0.9822 ]
Such an interesting performance of Spoke ‘return’ will be judiciously used to solve coupled quations, as can be seen later.
■ Bisection Method. For an interval of a≤x≤b it is evident that the equation fx=0 has at least one root inside the interval, if fafb<0 This is the premise of the bisection method. For the same problem discussed above, let us try to apply the bisection method over an interval 1≤x≤2 by
#> solve.x ( exp(x)-3 = sin(x) ) .bisect( 1,2 ).peeptail;
which results in
------bisection method ------
iter x f(x) xa xb
23 1.3818344 7.724e-008 ( 1.3818343 1.3818345 )
ans = 1.3818344
A total of 23 iterations are necessary to find the root x=1.3818344.
■ Newton Method. The Newton method utilizes the derivative of function to find a root of equation ex-3=sinx
xn+1=xn-fnfn'
where fn=f(xn) and fn'=f'(xn). The performance of the Newton method is examined by
#> solve.x { 2 } ( exp(x)-3-sin(x) = 0 ) .newton( exp(x)-cos(x) ).peep;
which yields
------Newton method ------
iter x f(x)
0 2 3.48
1 1.5541745 0.7313
2 1.3990555 0.06608
3 1.3820259 0.0007269
4 1.3818344 9.105e-008
ans = 1.3818344
Since the convergence rate of the Newton method is of second-order, only 4 iterations are required as above.
■ Secant Method. The secant method selects two points, x1 and x2. Then, the slope f2-f1/(x2-x1) plays a similar role of the tangent f'(x) in the Newton method, i.e.
xn+1=xn-fnxn-xn-1fn-fn-1
is repeated until convergence. When f1<|f2|, the initial direction is reversed, i.e. x1 and x2 are interchanged. The secant method with two points x1=1, x2=2 is tested by
#> solve.x ( exp(x)-3-sin(x) = 0 ) .secant( 1,2 ).peep;
------secant method ------
iter x f(x)
0 1 -1.123
0 2 3.48
1 1.2440152 -0.4776
2 1.4245115 0.1665
3 1.377849 -0.01508
4 1.3817247 -0.0004161
5 1.3818347 1.087e-006
6 1.3818344 -7.796e-011
ans = 1.3818344
■ Muller Method. The Muller method employs three different points x=x1,x2,x3, and determines a parabola passing through three points
x1,f1, x2,f2, (x3,f3)
Then, two roots of the parabola are examined to discard one of them, and this procedure is repeated. Similarly, the performance test is done by
#> solve.x ( exp(x)-3 = sin(x) ) .muller( 0,1,2 ).peep;
------Muller method ------
iter x f(x)
0 0 -2
0 1 -1.123
0 2 3.48
1 1.3340585 -0.1757
2 1.3797657 -0.007839
3 1.3818408 2.445e-005
4 1.3818344 -4.388e-010
ans = 1.3818344
Note that the convergence rate is similar to that of the Newton method.
■ Marching Method. The marching method starts from x=a and marches by an increment dx until x=b. During this process, the bisection method is applied whenever sign change is noticed. Otherwise, not-a-number (NaN) is returned. The performance test is as follows
#> solve.x ( exp(x)-3 = sin(x) ) .march( 1,2, 0.1 ).peep;
------marching method ------
step x f(x)
0 1 -1.123
1 1.1 -0.887
2 1.2 -0.6119
3 1.3 -0.2943
4 1.4 0.06975
ans = 1.3818343
Note that sign change has occurred over an interval 1.3≤x≤1.4, and thus the bisection method is applied to find a root.
■ Multiple-Roots Method. The multiple-roots method is applied to the equation of type
Fx=x-cngx=0
Since Fx/F'x=x-chx has a sinlge root, the derivative is delivered through the Spoke ‘multiple’. Subsequently, the secant method is applied to find a single root of the modified equation. Consider an equation which has multiple roots
Fx=x-34=x4-12x3+54x2-108x+81=0
It should be noted that the bisection method fails to find a root since Fx≥0 always. Using Spoke ‘multiple’ enables to solve this equation near x=2
#> solve.x { 2 } ( 81 -108*x +54*x*x -12*x^3 +x^4 = 0 )
.multiple( -108 +108*x -36*x*x +4*x^3 ).peep;
------Multiple method ------
iter x f(x)
0 2 -0.25
1 3 -2.842e-014
ans = 3.0000000
Surprisingly, only one iteration is sufficient to find a root of even multiplicity.
Section 5-2 Collection of Roots and Complex Root
■ Collection of Roots. For an equation fx=0, there can exist more than one root inside an interval a≤x≤b (e.g. sinx=0.5, 0≤x≤5π). Under this circumstance, it is desired to find them all, if possible. Though not deliberate, the interval can be divided into many subintervals and the bisection method can be applied to each subinterval. This approach may require huge amount of computational cost, but serves as a convenient tool.
Collection of roots over an interval of a≤x≤b is possible by Spoke ‘span’ which implies spanning subintervals over an interval of a≤x≤b and applying the bisection method. The syntax is
solve .x ( <opt>, f(x) = g(x) ) .span [n=100,g=1] (a,b)
It should be remembered that the return value is a matrix in this case, while all other methods return by default a scalar. Under favorable conditions, it is straightforward to locate collection of roots. For example, the following equation
sinx=0.5, 0≤x≤5π
can be solved and plotted by
#> solve.x ( sin(x) = 0.5 ) .span(0,5*pi) /pi*6;
#> plot.x(0,5*pi) ( sin(x), 0.5 );
The plot is shown in Figure 1, and the returned matrix is
ans =
[ 1 5.0055e-008 ]
[ 5 -1.0602e-016 ]
[ 13 -6.0066e-007 ]
[ 17 -6.0066e-007 ]
[ 25 1.0602e-015 ]
[ 29 -1.2013e-006 ]
Since the roots were divided by π/6, the real roots can be listed as
x=π6[1,5,13,17,25,29]
and agree with the exact solutions.
Figure 1 Plot of more than one root
■ Discontinuities. When a given function includes discontinuities, Spoke ‘span’ returns n×4 matrix the last two columns of which designate the points of discontinuities. This can be confirmed by
#> solve.x ( tan(x) = x ) . span (0,5*pi) /pi*2;
ans =
[ 0 0 1 -7.012e+006 ]
[ 2.8606 -3.731e-006 3 3.506e+006 ]
[ 4.918 1.7415e-005 5 -1.0518e+006 ]
[ 6.9418 -3.3665e-005 7 5.259e+006 ]
[ 8.9548 -5.7542e-005 9 2.6295e+006 ]
Note that the first column is indeed representing the roots (but divided by π/2). The second column is the corresponding residues of function. Large values in the fourth column designate that the third column are discontinuities. In the above, a total of 5 points
xdiscontinuties=π2[ 1,3,5,7,9 ]
indeed turn out to be discontinuities, as expected. The actual roots are
xroot=π2[ 0, 2.8606, 4.918, 6.9418, 8.9548 ]
■ Difficulty in Finding Roots. Let us consider a curve including the Bessel function
fx=J03x-2sin4x2
over an interval of 0≤x≤5. The plot is shown in Figure 2 by the following command
#> plot.x[1001](0,5) ( .J_0(3*x)-2*sin(4*x)^2, 0 );
Figure 2 Function with ill-behaved roots
Although a total of 9 roots exist inside the interval 0≤x≤5 it would be readily understandable how difficult it is to find roots of fx=0 from Figure 2. Two consecutive roots are too close to be identified near x=4. The best way of finding such ill-behaved roots is of course plotting the curve and examines the possible range of roots.
In the below, the difference between two grid resolutions can be obviously perceived. The former locates 7 roots, while the latter locates 9 roots including ill-behaved roots near x=4.
#> solve.x ( .J_0(3*x)-2*sin(4*x)^2 ) . span[101](0,5);
ans =
[ 0.18672 7.0457e-008 ]
[ 0.72046 -4.8847e-008 ]
[ 0.79823 -8.1143e-008 ]
[ 2.2583 -8.5721e-007 ]
[ 2.4526 -9.7138e-007 ]
[ 4.6356 -1.2173e-006 ]
[ 4.7744 1.2927e-006 ]
#> solve.x ( .J_0(3*x)-2*sin(4*x)^2 ) . span[1001](0,5);
ans =
[ 0.18672 -5.5199e-008 ]
[ 0.72046 3.8437e-008 ]
[ 0.79823 9.7982e-009 ]
[ 2.2583 2.6107e-007 ]
[ 2.4526 -2.9943e-007 ]
[ 3.9314 -1.2532e-008 ]
[ 3.9444 3.0044e-008 ]
[ 4.6356 -8.9246e-007 ]
[ 4.7744 -3.621e-007 ]
■ Complex Root. It is also possible to find a complex root by using Umbrella ‘solve’. Since no inequality is defined for complex numbers, only the secant method will be adopted in finding complex roots. A proper syntax is
solve .z { complex number } ( <opt>, f(z) = g(z) ) // default only
This means that, once a complex number is specified inside the Hub field, Spokes describing root-finding methods cannot be used.
Consider an equation in terms of a complex variable
fz=ez+1=0, z=x+iy
It is evident that fz=0 has no real root since ez>0 for all z. Therefore, let us find a complex root starting from z=1+10i by the following commands
#> solve.z { 1+10! } ( exp(z)+1 = 0 ).peep / pi;
The results include the iteration procedure due to the use of Spoke ‘peep’ as follows.
------secant method for complex root ------
iter real imag |f(z)|
0 ( 1 10 ) 1.956
0 ( 1.011 10.03 ) 2.008
1 ( 0.65476666 9.9056557 ) 1.136
2 ( 0.42414361 9.8163938 ) 0.7145
…
25 ( 3.5603844e-008 9.4247781 ) 1.403e-007
26 ( 1.7801913e-008 9.424778 ) 7.013e-008
ans = 5.66653e-009 + 3!
In the above, note that the returned complex number has been divided by π. Thus the final root is found to be