CHAPTER 4: MATHEMATICAL MODELING WITH MATLAB

** **

**Lecture 4.4: Root finding
algorithms for nonlinear equations**

__ __

**Problem:** Given a nonlinear equation*:*

*f(x) = 0*

Find
all zeros (roots) ** x = x_{*}** of the function

**Example:
**Consider a
tunnel diode that has the following voltage ** (V)** - current

*I = I(V) = V ^{3} -
1.5 V^{2} + 0.6 V *

The
tunnel diode is connected with a resistor ** (R)** and a voltage
source

For given values of ** E**
and

% simple
example of root finding algorithms for polynomial functions

R = 19.2; E = 1.4; c = [ 1,-1.5,0.6+1/R,-E/R];

% coefficients of the
polynomial function f(x) = 0

p = roots(c)' % three steady currents throguh the tunnel
diode

E = 1.8; c = [ 1,-1.5,0.6+1/R,-E/R]; p = roots(c)' % one (real) steady current

p = 0.8801 0.3099 - 0.1023i 0.3099 + 0.1023i

**Methods
of solution: **

·
Bracketing
(bisection) methods

·
Iterative
(contraction mapping, Newton-Raphson, and secant) methods

·
MATLAB
root finding algorithms

__Iterative
(contraction mapping) method:__

** **

Construct
the following mapping from ** x_{k }**to

*x _{k+1 }= g(x_{k})
= x_{k }+ f(x_{k})*

Starting
with an initial approximation ** x_{0},** use the mapping
above to find

* *

% Graphical
solution of the nonlinear equation: x*sin(1/x) = 0.2*exp(-x)

x = 0.0001 : 0.0001 : 2; f1 = x.*sin(1./x); f2 = 0.2*exp(-x); plot(x,f1,'b',x,f2,'g');

% Contraction
mapping method for finding the root of the nonlinear equation:

% f(x) = x*sin(1/x)-0.2*exp(-x)

x = 0.364; % starting approximation for the root

eps = 1; tol = 0.0001; total = 8; k = 0;

while ((eps > tol) & (k < total)) % two termination
criteria

xx = x + x*sin(1/x)-0.2*exp(-x); % a new approximation for the
root

eps = abs(xx-x); x = xx; k = k+1;

fprintf('k = %2.0f, x = %6.4f\n',k,x);

end

stability = abs(1 + sin(1/x)-cos(1/x)/x+0.2*exp(-x));

if (stability >= 1)

