Prepared by Asif Ali Shaikh/Sania Qureshi Assistant Professors Department of BSRS Mehran, UET, Jamshoro

Note: Detailed discussion of following topics is being carried out in the classes. It is highly recommended that you may follow textbooks suggested by your mentors. These handouts can never ever replace textbooks.

MEHRAN UNIVERSITY OF ENGINEERING AND TECHNOLOGY, JAMSHORO. FIRST SEMESTER FIRST YEAR (M.Phil.) 15BATCH PROGRAM: APPLIED MATHEMATICS Subject: Scientific Computing LECTURE 03 Conducted on: 01/04/2015 Timings: 3pm to 4.30pm

NUMERICAL SOLUTION OF ALGEBRAIC AND TRANSCENDENTAL EQUATIONS Objectives: o Polynomials and Transcendental Functions o Root and Zero o Simple Root and Root with Multiplicity o Descarte’s Rule of Signs o Mean Value Theorem o Rate of Convergence or Asymptotic Error Constant o Order of Convergence o Incremental Search Method o Graphical Method o Bisection Method In this document, our objective is to find root of a polynomial (preferred to be of some high degree) or to find zero of a transcendental function of the form f x   0 and to analyze iterative techniques exist to serve our purpose. Discussion for Rate of Convergence and computational complexities associated with such techniques will be our priority which will be carried out by comparing various techniques to each other. Based upon details for function of the form f x   0 , it may have only one root, several roots or infinite number of roots and can have no root at all. Moreover, roots can be real or complex. For example, the function f x  sin x  x hase a single root, x  0 . But, the function f x  tanx  x has an infinite number of roots x  0,4.493,7.725 Motivation: Given fig(a) shows a uniform beam subject to a linearly increasing distributed load. The equation for the resulting elastic curve

1

Prepared by Asif Ali Shaikh/Sania Qureshi Assistant Professors Department of BSRS Mehran, UET, Jamshoro

Note: Detailed discussion of following topics is being carried out in the classes. It is highly recommended that you may follow textbooks suggested by your mentors. These handouts can never ever replace textbooks.

is (see fig(b)):

y x  

 x 120 EIL 0

5

 2L2 x3  L4 x



Use an iterative technique to determine the point of maximum deflection (i.e., the value of x where dy dx  0 ). Then substitute this value into the given model to determine the value of the

maximum

deflection.

Use

the

following

parameters

in

your

computation:

L  600 cm, E  50,000 kN / cm , I  30,000 cm , and 0  2.5kN / cm . (Hint: dy dx  0 2

4

changes the sign between 267 and 269). This problem will later be solved by one of the iterative techniques developed for f x   0 . Before to dive into the sea of such techniques, let’s discuss some of the important terms. Polynomials Polynomials are special types of nonlinear algebraic equation of the general form:

f n  x   an x n  an 1 x n 1  an  2 x n 2    a1 x  a0 ; where n  order Coefficients a' s can be real and/or complex. In general, an nth order polynomial will have n roots. Polynomials are extremely important in various fields of science and engineering such as Curve Fitting. Descarte’s Rule of Signs The equation f x   0 cannot have more positive roots than the changes of signs in f x  and cannot have more negative roots than the changes of signs in f  x  . For example, f x   7 x8  4 x 5  9 x  12

 ve to  ve  ve

to  ve

 ve to  ve Thus, f x   0 cannot have more than 3 positive roots. Similarly f  x   7 x8  4 x 5  9 x  12

 ve to  ve  ve to  ve

 ve to  ve 2

Prepared by Asif Ali Shaikh/Sania Qureshi Assistant Professors Department of BSRS Mehran, UET, Jamshoro

Note: Detailed discussion of following topics is being carried out in the classes. It is highly recommended that you may follow textbooks suggested by your mentors. These handouts can never ever replace textbooks.

Thus, f x   0 cannot have more than 1 negative root. Generally: If an equation f x   0 of nth degree has at the most   p positive roots and q negative roots then the equation has at least   :

