BSCIT-403 Computer Oriented Numerical Methods.pdf

___________ BAOU Dr. Babasaheb Ambedkar Educa on for All Open University (Established by Government of Gujarat) COMPUTER ORIENTED...

___________ BAOU Dr. Babasaheb Ambedkar Educa on for All Open University (Established by Government of Gujarat) COMPUTER ORIENTED NUMERICAL METHOD BSCIT-403 Bachelor of Science (Hons.) SEM-4 Informa on Technology 2022 Computer Oriented Numerical Methods Dr. Babasaheb Ambedkar Open University Computer Oriented Numerical Methods Expert Committee Prof. (Dr.) Nilesh K. Modi Professor and Director, School of Computer Science, (Chairman) Dr. Babasaheb Ambedkar Open University, Ahmedabad Prof. (Dr.) Ajay Parikh Professor and Head, Department of Computer Science (Member) Gujarat Vidyapith, Ahmedabad Prof. (Dr.) Satyen Parikh Dean, School of Computer Science and Application (Member) Ganpat University, Kherva, Mahesana M. T. Savaliya Associate Professor and Head (Member) Computer Engineering Department Vishwakarma Engineering College, Ahmedabad Mr. Nilesh Bokhani Assistant Professor, School of Computer Science, (Member) Dr. Babasaheb Ambedkar Open University, Ahmedabad Dr. Himanshu Patel Assistant Professor, School of Computer Science, (Member Secretary) Dr. Babasaheb Ambedkar Open University, Ahmedabad Course Writer Dr. Ramesh Kataria Assistant Professor, Som-Lalit Institute of Computer Application - Ahmedabad Content Editors Prof. (Dr.) Nilesh K. Modi Professor and Director, School of Computer Science, Dr. Babasaheb Ambedkar Open University, Ahmedabad Mr. Nilesh Bokhani Assistant Professor, School of Computer Science, Dr. Babasaheb Ambedkar Open University, Ahmedabad ISBN – Printed and published by: Dr. Babasaheb Ambedkar Open University, Ahmedabad While all efforts have been made by editors to check accuracy of the content, the representation of facts, principles, descriptions and methods are that of the respective module writers. Views expressed in the publication are that of the authors, and do not necessarily reflect the views of Dr. Babasaheb Ambedkar Open University. All products and services mentioned are owned by their respective copyrights holders, and mere presentation in the publication does not mean endorsement by Dr. Babasaheb Ambedkar Open University. Every effort has been made to acknowledge and attribute all sources of information used in preparation of this learning material. Readers are requested to kindly notify missing attribution, if any.  Dr. Babasaheb BSCIT-403 Ambedkar Open University Computer Oriented Numerical Methods BLOCK-1: COMPUTER ARITHMETIC AND SOLVING NON-LINEAR EQUATIONS UNIT-1 COMPUTER ARITHMETIC: INTRODUCTION 04 UNIT-2 ERRORS 12 UNIT-3 SOLVING NON-LINEAR EQUATIONS 17 BLOCK-2: SOLVING SIMULTANEOUS LINEAR ALGEBRAIC EQUATIONS AND INTERPOLATION UNIT-1 MATRICES 39 UNIT-2 DETERMINANTS 47 UNIT-3 SOLVING SIMULTANEOUS LINEAR ALGEBRAIC EQUATIONS 56 UNIT-4 INTERPOLATION 82 UNIT-5 SPLINE INTERPOLATION 99 BLOCK-3: CURVE FITTING UNIT-1 METHOD OF LEAST SQUARE 108 UNIT-2 CURVE FITTING 118 BLOCK-4: NUMERICAL DIFFERENTIATION AND INTEGRATION UNIT-1 NUMERICAL DIFFERENTIATION 131 UNIT-2 NUMERICAL INTEGRATION 135 UNIT-3 NUMERICAL SOLUTION OF ORDINARY DIFFERENTIAL EQUATIONS 142 UNIT-4 NUMERICAL SOLUTION OF HIGHER ORDER DIFFERENTIAL EQUATIONS 150  Computer Oriented Numerical Method BLOCK1: Computer Arithmetic and Solving Non-Linear Equations UNIT 1 COMPUTER ARITHMETIC 2 UNIT 2 ERRORS 14 UNIT 3 SOLVING NON-LINEAR EQUATIONS 26 1 BLOCK 1: COMPUTER ARITHMETIC AND SOLVING NON-LINEAR EQUATIONS Block Introduction Science and Engineering often require numerical results involving all types of numbers, its arithmetic. Numbers and numerical arithmetic or computation are the subjects of algebra, a branch of Mathematics. There are four basic arithmetic operations: addition, subtraction, multiplication, and division. Computer arithmetic involves the representation of integers and real numbers in computer systems, as well as the manipulation of those numbers through hardware circuits or software programs. In this block, we will introduce different types of the number presentation of numbers, arithmetic operations on it. Concept of errors and numerical methods to solve the non-linear equations. Unit:2 focuses on floating-point presentation of numbers and due to limitation of the systems errors are introduced in the numerical calculations, are also discussed here. In Unit: 3, different types of equations and some numerical methods are discussed to obtain the solution/s of the non-linear equations. Block Objective This block aims to make students aware about arithmetic operations in computer systems, different presentations of the numbers, storage of real numbers and errors in the numerical calculation. Finally, the block will clear the concept of the roots/solutions of an equation, graphical meaning of the root/solution of equations, different types of equations and various iterative methods to solve non-linear equations. 2 Block Structure BLOCK1: UNIT1 Computer Arithmetic Objectives, Mathematical Background, Number presentation, Floating point Arithmetic, Let us Sum Up UNIT2 Errors Objectives, Errors, Different type’s errors. Let us Sum Up UNIT3 Solving Non-Linear Equations Objectives, Introduction, Types of non-linear equations, Methods of solving non-linear equations, Convergence of iterative methods, let us Sum Up 3 UNIT 1 COMPUTER ARITHMETIC: INTRODUCTION Unit Structure 1.0 Learning Objectives 1.1 Number Systems 1.2 Number Representation 1.2.1 Floating Point Arithmetic 1.2.2 Addition Operation 1.2.3 Subtraction Operation 1.2.4 Multiplication Operation 1.2.5 Division Operation 1.3 Limitation of Floating-Point Representation 1.4 Let Us Sum Up 1.5 Suggested Answer for Check Your Progress 1.6 Glossary 1.7 Assignment 1.8 Activities 1.9 Further Readings 4 1.0 LEARNING OBJECTIVES: After learning this unit, you should be able to: Understand the floating-point arithmetic Understand different arithmetic operations on floating-point numbers Learn errors due to arithmetic operations 1.1 NUMBER SYSTEMS: A number system is defined as a system of writing to express numbers. It is the mathematical notation for representing numbers of a given set by using digits or other symbols in a consistent manner. It provides a unique representation of every number and represents the arithmetic and algebraic structure of the figures. It also allows us to operate arithmetic operations like addition, subtraction and division. The value of any digit in a number can be determined by: The digit Its position in the number The base of the number system 1.2 NUMBER REPRESENTATION: There are many potential sources of error in mathematical computations. To try to acquire an application of how errors arise, it is first useful to observe how number are represented in computers. Now, let us discuss how many different categories are there for different types of Data-types in details: 1.2.1 Floating point arithmetic Fractional Quantities are generally represented in computers by floating point form. In floating point form a number is represented in three parts: One for Sign One for Fractional part One for Exponent Thus, the representation of a float is of the following general form: Float + sign * mantissa * bexponent For example: The Numbers (i) 0.7292 * 104 (ii) 0.1516 *10-13 5  Also written as (i) 0.7292 E 04 (ii) 0.1516 E -13 A floating-point number is said to be normalised if the first digit of the mantissa is always 1. The shifting of mantissa to the left till it’s most significant digit is non-zero is known as normalisation. Zero is a special case, because its fractional part has all zeros and a zero exponent, so zero can never be normalised. 1.2.2 Addition operation Addition operation is performed on normalized floating–point number if the exponents of the numbers is equal. If the exponents are not equal the exponent of the number with smaller exponent is made equal to larger exponent and its mantissa is modified. For example, if 0.7253E2 and 0.5467E5 are to be added, then the decimal point of mantissa of 0.7253E2 is shifted by 3 positions (i.e. 5-2) to left hand side and the exponent is increased by 3. Due to which exponent of both numbers become equal and addition operation can be operated on them by adding their mantissas. Example 1.1 Add 0.8475E6 to 0.4376E3 Solution: The decimal point of the mantissa of 0.4376E3 is shifted three position to left. It becomes 0.0004, whereas the digits 6, 7, and 3 are chopped off i.e. truncated. The exponent is incremented by 3. Therefore, the number after normalization becomes 0.0004E6. These numbers can be added now as shown below: Addend 0.0004E6 Augend 0.8475E6 Sum 0.8479E6 Example 1.2 Add 0.8475E2 to 0.4376E2 Solution: Since The exponents are already equal, they can be directly added as: Addend 0.4376E2 Augend 0.8475E2 Sum 1.2851E2 In this case, the mantissa of the sum is greater than 1.0, so the decimal point is shifted to the left by one position and the exponent is increased by 1. The last digit of mantissa of sum is truncated giving the normalized sum as 0.1285E3. 6 Example 1.3 Add 0.8475E99 to 0.4376E99 Solution: Since the exponents are already equal, they can be directly added as: Addend 0.4376E99 Augend 0.8475E99 Sum 1.2851E99 In this case, the mantissa of the sum is greater than 1.0, so the decimal point is shifted to the left by one position and the exponent is increased by 1. The last digit of the mantissa of sum is chopped off, giving the normalized sum as 0.1285E100. Since this number is greater than the largest number which our hypothetical computer can handle, it is a case of overflow and the system will indicate this condition. 1.2.3 Subtraction operation Like addition operation, the Subtraction operation is performed on normalized floating–point number if the exponents of the numbers is equal. If the exponents are not equal the exponent of the number with smaller exponent is made equal to larger exponent and its mantissa is modified. During normalization if exponent of difference becomes less than -99 than the system will show an error because the number is smaller than the smallest number which our hypothetical computers can store. And this situation is known as underflow. For example, if 0.7253E2 and 0.5467E5 are to be subtracted, then the decimal point of mantissa of 0.7253E2 is shifted by 3 positions (i.e. 5-2) to left hand side and the exponent is increased by 3. Due to which exponent of both numbers become equal and subtraction operation can be operated on them by subtracting their mantissas. Example 1.4 Subtract 0.8475E3 to 0.4376E6 Solution: The number 0.8475E3 is normalized, thus resulting in number 0.0008E6. Now the number 0.0008E6 an be subtracted from 0.4376E6 as: Minuend 0.4376E6 Subtrahend 0.0008E6 Difference 0.4368E6 Example 1.5 Subtract 0.8475E6 to 0.8476E6 Solution: Since the exponents are already equal, the number 0.8475E6 can be directly subtracted from 0.8476E6 as: Minuend 0.8476E6 Subtrahend 0.8475E6 7  Difference 0.0001E6 In this case, the mantissa of the difference is less than 0.1, so the decimal point is shifted three positions to the right and the exponent is decreased by 3. The resultant difference becomes 0.1E3. Check Your Progress-1 1. Add x1 = 264.9 to x2 = 26.00 [A] 0.2909×103 [B] 0.2654×103 [C] 0.6754×10-3 [D] 0.7459×10-3 2. Add 1.37 and 0.0269 [A] 0.139×101 [B] 0.153×102 [C] 0.757×10-1 [D] None 3. Subtract 5481 from 7482 [A] 0.2001 [B] 0.2001×103 [C] 0.2001×10-3 [D] None 1.2.4 Multiplication operation In multiplication operation, the mantissas are multiplied and the exponents are added. In case of multiplication operation, the mantissa of product will not be in normalized form. For example, if 0.1234E5 and 0.1111E13 are to be multiplied the we directly multiply the mantissa and add the exponents. i.e. 0.1234E5 * 0.1111E13 = (0.1234 * 0.1111) E (5+13) Example 1.6 Multiply 0.3754E5 by 0.8576E4 Solution: 0.3754E5 × 0.8576E4 = (0.3754 × 0.8576) E (5+4) = 0.32194304E9 Thus, we obtain 0.1475E8 after truncating the mantissa of the product to four digits. And this product is already in the normalized form. Example 1.7 Multiply 0.2345E5 by 0.1111E15 Solution: 0.2345E5 ×0.1111E15 = (0.2345 × 0.1111)E(5+15) =0.02605295E20 Thus, we obtain 0.0260E20 after truncating the mantissa of the product to four digits. And since the product is less than 0.1, the decimal point is shifted one position to the right and the exponent is decreased by 1. The resultant product becomes 0.2605E19. 8  1.2.5 Division operation In Division Operation, the mantissa of first number is divided by the mantissa of the second number and the exponent of second number is subtracted from the exponent of first number. Note, that the mantissa of the quotient will not be normalized. Further Note that the magnitude of the quotient may become greater than 1.0, but will always be less than 10.0. Therefore, at most the decimal point of mantissa of the quotient will be shifted one position to the left and the exponent will be increased by 1. And as a result of the increase, the exponent can increase up to +99. Now if the mantissa is negative, this results in underflow or overflow. Example 1.8 Divide 0.9876E-5 by 0.1231E99 Solution: 0.9876E-5 ÷ 0.1231E99 = (0.9876 ÷ 0.1231)E(-5-99) = 8.02274574E-104 The mantissa of the quotient is greater than 1.0 therefore the decimal point is shifted one position to the left and the exponent is increased by 1. The resultant quotient becomes 0.8022E-103. Example 1.9 Divide 0.9876E5 by 0.1231E-99 Solution: 0.9876E5 ÷ 0.1231E-99 = (0.9876 ÷ 0.1231)E(5-(-99)) = 8.02274574E104 The mantissa of the quotient is greater than 1.0 therefore the decimal point is shifted one position to the left and the exponent is increased by 1. The resultant quotient becomes 0.8022E105. Check Your Progress-2 0.0785 0.785 10−1 1. Compute = 3580 0.385 104 [A] 0.219×10-4 [B] 0.219×104 [C] 0.253×105 [D] 0.253×10-5 2. Divide 0.10000 E 99 to 0.6457 E -6 [A] 0.6457 E -105 [B] 0.6457 E 105 [C] 0.1548 E -105 [D] 0.1548 E 105 3. If we multiply 0.1237 E 51 to 0.5286 E 55 then result obtained would [A] Overflow [B] Underflow [C] Can’t be said [D]None of the Above 9  1.3 LIMITATION OF FLOATING-POINT REPRESENTATION: Given below are major limitations of floating-point representation: Floating point representation is basically complex representation as it uses two fields: mantissa and exponent. Length of registers for storing floating point numbers is large. 1.4 LET US SUM UP In this unit, we: Discussed number systems and floating-point representation Explained floating-point arithmetic 1.5 SUGGESTED ANSWERS FOR CHECK YOUR PROGRESS Check Your Progress-1 1. [A] 0.2909×103 2. [A] 0.139×101 3. [B] 0.2001×103 Check Your Progress-2 1. [A] 0.219×10-4 2. [A] 0.6457 E -105 3. [A] Overflow 1.6 GLOSSARY 1. Operation is a function which takes zero or more input values to obtain well-defined output value. 2. Arithmetic Operation is part of mathematics that involves the study of numbers and the operation of number that are useful in all other fields of mathematics. 3. Mantissa is the part of logarithm after the decimal point. 4. Error is the difference between the true value and the approximation of that value. 10  1.7 ASSIGNMENT 1. What is normalized floating-point representation? Express the number 0.000008754 in the normalized floating-point representation. 2. (i) Divide 0.9854 E 5 by 0.1987 E 3 (ii) Multiply 0.6585 E 56 by 0.4375 E -53 (iii) Subtract 0.3846 E 4 by 0.1736 E 4 (iv) Add 0.7364 E 2 by 0.7686 E 4 1.8 ACTIVITY 1. Find the value of ( 6 × 4 ) ÷ 12 + 72 ÷ 8 – 9 2. Find the value of 108 ÷ 3 + (5 × 2) - 36 1.9 FURTHER READING M.K. Jain, S.R.K. Iyengar and R.K. Jain, Numerical Methods for Scientific and Engineering Computation, 6th Ed., New age International Publisher, India, 2007. S.S. Sastry, Introductory Method of Numerical Analysis, 5th Edition, PHI Learning Private Limited, New Delhi, 2012 11 UNIT 2: ERRORS Unit Structure 2.0 Learning Objectives 2.1 Errors 2.1.1 Errors in Computation 2.1.2 Absolute and Relative errors 2.2 Let Us Sum Up 2.3 Suggested Answer for Check Your Progress 2.4 Glossary 2.5 Assignment 2.6 Activities 2.7 Further Readings 12 2.0 LEARNING OBJECTIVES: Whenever a numerical method is applied to the problem, then the numerical result obtained is certainly approximate, i.e. it contains some errors. From this chapter students will learn and understand errors occurring in numerical computation. 2.1 ERRORS Whenever a numerical operation is applied to a problem, then the result obtained is certainly approximate value, i.e. the value obtained contains some errors. An error is defined as the difference between the actual value of the problem and the value obtained from the numerical operation. Consider x represents some quantity and xa is approximations to x. Then; Error = actual value – approximate value = x -xa 2.1.1 Errors in Computation Computational errors are mostly occurred due to following reasons: Normalized floating-point representation Truncation of infinite series expansion Inefficient algorithm As we have earlier seen in previous chapter that while performing different arithmetic operations using normalized floating-point representation, some digits have to be shorten in order to fit it in normalized float-point form. And due to such shortenings various types of errors occur depending upon size of the numbers and the type of operations used. Similarly, errors are introduced because of infinite series or process, such as sin(x), cos(x), log(x), e^x, etc. For example, the following infinite series x3 x5 x 7 sin( x) = x − + − +..... 3! 5! 7! Is used to find the value of sin(x). Since the series is infinite, it is not possible to use all the terms in computations. In such cases the process is executed till a finite number of terms and then terminated. And because of the terms omitted due to termination will create error in the computed result. 2.1.2 Truncation And Round-off errors Truncation Error: In truncation errors, some digits are excluded from the number. There are mainly two situations when i) Numbers are represented in normalized floating-point form. The mantissa can accommodate only a few digits, for example, only four in our hypothetical computer. In our hypothetical computer, 0.0356879 would take the form 0.3568E-1 in normalized floating-point representation. The digits 7 and 9 have been discarded. This truncation introduces errors in the input data. 13 ii) A number is converted from one system to another. For example, in binary, 13.1 corresponds to 1101.0001100110011... There is a repeating fraction in this number, which means the conversion must be terminated after some digits. This introduces the truncation error. Round-off Error: Round-off errors occur when floating-point numbers are chopped, rounded, or used in arithmetic operations. Computers can only represent real numbers with a fixed precision of mantissa. A computer's representation of true values may not always be accurate. This is called round-off error. It is possible that round-off errors are the most troublesome, since they can accumulate and cause significant problems. Check Your Progress-1 1. Let x = 0.00573735. Find the absolute error if x is truncated to the three decimal points. 2. The solution of a problem is given as 5.284. It is known that the absolute error in the solution is less than 0.01.Find the interval within which the exact value must lie. 2.1.3 Absolute and Relative errors Absolute Error: Absolute Error is the modulus of the difference between the actual value and the approximate value of the given problem. If x is the actual value and xa is the approximate value of the problem then the Absolute Error is given by: Absolute Error = | x – xa | Relative Error: Relative Error is the ratio of the error to the actual value of the problem. If x is the actual value and xa is the approximate value of the problem, then the relative error is given by: x − xa Relative Error = x Example 2.1 Let x = 0.00473817. Find the absolute error if x is truncated to three decimal point. Solution: Given, x = 0.00473819 = 0.473819×10-2 xa = 0.473×10-2 error = x – xa = 0.473819×10-2 - 0.473×10-2 Hence | x – xa | = 0.000473×10-2 = 0.473×10-5 -2-3 Which is less than 10. Example 2.2 14 Given the solution of a problem as x0 = 46.93 with the relative error in the solution at most 4%. Find, to four decimal digits, the range of values within which the exact value the solution must lie. Solution: Given maximum relative error = 0.04 Therefore x − x0 −0.04   0.04 x x − x0 If  0.04 x Then, x − x0  0.04 x  x(1 − 0.04)  xa x0 x (1 − 0.04) 46.93 = = 48.8854167 0.96 x − x0 However, if −0.04  x Then x − x0  −0.04 x  x(1 + 0.04)  x0 x0 x 1 + 0.04 46.93 = = 45.125 1.04 Hence, 45.125 < x < 48.8854167 Check Your Progress-2 1. Let the exact value of a number be 7.745. If it is round-off to two decimal digits, find the value of absolute and relative error and also write relative percentage. 2. Find the relative percentage error in approximate representation of 4/3 by 1.33. 2.2 LET US SUM UP In this unit, we: Learn how errors are computed and how it can be prevent Learned different types of errors and how to calculate them. 2.3 SUGGESTED ANSWERS FOR CHECK YOUR PROGRESS Check Your Progress-1 1. 0.573×10-5 15  2. (5.274, 5.294) Check Your Progress-2 1. 0.5% 2. 0.0645% 2.4 GLOSSARY 1. Truncation means to be shortened off. 2. Inefficient means not working or producing in best way possible 3. Error is the difference between the true value and the approximation of that value. 4. Modulus is the of the positive value of any problem. 2.5 ASSIGNMENT  1. Write down approximate value of correct to four significant digits and then find (a) absolute 4 error (b) relative error (c) relative percentage. 2. If  x = 0.006 and  y = 0.003 be the absolute errors in x = 2.56 and y = 5.24, find the relative error in computation of x + y. 3. Let x = 0.0083465, find the relative error if (i) x is rounded off to three decimal digits (ii) If x is chopped off to three decimal digits 2.6 ACTIVITY 1. A thermometer is calibrated 1500C to 2000C. The accuracy is specified within ±0.25%. What is the maximum static error? 2. A Calculator has a capacity to calculate numbers up to 12 digits. And the accuracy is specified within ±0.20%. What is the maximum error possible? 2.7 FURTHER READING M.K. Jain, S.R.K. Iyengar and R.K. Jain, Numerical Methods for Scientific and Engineering Computation, 6th Ed., New age international publisher, India, 2007. S.S. Sastry, Introductory Method of Numerical Analysis, 5th Edition, PHI Learning Private Limited, New Delhi, 2012 16 UNIT 3: SOLVING NON-LINEAR EQUATIONS Unit Structure 3.0 Learning Objectives 3.1 Introduction 3.2 Types of Non-Linear Equations 3.2.1 Polynomial equations 3.2.2 Transcendental equations 3.3 Methods of Solving Non-Linear Equations 3.3.1 Direct methods 3.3.2 Iterative methods 3.4 Iterative Methods 3.4.1 Bisection method 3.4.2 False position method 3.4.2 Secant method 3.4.3 Newton-Raphson method 3.5 Convergence of Iterative Methods 3.6 Let Us Sum Up 3.7 Suggested Answer for Check Your Progress 3.8 Glossary 3.9 Assignment 3.10 Activities 3.11 Further Readings 17 3.0 LEARNING OBJECTIVES In this unit we will learn variety of methods to find approximations for the roots of an equations. And we will also learn methods use to find solutions of algebraic and Transcendental Equations which are represented by f(x) =0. Iteration is the fundamental principle in computer science. As the name says iteration means that the process is repeated until an answer is achieved. This technique is used to roots of equation, solution of linear equation and nonlinear system of equations. 3.1 INTRODUCTION Finding solution of non-linear equations is of interest not only to mathematicians but also to scientist and engineers. Though most of the non-linear equation can be solved algebraically i.e. analytically. But there are many non-linear equations which cannot be solved algebraically. For example; 32x + 3x – ex = 0 It seems very simple but cannot be solved algebraically. The solution obtained by algebraic manipulation are known as algebraic solutions or analytical solutions. There are many examples in which algebraic solutions of non-linear equations does not exist but is extremely complicated and impractical for day-to-day purposes. In this chapter we will learn some methods to find solutions of non-linear equations, which will be numerical and not algebraic, and called numerical solutions. 3.2 TYPES OF NON-LINEAR EQUATIONS 3.2.1 Polynomial equations The polynomials are frequently occurring functions in science and engineering field. A polynomial equation can be written in following general form: anxn + an-1xn-1 + an-2xn-2 +.... + a2x2 + a1x + a0 = 0 where an ≠0 It is a nth degree polynomial in x and has n roots. These roots may be Real and Different Real and Equal Complex Note: Complex numbers appear in pairs and are in form of a+ib and a-ib where i=……., ia an imaginary number, and a and b are real numbers and represent real and imaginary parts of root. 18 3.2.2 Transcendental equations Non-polynomial equations are called transcendental equations. Some examples of transcendental equations are xex - xsinx = 0 excosx – 3x = 0 x e3 x − = 0 2 2 –x–3=0 x Note: A Transcendental equation may have finite/infinite number of real roots or may not have any real root at all. Check Your Progress-1 1. Find the approximate value of a real root of f(x) = xlog10x – 1.2 = 0 2. Find the approximate value of a real root of the equation sinx – x3 – 1 = 0 3.3 METHODS OF SOLVING NON-LINEAR EQUATIONS In general, there are two kind of methods to find solution of non-linear equations. These are: 1. Direct methods 2. Iterative methods 3.3.1 Direct methods Direct method gives the roots of non-linear equations in a few steps and this methos is also capable to give all the roots at the same time. For example, the root of quadratic equation: ax2 + bx + c = 0 Where a ≠ 0 are given by −b  b 2 − 4ac X 1,2 = 2a 19 3.3.2 Iterative methods Iterative methods are also known as trial-and-error methods. An iterative method involves repeating approximations in order to arrive at a result. A series of approximations are obtained by repeating a fixed sequence of steps until the solution is accurate enough. Generally, through Iterative method we only get one root at a time. And this method is very difficult and time consuming for solving the non-linear equation manually. However, it is suitable for use on computer because of following reasons: 1. Iterative methods can be concisely expressed as computational algorithms. 2. Round-off error are negligible in iterative methods as compared to direct methods. 3. It is possible to formulate algorithms which can handle class of similar problems. 3.4 ITERATIVE METHODS With this knowledge, Let’s discuss some very well-known iterative methods to find solution of the algebraic and transcendental equations. 3.4.1 Bisection method This is one of the simplest iterative methods. Bisection method is also known as Bolzano method or Interval method. To start with two initial approximations x0 and x1 such that f(x0)f(x1)< 0 which means f(x0) and f(x1) have opposite signs, which ensures that the next root lies between x0 and x1, then next root say x2, is the midpoint of the interval [x0, x1] is computed. Now, If f(x2) =0, then we get the root at x2. If f(x0) and f(x2) are of opposite sign, then the root lies between the interval (x0, x2). Thus, x1 is replaced by x2, and the new interval is formed, which is the half of the previous interval. And this interval is again bisected. If f(x1) and f(x2) are of opposite sign, then the root lies between the interval (x2, x1). Thus, x0 is replaced by x2, and the new interval is formed, which is the half of the previous interval. And this interval is again bisected. Similarly, this procedure is continued until a root of desired accuracy is obtained. 20  Figure3.1 Root approximation by approximation method Example 3.1 Find a root of an equation f(x)=2x3-2x-5 using Bisection method Solution: Here 2x3-2x-5=0 Let f(x)=2x3-2x-5 Here x 0 1 2 f(x) -5 -5 7 1st iteration: Here f(1)=-50 ∴ Now, Root lies between 1 and 2 x0 = (1+2)/2 = 1.5 f(x0) = f(1.5) = 2⋅1.53-2⋅1.5-5 = -1.250 3rd iteration : Here f(1.5) = -1.250 ∴ Now, Root lies between 1.5 and 1.75 21 x2 = (1.5+1.75)/2 = 1.625 f(x2) = f(1.625) = 2⋅1.6253-2⋅1.625-5 = 0.332>0 4th iteration : Here f(1.5) = -1.250 ∴ Now, Root lies between 1.5 and 1.625 x3 = (1.5+1.625)/2 = 1.5625 f(x3) = f(1.5625) = 2⋅1.56253-2⋅1.5625-5 = -0.49560 8th iteration : Here f(1.5938) = -0.09110 ∴ Now, Root lies between 1.5938 and 1.6016 x7 = (1.5938+1.6016)/2 = 1.5977 f(x7) = f(1.5977) = 2⋅1.59773-2⋅1.5977-5 = -0.0393 x1, then the problem is classified as a boundary value problem. In this book we will consider only the initial value problem. 3.1 EULER’S METHOD One of the simplest and oldest methods is Euler's method. Euler's method involves developing a piecewise linear approximation to the solution. For the initial value problem, the slope of the solution curve is given as well as the starting point. A solution curve is extrapolated using the specified step size based on this information. Take a look at the following first order ordinary differential equation dy = f ( x, y ) dx With the initial conditions as: y = y1 for x = x1 Using the mean value theorem we can find the solution of the above problem And, the mean value theorem states that, “If a function is continuous and differentiable between two points A(x1,y1) and B(x2,y2), then the slope of the line joining these points is equal to the derivative of the function at least at one other point c between these two points” y ( x2 ) − y ( x1 ) y '(c) = x2 − x1 Substituting c = x1 and h = x2 – x1, in the above equation we get y ( x2 ) = y ( x1 ) + hf ( x1 , y1 ) And, we know that y '( x1 ) = f ( x1 , y1 ) So,  hf ( x1 , y1 ) = y ( x2 ) − y ( x1 ) 143   y ( x2 ) = y ( x1 ) + hf ( x1 , y1 )  y3 = y2 + hf ( x2 , y2 ) Similarly, if we take (x2,y2) as the starting point,  y3 = y2 + hf ( x2 , y2 ) And we can write generalized form of the above equations as;  yi +1 = yi + hf ( xi , yi ) Example 3.1 Find y(0.5) for y′=-2x+3y, x0=0,y0=-1, with step length 0.1 using Euler method. Solution: Given y′=-2x+3y,y(0)=-1,h=0.1,y(0.5)=? Euler method y1 = y0 + hf(x0,y0) = -1+(0.1)f(0,-1) = -1+(0.1)⋅(-3) = -1+(-0.3)= -1.3 y2 = y1+hf(x1,y1) = -1.3+(0.1)f(0.1,-1.3) = -1.3+(0.1)⋅(-4.1) = -1.3+(-0.41)= -1.71 y3 = y2+hf(x2,y2) = -1.71+(0.1)f(0.2,-1.71) = -1.71+(0.1)⋅(-5.53) = -1.71+(-0.553) = -2.263 y4 = y3+hf(x3,y3)= -2.263+(0.1)f(0.3,-2.263) = -2.263+(0.1)⋅(-7.389) = -2.263+(-0.7389) = -3.0019 y5 = y4+hf(x4,y4) = -3.0019+(0.1)f(0.4,-3.0019) = -3.0019+(0.1)⋅(-9.8057) = -3.0019+(-0.9806) = - 3.9825 ∴y(0.5) = -3.9825 Example 3.2 Find y(0.7) for y′=7x-9y/3, x0=0,y0=1, with step length 0.2 using Euler method. Solution: Given y′=7x-9y3,y(0)=1,h=0.2,y(0.7)=? Euler method y1 = y0 + hf(x0,y0) = 1+(0.2)f(0,1) = 1+(0.2)⋅(-3) = 1+(-0.6) =0.4 y2 = y1+hf(x1,y1) = 0.4+(0.2)f(0.2,0.4) = 0.4+(0.2)⋅(-0.7333) = 0.4+(-0.1467) = 0.2533 y3 = y2+hf(x2,y2) = 0.2533+(0.2)f(0.4,0.2533) = 0.2533+(0.2)⋅(0.1733) = 0.2533+(0.0347) = 0.288 ∴y(0.6)=0.288 Check Your Progress-1 1. Find y(0.5) for y’ = -2x-y, y(0) = -1, with step length 0.1 x− y 2. Find y(2) for y ' = , y(0)=1, with step length 0.2 2 144 3.2 RUNGE-KUTTA SECOND ORDER METHODS Each of the Runge-Kutta second order methods matches the Taylor series method up to the second degree terms in h, where h is the step size. For these methods, the interval [xi, xf] is divided into subintervals and all derivatives (slopes) within those intervals are averaged to determine the dependent variable. In theory, like Euler's method, these methods are one-step methods, in which information at the preceding (xi, yi) point is all we need to evaluate yi+1. Consider the following differential equation dy = f ( x, y ) dx With the initial conditions as: y = y1 for x = x1 At the starting point, Let the slope be s1 = f ( x1 , y1 ). Now let the second point be x2 = x1 + h where y2 = y1 + s1h and let the slope be s2 = f ( x2 , y2 ). So, the average value of the slope would be s = (s1 +s2)/2. And then the value of the solution would be calculated using this new average value of slope as y2 = y1 +sh Similarly, general form of formula for third, fourth terms. We can write the value of y at the (i+1)th position as yi+1 = yi +sh si + si +1 where, s = 2 si = f ( xi , yi ) si +1 = f ( xi + h, yi + si h) Example 3.3 Find y(0.7) for y′=x2y, x0 = 0,y0 = -1, with step length 0.1 using Runge-Kutta 2 method Solution: Given y′=x2y , y(0) = -1, h=0.1, y(0.7)=? k1 = hf(x0,y0) = (0.1)f(0,-1) = (0.1)⋅(0) = 0 k2 = hf(x0+h,y0+k1) = (0.1)f(0.1,-1) = (0.1)⋅(-0.01) = -0.001 y1=y0+k1+k2/2 = -1-0 = -1.0005 ∴y(0.1) = -1.0005 Again taking (x1,y1) in place of (x0,y0) and repeat the process k1 = hf(x1,y1) = (0.1)f(0.1,-1.0005) = (0.1)⋅(-0.01) = -0.001 k2= hf(x2+h,y1+k1) = (0.1)f(0.2,-1.0015) = (0.1)⋅(-0.0401) = -0.004 145 y2 = y2+k2+k2/2 = -1.0005-0.0025 = -1.003 ∴y(0.2) = -1.003 Again taking (x2,y2) in place of (x1,y1) and repeat the process k1 = hf(x2,y2) = (0.1)f(0.2,-1.003) = (0.1)⋅(-0.0401) = -0.004 k2 = hf(x2+h,y2+k1) = (0.1)f(0.3,-1.007) = (0.1)⋅(-0.0906) =-0.0091 y3 = y2+k1+k2/2 = -1.003-0.0065 = -1.0095 ∴y(0.3) = -1.0095 Again taking (x3,y3) in place of (x2,y2) and repeat the process k1 = hf(x3,y3) = (0.1)f(0.3,-1.0095) = (0.1)⋅(-0.0909) = -0.0091 k2 = hf(x3+h,y3+k1) = (0.1)f(0.4,-1.0186) = (0.1)⋅(-0.163) = -0.0163 y4 = y3+k1+k2/2 = -1.0095-0.0127 = -1.0222 ∴y(0.4) = -1.0222 Again taking (x4,y4) in place of (x3,y3) and repeat the process k1 = hf(x4,y4) = (0.1)f(0.4,-1.0222) = (0.1)⋅(-0.1636) = -0.0164 k2 = hf(x4+h,y4+k1) = (0.1)f(0.5,-1.0386) = (0.1)⋅(-0.2596) = -0.026 y5 = y4+k1+k2/2 = -1.0222-0.0212 = -1.0434 ∴y(0.5) = -1.0434 Again taking (x5,y5) in place of (x4,y4) and repeat the process k1 = hf(x5,y5) = (0.1)f(0.5,-1.0434) = (0.1)⋅(-0.2608) = -0.0261 k2 = hf(x5+h,y5+k1) = (0.1)f(0.6,-1.0695) = (0.1)⋅(-0.385) = -0.0385 y6 = y5+k1+k2/2 = -1.0434-0.0323 = -1.0757 ∴y(0.6) = -1.0757 Again taking (x6,y6) in place of (x5,y5) and repeat the process k1 = hf(x6,y6) = (0.1)f(0.6,-1.0757) = (0.1)⋅(-0.3872) = -0.0387 k2 = hf(x6+h,y6+k1) = (0.1)f(0.7,-1.1144) = (0.1)⋅(-0.5461) = -0.0546 y7 = y6+k1+k2/2 = -1.0757-0.0467 = -1.1224 ∴y(0.7) = -1.1224 Example 3.4 Find y(0.9) for y′=xy, x0=0,y0=-1, with step length 0.3 using Runge-Kutta 2 method. Solution: Given y′=xy,y(0)=-1,h=0.3,y(0.9)=? k1 = hf(x0,y0) = (0.3)f(0,-1) = (0.3)⋅(0) = 0 k2 = hf(x0+h,y0+k1) = (0.3)f(0.3,-1) = (0.3)⋅(-0.3) = -0.09 y1 = y0+k1+k2/2 = -1-0.045 = -1.045 ∴y(0.3) = -1.045 Again taking (x1,y1) in place of (x0,y0) and repeat the process k1 = hf(x1,y1) = (0.3)f(0.3,-1.045) = (0.3)⋅(-0.3135) = -0.094 k2 = hf(x1+h,y1+k1) = (0.3)f(0.6,-1.139) = (0.3)⋅(-0.6834) = -0.205 y2 = y1+k1+k2/2 = -1.045-0.1495 = -1.1945 ∴y(0.6) = -1.1945 Again taking (x2,y2) in place of (x1,y1) and repeat the process k1 = hf(x2,y2) = (0.3)f(0.6,-1.1945) = (0.3)⋅(-0.7167) = -0.215 k2 = hf(x2+h,y2+k1) = (0.3)f(0.9,-1.4096) = (0.3)⋅(-1.2686) = -0.3806 146 y3 = y2+k1+k2/2 = -1.1945-0.2978 = -1.4923 ∴y(0.9) = -1.4923 Check Your Progress-2 1. Find y(0.3) for y’ = -(xy2 + y), y(0) = 1, with step length 0.1 2. Find y(0.2) for y’ = -y, y(0) = 1, With step length 0.1 3.3 RUNGE-KUTTA FOURTH ORDER METHODS Consider the following differential equation dy = f ( x, y ) dx With the initial conditions as: y = y1 for x = x1 Let the starting point be x1 at which y1 and let the slope be s1 = f(x1,y1) Let the second point be at x2 = x1 + (h/2) at which y2 = y1 +s1(h/2) and let the slope be s2 = f(x2,y2). Let the third point be at x3 = x1 + (h/2) at which y3 = y1 + s2(h/2) and let the slope be s3 = f(x3,y3). Let the third point be at x4 = x1 + h at which y4 = y1 + s3h and let the slope be s4 = f(x4,y4). Let the average of the slope be given by s = (s1 + 2s2 + 2s3 + s4)/6 Thus the value of the dependent variable y is computed as y2 = y1 + hs In general the value of y at (i+1)th point of solution curve is obtained from ith point using the following equation. yi+1 = yi +hs Example 3.5 Find y(0.8) for y′=x+5y, x0=0,y0=-1, with step length 0.2 using Runge-Kutta 4 method Solution: Given y′=x+5y, y(0)=-1, h=0.2, y(0.8)=? k1 = hf(x0,y0)=(0.2)f(0,-1)=(0.2)⋅(-5)=-1 k2 = hf(x0+h/2,y0+k1/2) = (0.2)f(0.1,-1.5) = (0.2)⋅(-7.4) = -1.48 k3 = hf(x0+h/2,y0+k2/2) = (0.2)f(0.1,-1.74) = (0.2)⋅(-8.6) = -1.72 k4 = hf(x0+h,y0+k3) = (0.2)f(0.2,-2.72) = (0.2)⋅(-13.4) = -2.68 y1 = y0 + 16(k1+2k2+2k3+k4) y1 = -1 + 16[-1+2(-1.48)+2(-1.72)+(-2.68)] y1=-2.68 ∴y(0.2)=-2.68 Again taking (x1,y1)in place of (x0,y0) and repeat the process k1 = hf(x1,y1) = (0.2)f(0.2,-2.68) = (0.2)⋅(-13.2) = -2.64 k2 = hf(x1+h2,y1+k1/2) = (0.2)f(0.3,-4) = (0.2)⋅(-19.7) = -3.94 k3 = hf(x1+h2,y1+k2/2) = (0.2)f(0.3,-4.65) = (0.2)⋅(-22.95) = -4.59 147 k4 = hf(x1+h,y1+k3) = (0.2)f(0.4,-7.27) = (0.2)⋅(-35.95) = -7.19 y2 = y1 + 16(k1+2k2+2k3+k4) y2 = -2.68 + 16[-2.64+2(-3.94)+2(-4.59)+(-7.19)] y2=-7.1617 ∴y(0.4)=-7.1617 Again taking (x2,y2) in place of (x0,y0) and repeat the process k1 = hf(x2,y2) = (0.2)f(0.4,-7.1617) = (0.2)⋅(-35.4083) = -7.0817 k2 = hf(x2+h2,y2+k1/2) = (0.2)f(0.5,-10.7025) = (0.2)⋅(-53.0125) = -10.6025 k3 = hf(x2+h2,y2+k2/2) = (0.2)f(0.5,-12.4629) = (0.2)⋅(-61.8146) = -12.3629 k4 = hf(x2+h,y2+k3) = (0.2)f(0.6,-19.5246) = (0.2)⋅(-97.0229) = -19.4046 y3 = y2 + 16(k1+2k2+2k3+k4) y3 = -7.1617 + 16[-7.0817+2(-10.6025)+2(-12.3629)+(-19.4046)] y3=-19.2312 ∴y(0.6)=-19.2312 Again taking (x3,y3) in place of (x0,y0) and repeat the process k1 = hf(x3,y3) = (0.2)f(0.6,-19.2312) = (0.2)⋅(-95.5559) = -19.1112 k2 = hf(x3+h2,y3+k1/2) = (0.2)f(0.7,-28.7868) = (0.2)⋅(-143.2339) = -28.6468 k3 = hf(x3+h2,y3+k2/2) = (0.2)f(0.7,-33.5546) = (0.2)⋅(-167.0728) = -33.4146 k4 = hf(x3+h,y3+k3) = (0.2)f(0.8,-52.6457) = (0.2)⋅(-262.4287) = -52.4857 y4 = y3 + 16(k1+2k2+2k3+k4) y4 = -19.2312 + 16[-19.1112+2(-28.6468)+2(-33.4146)+(-52.4857)] y4 = -51.8511 ∴y(0.8)=-51.8511 Check Your Progress-3 1. Find y(0.2) for y′=x2-y, x0=0,y0=1, with step length 0.1 using Runge-Kutta 4 method. 2. Find y(0.2) for y′=x2y -1, x0=0,y0=1, with step length 0.1 using Runge-Kutta 4 method 3.4 LET US SUM UP In this chapter we learnt different methods like Euler’s Method Runge-Kutta Second Order Method Runge-Kutta Fourth Order Method to obtain values of the function at any given point. 3.5 SUGGESTED ANSWERS FOR CHECK YOUR PROGRESS Check Your Progress-1 1. -0.7715 2. 0.9075 Check Your Progress-2 148  1. 0.7143 2. 0.819 Check Your Progress-3 1. 0.9052 2. 0.9003 3.6 GLOSSARY 1. Extrapolates means to extend by inferring unknown values from trends in the known data. 2. Erected means to put together and set upright. 3.7 ASSIGNMENT 1. Find y(0.2) for y’ = -(x2y2 + y), y(0) = 1.2, with step length 0.6. 2. Find y(0.4) for y’ = -xy2 , y(0) = 1, with step length 0.5. 3.8 ACTIVITY 1. Find y(0.3) for y’ = -y2 , y(0) = 1.9, with step length 0.4 3.9 FURTHER READING M.K. Jain, S.R.K. Iyengar and R.K. Jain, Numerical Methods for Scientific and Engineering Computation, 6th Ed., New age international publisher, India, 2007. S.S. Sastry, Introductory Method of Numerical Analysis, 5th Edition, PHI Learning Private Limited, New Delhi, 2012 149 UNIT 4: NUMERICAL SOLUTION OF HIGHER ORDER DIFFERENTIAL EQUATIONS Unit Structure 4.0 Higher Order Differential Equation 4.1 Let Us Sum Up 4.2 Suggested Answer for Check Your Progress 4.3 Glossary 4.4 Assignment 4.5 Activities 4.6 Further Readings 150 4.0 HIGHER ORDER DIFFERENTIAL EQUATION: A system of first order differential equations can be derived from a system of higher order differential equations. Consider the following second order differential equation as an example. d2y dy 2 = f ( x, y , ) dx dx With the initial conditions: dy y( x1 ) = y1 and ( ) x1 = y11. dx Now substitute, dy z= dx And; dz d 2 y = dx dx 2 d2y dy Hence the equation 2 = f ( x, y, ) can also be written as the system of first order dx dx differential equations; dy =z dx dy = g ( x, y , z ) dx With the initial conditions; y( x1 ) = y1 z ( x1 ) = y11 Example 4.1 Find y(0.1) for y′′=1+2xy-x2z, x0=0,y0=1,z0=0, with step length 0.1 using Runge-Kutta 2 method (2nd order derivative) Solution: Given y′′=1+2xy-x2z, y(0)=1, y′(0)=0, h=0.1, y(0.1)=? dy d 2 y dz put = z and differentiate w.r.t. x, we obtain = dx dx 2 dx We have system of equations 151 dy = z = f ( x, y , z ) dx dz = 1 + 2 xy − x 2 z = g ( x, y, z ) dx Method-1 : Using formula k2 = hf ( x0 + h, y0 + k1 ) Second order R-K method for second order differential equation k1 = h ⋅ f(x0,y0,z0) = (0.1) ⋅ f(0,1,0) = (0.1) ⋅ (0) = 0 l1 = h ⋅ g(x0,y0,z0) = (0.1) ⋅ g(0,1,0) = (0.1) ⋅ (1) = 0.1 k2 = h ⋅ f(x0+h,y0+k1,z0+l1) = (0.1) ⋅ f(0.1,1,0.1) = (0.1) ⋅ (0.1) = 0.01 l2 = h ⋅ g(x0+h,y0+k1,z0+l1) = (0.1) ⋅ g(0.1,1,0.1) = (0.1) ⋅ (1.199) = 0.1199 k1 + k2 y1 = y0 + = 1 + 0.005 = 1.005 2 ∴y(0.1)=1.005 Check Your Progress – 1 d2y dy 1. Given 2 +2 + y =0 dx dx dy with initial conditions y(0)=0, and |x = 0 = 1 dx Find the solution of the above system of equations in the interval [0,0.4] with interval size h=0.1 using the second order Runge-Kutta method. 2. Find y(0.2) for y” = xz2 – y2, x0 = 1, y0 = 1, z0 = 1, with step length 0.2 using Runge- Kutta 2 method (2nd order derivative). 4.1 LET’S SUM UP: In this chapter we learnt how do we find solution for second order derivatives. Generally by using second-order derivative. 4.2 Suggested Answer for Check Your Progress: Check Your Progress-1 1. 152  i 1 2 3 4 5 Xi 0.0 0.1 0.2 0.3 0.4 Yi 0.0 0.09 0.163 0.221 0.267 2. y(0.2) = 0.98 4.3 GLOSSARY: 1. Numerical Integration over more than one dimension is descried as cubature. 2. A function which is to be integrated is known as integrand. 4.4 ASSIGNMENT: 1. Find the solution of the following differential equation using the second order Runge-Kutta method d2y dy 2 − 2 − y = 2e x dx dx With y(3) = 1 for x = 3.1, 3.2, 3.3, 3.4 2. Find the solution of the following differential equation using the second order Ruge- kutta method d 2 y dy 3 2 − − 2 y = x2 dx dx 4.5 ACTIVITY: 1. Find y(0.2) for y” = x2z2 – y2,x0 = 0, y0 = 0, z0 = 0, with step length 0.2 using Runge-Kutta 2 method (2nd order derivative) 4.6 FURTHER READING: M.K. Jain, S.R.K. Iyengar and R.K. Jain, Numerical Methods for Scientific and Engineering Computation, 6th Ed., New age International Publisher, India, 2007. S.S. Sastry, Introductory Method of Numerical Analysis, 5th Edition, PHI Learning Private Limited, New Delhi, 2012. ………………….. 153