fprintf('The contraction
mapping method is UNSTABLE');

else

fprintf('The contraction
mapping method is STABLE\n');

fprintf('The numerical
approximation for the root: x = %6.4f',x);

end

k = 2, x = 0.3684

k = 3, x = 0.3827

k = 4, x = 0.4391

k = 5, x = 0.6442

k = 6, x = 1.1832

k = 7, x = 2.0070

k = 8, x = 2.9393

The contraction mapping method is UNSTABLE

The
main issue of the iterative method is to ensure that the sequence ** {x_{k}}
**really converges to a fixed point

The
contraction mapping method is stable under the condition that

*|g'(x _{*})|
= |1 + f'(x_{*})| < 1*

When
the method is stable, the absolute error ** e_{k }= | x_{k }–
x_{*} | **reduces with the number of iterations

** e_{k+1}
= c e_{k}, **where

x = 0.4; % starting approximation for the root

eps = 1; tol = 0.000001; total = 10000; k = 0;

while ((eps > tol) & (k < total)) % two termination
criteria

xx = x + x^3-1.5*x^2+0.6*x-E/R+x/R; % a new approximation for the
root

eps = abs(xx-x); x = xx; k = k+1; e(k) = x-0.880133;

end

stability = abs(1 + 3*x^2-3*x+0.6+1/R);

if (stability >= 1)

fprintf('The contraction
mapping method is UNSTABLE');

else

fprintf('The contraction
mapping method is STABLE\n');

fprintf('The numerical
approximation for the root: x = %8.6f',x);

end

ee = e(2:k)./e(1:k-1); plot(ee);

The contraction mapping method is STABLE

The numerical approximation for the root: x = 0.532249

__Iterative
(Newton-Raphson) method:__

** **

Construct
the following mapping from ** x_{k }**to

* *

*x _{k+1 }= g(x_{k})
= x_{k }– _{}*

* *

Starting
with an initial approximation ** x_{0},** use the mapping
above to find

*y =
f(x _{k}) + f'(x_{k}) (x-x_{k})*

The
estimate of the root ** x = x_{k+1} **is obtained by setting:

* *

The
Newton-Raphson method is always stable. It converges when ** x_{0}**
is sufficiently close to

%
Newton-Raphson method for finding the root of the nonlinear equation:

% f(x) = x*sin(1/x)-0.2*exp(-x)

x = 0.55; % starting approximation for the root

eps = 1; tol = 10^(-14); total = 100; k = 0; format long;

while ((eps > tol) & (k < total)) % two termination
criteria

f = x*sin(1/x)-0.2*exp(-x);
% the function at x = x_k

f1 =
sin(1/x)-cos(1/x)/x+0.2*exp(-x); % the derivative at x = x_k

xx = x-f/f1; % a new approximation for the root

eps = abs(xx-x); x = xx; k = k+1;

fprintf('k = %2.0f, x = %12.10f\n',k,x);

end

k = 2, x = 0.3717802027

k = 3, x = 0.3636237820

k = 4, x = 0.3637156975

k = 5, x = 0.3637157087

k = 6, x = 0.3637157087

% Example of
several roots of the nonlinear equations:

% Newton-Raphson method converges to one of the roots,

% depending on the initial approximation

% f(x) = cos(x)*cosh(x) + 1 = 0 (eigenvalues of a uniform beam)

eps = 1; tol = 10^(-14); total = 100; k = 0; x = 18;

while ((eps > tol) & (k < total))

f = cos(x)*cosh(x)+1; f1 =
sinh(x)*cos(x)-cosh(x)*sin(x);

xx = x-f/f1; eps = abs(xx-x); x = xx; k = k+1;

end

fprintf('k = %2.0f, x = %12.10f\n',k,x);

k = 7, x = 17.2787595321

while ((eps > tol) & (k < total))

f = cos(x)*cosh(x)+1; f1 = sinh(x)*cos(x)-cosh(x)*sin(x);

xx = x-f/f1; eps = abs(xx-x); x = xx; k = k+1;

end

fprintf('k = %2.0f, x = %12.10f\n',k,x);

while ((eps > tol) & (k < total))

f = cos(x)*cosh(x)+1; f1 =
sinh(x)*cos(x)-cosh(x)*sin(x);

xx = x-f/f1; eps = abs(xx-x); x = xx; k = k+1;

end

fprintf('k = %2.0f, x = %12.10f\n',k,x);

while ((eps > tol) & (k < total))

f = cos(x)*cosh(x)+1; f1 =
sinh(x)*cos(x)-cosh(x)*sin(x);

xx = x-f/f1; eps = abs(xx-x); x = xx; k = k+1;

end

fprintf('k = %2.0f, x = %12.10f\n',k,x);

__Iterative
(secant) method:__

** **

The
slope of the curve ** y = f(x)** at the point

* *

*f'(x _{k}) _{}_{}*

With
the backward approximation, the Newton-Raphson method becomes the secant
method:

*x _{k+1 }= g(x_{k})
= x_{k }– _{}*

* *

Starting
with two initial approximations ** x_{0},x_{1},** use
the mapping above to find

* *

In
order to start the secant method, two starting points ** x_{0}**
and

% Secant
method for finding the root of the nonlinear equation:

% f(x) = x*sin(1/x)-0.2*exp(-x)

eps = 1; tol = 10^(-14); total = 100; k = 1; format long;

x = 0.55; f = x*sin(1/x)-0.2*exp(-x); % starting approximation for
the root

x0 = x; f0 = f; x = x + 0.01;

while ((eps > tol) & (k < total)) % two termination
criteria

f = x*sin(1/x)-0.2*exp(-x);
% the function at x = x_k

xx = x-f*(x-x0)/(f-f0); % a
new approximation for the root

eps = abs(xx-x); x0 = x; f0 = f; x = xx; k = k+1;

fprintf('k = %2.0f, x = %12.10f\n',k,x);

end

k = 3, x = 0.3878094689

k = 4, x = 0.3650168333

k = 5, x = 0.3636699429

k = 6, x = 0.3637157877

k = 7, x = 0.3637157087

k = 8, x = 0.3637157087

k = 9, x = 0.3637157087

__MATLAB
root finding algorithms:__

__ __

·
**fzero**: searches for a zero of a nonlinear function
starting with an initial approximation

·
**fminbnd: **searches for a minimum of a nonlinear function in a finite interval

·
**fminsearch:** searches for a minimum of a nonlinear function with an initial
approximation (using the Nelder-Mead simplex method of direct search)

% MATLAB
algorithn for finding the root of the nonlinear equation:

% f(x) = x*sin(1/x)-0.2*exp(-x)

x0 = 0.55; f = 'x*sin(1/x)-0.2*exp(-x)';

x = fzero(f,x0);

fprintf('The root is found at x = %12.10f\n',x);

The root is found at x = 0.3637157087

__ __

% Search of zero of f(x) is the same as the
search of minimum for [f(x)]^2

x = fminbnd('(x*sin(1/x)-0.2*exp(-x))^2',0.2,0.5)

x = fminsearch('(x*sin(1/x)-0.2*exp(-x))^2',0.3)

opt = optimset('TolX',10^(-10)); % termination tolerance on the
value of x

x = fminbnd('(x*sin(1/x)-0.2*exp(-x))^2',0.2,0.5,opt)

x = 0.36369140625000

x = 0.36371570742172

% The search
fails when no roots are found

% The function y = 1 + x^2 has no zeros

x = fzero('1 + x^2',1)

Exiting fzero: aborting search for an interval containing a sign
change

because NaN or Inf
function value encountered during search

(Function value at -1.716199e+154 is Inf)

Check function or try again with a different starting value.

x = NaN

__ __