n   p  q complex roots. The equation MATLAB as:

f x   x 4  4 x 3  7 x 2  6 x  3 has only complex roots computed using

-1.5000 + 0.8660i -0.5000 + 0.8660i

-1.5000 - 0.8660i -0.5000 - 0.8660i

The equation f x   x 4  8 x3  14 x 2  8 x  15 has 3 positive roots and 1 negative root computed using MATLAB as: 5.0000

3.0000

-1.0000

1.0000

Root with Multiplicity: A root p of the equation f x   0 is said to be a root of Multiplicity m if f x  can be written in the form: f  x    x  p  q  x  where lim qx   0 m

x p

Note that a root with multiplicity 1 is called a simple root. Example: f x   x 6  x 5  12 x 4  2 x 3  41x 2  51x  18 The given polynomial is written as f  x    x  1  x  3  x  2  3

2

Thus, it has multiplicity 3 at x  1 , multiplicity 2 at x  3 and x  2 is a simple root. Note: Using MATLAB, factorization is performed by the command factor(polynomial) and if factors are given then the command expand(factors) can be used. Example: Non-polynomial (transcendental) equation: 1 x  f  x   2 x  ln   1 x  3

Prepared by Asif Ali Shaikh/Sania Qureshi Assistant Professors Department of BSRS Mehran, UET, Jamshoro

Note: Detailed discussion of following topics is being carried out in the classes. It is highly recommended that you may follow textbooks suggested by your mentors. These handouts can never ever replace textbooks.

Clearly, x  0 is the zero of the function but what about multiplicity? Theorem: Let f x  be a continuous function with m continuous derivatives. The equation

f x   0 has a root of multiplicity m at x  p if and only if

f  p   f  p   f  p   f  p     f m1  p   0 but f m   p   0 . Thus, for

1 x  f  x   2 x  ln   1 x   1  x     f  x   21    f x   4 and  2  1 x2  1 x  



  2  



f 0  f 0  f 0  0 but f 0  0 Hence, f x   0 has multiplicity 3 at x  0 . Intermediate Value Theorem If f x   0 is continuous (we can draw it without lifting our pencil from the paper) in a closed interval a, b and f a  f b  0 , then f x   0 will have at least one real root within the closed interval. Note: If there is more than 1 roots in a, b then they will always be odd in numbers. Consider: f x  sin 10x  cos3x;

1,1.5

Look at its graph and count number of zeros:

Rate of Convergence Suppose we have a sequence of real numbers that converge to some point x:

x0 , x1, x2 ,  xn 4

Prepared by Asif Ali Shaikh/Sania Qureshi Assistant Professors Department of BSRS Mehran, UET, Jamshoro

Note: Detailed discussion of following topics is being carried out in the classes. It is highly recommended that you may follow textbooks suggested by your mentors. These handouts can never ever replace textbooks.

Is there a way to talk about how fast the numbers are converging to x? Definition: If

xn 1  x

lim

n 

xn  x

 

then the sequence converges linearly to x. The constant   0 is called the asymptotic error or rate of convergence. Similarly, if lim

n

xn1  x xn  x

Generally, if lim

n

2

    , then convergence is quadratic.

xn1  x xn  x



    , then order of convergence of the sequence is determined

to be  .

1 for some fixed k  0 . Then we know that this sequence converges to nk 0. Its order of convergence is determined as: Example: Let xn 

lim

xn 1  0

n 

xn  0



 lim

n 

xn 1 xn



 lim

n 

1 n  1k 1 nk

k

 n   .  lim  n  n  1   



Thus, it will be convergent if  is set to be 1 (linear convergence). There are several possible ways to solve nonlinear mathematical model of the form f x   0 .  Analytical Solutions  possible for special equations only  Graphical Illustration  Useful for providing initial guesses for other methods  Numerical Solutions  Open methods  Bracketing methods

