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 Odd Semester 2015 15 Batch 1st Semester M.E. in Structural Engineering Subject: Computer Applications in Structural Engineering From 24/03/2015 to 12/07/2015 LECTURE 04 Conducted on: 07/04/2015 Timings: 4.30pm to 6.00pm
Computer Applications in Structural Engineering
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 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.
single root, x 0 . But, the function f x tanx 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 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 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 3 positive roots. Similarly 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 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 qx 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
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.
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
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 m1 p 0 but f m p 0 .
Thus, for
1 x f x 2 x ln 1 x
x 1 f x 21 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 cos3x;
1,1.5
Look at its graph and count number of zeros:
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.
Rate of Convergence Suppose we have a sequence of real numbers that converge to some point x:
x0 , x1, x2 , xn 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
xn1 x xn x
Generally, if lim
n
2
, then convergence is quadratic.
xn1 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
n
xn 1 0 xn 0
lim
n
xn 1 xn
lim
n
1 n 1k 1 nk
k
n . lim n n 1
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.
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 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:
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 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:
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 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.
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…………… 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.
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.
Certain singularities (poles) of f(x) can be mistaken for roots. For example, 1 f x tanx changes sign at x n , n 1,3,5, . Look at this: 2 3 At x 0.01 , tanx 99.996666644444531 2 3 At x 0.01 , tanx - 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, vt , can be expressed as: vt 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. 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:
9
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.
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
1.1288 0 , using the Incremental h
Search Method with x0 1.0 and x 0.1 .
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
x 1.400000000000000 1.410000000000000 1.420000000000000 1.430000000000000 1.440000000000000 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.
1.450000000000000 1.460000000000000 1.470000000000000 1.480000000000000 1.490000000000000 1.500000000000000
0.015735714188432 0.011695080879372 0.007698434952706 0.003745031094040 -0.000165858936611 -0.004034946800313
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:
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:
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.
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 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 . Solution: From
s 1 sinh , we have L 1100 1 sinh 1000 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.
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
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
L , we get 2 0
L 77 10 3 1000 0 2.93909385 0000000e + 07 2 20.763401 Thus, the maximum tensile stress in the cable is:
max 0 cosh 2.93909385 0000000e + 07 cosh 0.763401 =3.837926757626018e+07.
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.
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: yx
0 120 EIL
x
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
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: 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.
Thus, approximate root at 15th iteration is 268.32819. But
dy -2.5270092041 01563 e+04 dx x268.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; 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.
>> 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
ba 1
2
ba 2
3
2nd iteration will be
nth iteration will be
ba 22 ba 2n
.
Error can be bounded for each iterate mid n for the bisection method as: mid n x
ba 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 ba 2n
s
ba s
2n
ba ln s
ln 2n
ba ln s
n ln 2
ba ln s ln 2
n 17
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.
ba 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 . Convergence From the above discussion, it is found that if initial interval size b a is denoted by I 0 , then after first bisection, the interval size becomes I1
I0 I I , after second bisection I 2 1 0 2 4 2
and the story is generalized as after nth bisection,
In
I n1 I 0 n 2 2
So if the tolerance level is s , then number of bisections required to reduce the interval width to s is given by
ba ln n s ln 2
Or in other words, if we think of the current interval size as the error, the error in the nth step is proportional to the error in the (n-1)th step. This is characteristic of a linearly convergent algorithm. We say that a method is linearly convergent if the error in the nth step is proportional to the error in the (n-1)th step raised to the first power, quadratically convergent if the error in the nth step is proportional to the error in the (n-1)th step raised to the second power, quartically convergent if raised to the fourth power, and so forth. 18