5

Prepared by Asif Ali Shaikh/Sania Qureshi Assistant Professors Department of BSRS Mehran, UET, Jamshoro

Note: Detailed discussion of following topics is being carried out in the classes. It is highly recommended that you may follow textbooks suggested by your mentors. These handouts can never ever replace textbooks.

We will start our discussion of Root Location Numerical Methods from Graphical Method which is better to get an initial insight about possibility for zero of given nonlinear mathematical model. It also gives us an initial guess to start an iterative technique. Graphical Method: Step 01. Write the given mathematical model in the form f x  g x . Step 02. Plot both of these functions on the same graphing window. Step 03. Locate approximately their point of intersection using data cursor available in MATLAB. This point of intersection is the approximate solution of the given mathematical model. It is to be noted that graphical method is in no way a recommended method for solving practical problems. However, it is helpful to get an initial guess for iterative techniques to be discussed soon. Example. Use graphical method to find a real root of x3  x 1  0 ? The given equation can be written as: x3  x  1 . Here, f x   x 3 and g x   x  1. Now plot both of these functions on MATLAB as given below:

Thus, the approximate root is found to be 1.325. Example. Use graphical method to find a real root of cos x  3x  1  0 ? The given equation can be written as: cos x  3x  1. Here, f x  cos x and g x  3x  1 . Now plot both of these functions on MATLAB as given below:

6

Prepared by Asif Ali Shaikh/Sania Qureshi Assistant Professors Department of BSRS Mehran, UET, Jamshoro

Note: Detailed discussion of following topics is being carried out in the classes. It is highly recommended that you may follow textbooks suggested by your mentors. These handouts can never ever replace textbooks.

Thus, the approximate root is found to be x = 0.6062 at which y = 0.0032. Order of Convergence The number of iterations required to reach the root depends largely on the intrinsic order of convergence of the method. Letting Ek be the error in the computed root after the kth iteration, an approximation of the error after the next iteration has the form Ek 1  cEkm

c 1

The order of convergence is determined by m. If m = 1, then the method is said to converge linearly. Methods for which m > 1 are called super-linearly convergent. The best methods display quadratic convergence (m = 2). There are various iterative schemes to find solution of nonlinear models with one independent variable and one scheme is given preference over the other for many reasons discussed in the class. Having learnt all of these schemes, it would be extremely useful for you to do a comparative study about the schemes from different aspects which will certainly increase your knowledge and sense of comprehension. Some of the commonly used techniques for f x   0 are mentioned below:             

Incremental Search Method Bisection or Interval-Halving Method Regula-Falsi Method Ridder’s Method Successive Approximation Method Newton-Raphson Method Generalized Newton-Raphson Method Modified Newton-Raphson Method Secant Method Muller’s Method Horner’s Method Graeffe’s Root-Squaring Method Brent’s Method

 Aitken’s 2 Method  Many more…………… 7

Prepared by Asif Ali Shaikh/Sania Qureshi Assistant Professors Department of BSRS Mehran, UET, Jamshoro

Note: Detailed discussion of following topics is being carried out in the classes. It is highly recommended that you may follow textbooks suggested by your mentors. These handouts can never ever replace textbooks.

INCREMENTAL SEARCH METHOD All methods for finding roots are iterative procedures that require a starting point, that is, an estimate of the root. This estimate can be crucial; a bad starting value may fail to converge, or it may converge to the “wrong” root (a root different from the one sought). There is no universal recipe for estimating the value of a root. If the equation is associated with a physical problem, then the context of the problem (physical insight) might suggest the approximate location of the root. Otherwise, the function must be plotted (a rough plot is often sufficient), or a systematic numerical search for the roots can be carried out. It is highly advisable to go a step further and bracket the root (determine its lower and upper bounds) before passing the problem to a root finding algorithm. The approximate locations of the roots are best determined by plotting the function. Often a very rough plot, based on a few points, is sufficient to give us reasonable starting values. Another useful tool for detecting and bracketing roots is the incremental search method. It can also be adapted for computing roots, but the effort would not be worthwhile, since other methods described in rest of the document and discussed in the class are more efficient for that. The basic idea behind the incremental search method is simple: if f (x1) and f (x2) have opposite signs, then there is at least one root in the interval (x1, x2). If the interval is small enough, it is likely to contain a single root. Thus the zeros of f (x) can be detected by evaluating the function at intervals x and looking for change in sign. There are a couple of potential problems with the incremental search method:  It is possible to miss two closely spaced roots if the search increment x is larger than the spacing of the roots.  A double root (two roots that coincide) will not be detected.  Certain singularities (poles) of f(x) can be mistaken for roots. For example, 1 f x  tanx changes sign at x   n , n  1,3,5,  . Look at this: 2 3 At x    0.01 , tanx  99.996666644444531 2 3 At x    0.01 , tanx  - 99.996666644448197 2 Example: Water is discharged from a reservoir through a long pipe as shown. By neglecting the change in the level of the reservoir, the transient velocity of the water flowing from the pipe, vt  , can be expressed as: vt   t   tanh  2 gh , 2 gh  2L 

where h is the height of the fluid in the reservoir, L is the length of the pipe, g is the acceleration due to gravity, and t is the time elapsed from the beginning of the flow.

8

Prepared by Asif Ali Shaikh/Sania Qureshi Assistant Professors Department of BSRS Mehran, UET, Jamshoro

Note: Detailed discussion of following topics is being carried out in the classes. It is highly recommended that you may follow textbooks suggested by your mentors. These handouts can never ever replace textbooks.

Find the value of h necessary for achieving a velocity of v = 5m/s at time t = 3s when L = 5m and g = 9.81m/s2. Solution: Substituting the value of given parameters into the mathematical model, we get



1.1288  tanh 1.3288 h h



By drawing both sides on same graphing window, we have got a rough estimate of the solution, that is; h  1.45 .

Let’s apply Incremental Search Method on the above problem but it should be remembered that the method is the most basic automated numerical method to solve nonlinear equations. Procedure of the method is described below:  Pick a starting point x0 and a step size x . Use a positive x if you want to search to the right, and a negative x if you want to search to the left.  Let x1  x0  x and calculate f x0  and f x1  .

 If the sign of f x  changes between x0 and x1 , it is assumed that a root of f x  exists in the interval .  If the sign of f x  does not change between x0 and x1 , let x2  x1  x and repeat the procedure.





Example: Find root of the equation tanh 1.3288 h  Search Method with x0  1.0 and x  0.1 .

9

1.1288  0 , using the Incremental h

Prepared by Asif Ali Shaikh/Sania Qureshi Assistant Professors Department of BSRS Mehran, UET, Jamshoro

Note: Detailed discussion of following topics is being carried out in the classes. It is highly recommended that you may follow textbooks suggested by your mentors. These handouts can never ever replace textbooks.

Solution: We will keep evaluating the given equation until it changes sign as given below:

f x 

x

x0  1.0 x1  1.1 x2  1.2

x3  1.3 x4  1.4

x5  1.5

f x0   0.25984425 9672842

f x1   0.192295658863906

f x2   0.133648679802203

f x3   0.08218929 7694172

f x4   0.036626039042872

f x5   -0.0040349 46800313

Observation for the above table tells us that there is a root between 1.4 and 1.5 as the sign of the function is changed. Repeating the above procedure with x0  1.4 and x  0.01 would produce more accurate estimate of the root.

f x  0.036626039042872 0.032353654449872 0.028129258592595 0.023952012684799 0.019821097782423 0.015735714188432 0.011695080879372 0.007698434952706 0.003745031094040 -0.000165858936611 -0.004034946800313

x 1.400000000000000 1.410000000000000 1.420000000000000 1.430000000000000 1.440000000000000 1.450000000000000 1.460000000000000 1.470000000000000 1.480000000000000 1.490000000000000 1.500000000000000

It is to be depicted from above table that the required root lies between 1.48 and 1.49. Repeating the procedure with x0  1.48 and x  0.001 would give us even more accurate root. BISECTION METHOD This method is iterative in nature and also known as Interval Halving or Binary Chopping or Bolzano’s Method. The method comes under the category of Closed (Bracketing) Methods. The function f x   0 whose zero is required must be continuous on a, b and f a  f b  0 . It depends on the idea of halving the interval a, b within which root is assumed to be present. After that, function is evaluated at the midpoint and the new interval is decided using the concept of Intermediate Value Theorem. The iterations will remain continue until stopping criterion is met. The entire process is summarized below:

10

Prepared by Asif Ali Shaikh/Sania Qureshi Assistant Professors Department of BSRS Mehran, UET, Jamshoro

Note: Detailed discussion of following topics is being carried out in the classes. It is highly recommended that you may follow textbooks suggested by your mentors. These handouts can never ever replace textbooks.

Smaller the interval containing the root; fewer will be the iterations to approximate the root. Also note that Bisection Method is certain to converge to the root once the root is bracketed but the method is slow convergent as compared to its other counterparts discuss later. Graphical representation of the method is shown below:

Example: A steel cable of length s is suspended as shown. The maximum tensile stress in the cable, which occurs at the supports, is

 max   0 cosh  where



L 2 0

 0  tensile stress in the cable at O.

  weight of the cable per unit volume L  horizontal span of the cable 11

Prepared by Asif Ali Shaikh/Sania Qureshi Assistant Professors Department of BSRS Mehran, UET, Jamshoro

Note: Detailed discussion of following topics is being carried out in the classes. It is highly recommended that you may follow textbooks suggested by your mentors. These handouts can never ever replace textbooks.

The length to span ratio of the cable is related to  by s 1  sinh  L 

Find if   77  103 N / m3 steel  , L  1000 m and s  1100 m . s 1  sinh  , we have L 

Solution: From

1100 1  sinh  1000 

1



sinh   1.1

sinh   1.1 sinh   1.1  0

This last equation is a nonlinear model that we will solve using bisection method. A rough estimate of the root using graphical method is obtained as 0.5, 1 .

f 0.5  -0.028904694506253 and f 1  0.075201193643801 First Iteration:

mid 

0.5  1  0.75 and f 0.75  -0.002683268064170 2

0.75 0.5 negative negative Thus, the root lies between 0.75 and 1.

1 positive

Second Iteration:

mid 

0.75  1  0.875 and f 0.875  0.028506637144295 2

0.75 0.875 negative positive Thus, the root lies between 0.75 and 0.875.

 per 

1 positive

0.875  0.75 100 %  14.2857142 85714285% 0.875 12

Prepared by Asif Ali Shaikh/Sania Qureshi Assistant Professors Department of BSRS Mehran, UET, Jamshoro

Note: Detailed discussion of following topics is being carried out in the classes. It is highly recommended that you may follow textbooks suggested by your mentors. These handouts can never ever replace textbooks.

Third Iteration:

mid 

0.75  0875  0.8125 and f 0.8125  0.011143738566064 2

0.75 0.8125 negative positive Thus, the root lies between 0.75 and 0.8125.

 per 

0.875 positive

0.8125  0.875  100 %  7.69230769 2307693% 0.8125

Iterations and percentile errors are shown in the following table:

13

Prepared by Asif Ali Shaikh/Sania Qureshi Assistant Professors Department of BSRS Mehran, UET, Jamshoro

Note: Detailed discussion of following topics is being carried out in the classes. It is highly recommended that you may follow textbooks suggested by your mentors. These handouts can never ever replace textbooks.

Thus, approximate value of  is 0.763401 which can be visualized from the graph given below. From the equation  

0 

L , we get 2 0

L 77 10 3 1000    2.93909385 0000000e + 07 2 20.763401 

Thus, the maximum tensile stress in the cable is:

 max   0 cosh   2.93909385 0000000e + 07 cosh 0.763401  =3.837926757626018e+07.

Finally, you must be analyzing last column of the above table that ratio of errors, that is, next error divided by the previous is about 1/2. This is asymptotic error or rate of convergence of the Bisection Method. Also note that order of convergence of the method is 1. Thus, it is linearly convergent method. This method is a slow one requiring large number of iterations to get a good estimate of the root. In addition to this, the method fails for roots with multiplicity m. Try the following:

f x   x 6  10 x5  40 x 4  82 x3  91x 2  52 x  12 Motivation. Let’s solve the problem given in start of this lecture:

yx  

 x 120 EIL 0

5

 2L2 x3  L4 x



Substituting the given parameters’ values, we get



yx   2.3148 e-14  x 5  720000 x 3  1.2960 e+11x

14



Prepared by Asif Ali Shaikh/Sania Qureshi Assistant Professors Department of BSRS Mehran, UET, Jamshoro

Note: Detailed discussion of following topics is being carried out in the classes. It is highly recommended that you may follow textbooks suggested by your mentors. These handouts can never ever replace textbooks.





dy  2.3148 e-14  5x 4  2160000 x 2  1.2960 e+11 dx





dy  0  2.3148 e-14  5x 4  2160000 x 2  1.2960 e+11  0 dx 5x4  2160000x2  1.2960e+11  0 f 267  1.026367605000000e+09 and f 269  -5.191883950000000e+08 Function changes the sign between 267 and 269. Look at the following table calculated using bisection method for absolute error tolerance  0.0001:

Thus, approximate root at 15th iteration is 268.32819. But 15

Prepared by Asif Ali Shaikh/Sania Qureshi Assistant Professors Department of BSRS Mehran, UET, Jamshoro

Note: Detailed discussion of following topics is being carried out in the classes. It is highly recommended that you may follow textbooks suggested by your mentors. These handouts can never ever replace textbooks.

dy  -2.5270092041 01563 e+04 dx x268.32819 More number of iterations are required to get a good estimate. MATLAB calculates: >> fzero('5*x^4-2160000*x^2+1.2960e+11',[267 269]) ans = 2.683281572999748e+02 >> x=2.683281572999748e+02; >> y=5*x.^4-2160000*x.^2+1.2960e+11 y= -4.577636718750000e-05 This is much better and acceptable. Hence, the point of maximum deflection is approximately 2.683281572999748e+02. Number of Iterations required in Bisection Method: Now we know that in bisection method, after each iteration the interval is halved. Size of the interval after 1st iteration will be



3rd iteration will be



ba 1

2

ba 2

3

2nd iteration will be



nth iteration will be



ba 22 ba 2n

.

Error can be bounded for each iterate mid n for the bisection method as: mid n  x 

ba 2n

Note that as n   , mid n  x . If you have to keep continue your iterations until size of given interval   s , then it means ba 2n

 s

16

Prepared by Asif Ali Shaikh/Sania Qureshi Assistant Professors Department of BSRS Mehran, UET, Jamshoro

Note: Detailed discussion of following topics is being carried out in the classes. It is highly recommended that you may follow textbooks suggested by your mentors. These handouts can never ever replace textbooks.

 ba     s

   2n  

 

 ba ln   s

   ln 2n  

 ba ln   s

   n ln 2   

 ba ln   s ln 2

  n

 ba ln   Or n   s ln 2

  

You can now easily determine approximate number of iterations using the information of interval and pre-specified error tolerance in the question. For example, how many iterations are required by bisection method for f x   3x  sin x  e x on 0,1 if  s  0.0001 ?

 1 0  ln   0.0001   or n  13.2877 n ln 2 Therefore, approximately 14 iterations are enough in bisection method in order to achieve value of the root containing error  0.0001 .

17