VIEWS: 53 PAGES: 48 CATEGORY: Engineering POSTED ON: 10/16/2012 Public Domain
1. Course No. & Title TEE-200 Computer Methods in Electrical Engineering 2. Credit break-up 2(2-2-0) 3. Pre-requisite TEE-150 Circuit Theory Catalog Description: Introduction to Computing, short round-up of digital computer architecture, flow chart, FORTRAN language, simple numerical methods and their applications to electrical problems, analog computers and their uses in electrical engineering problems. Syllabus: 1. Computer block diagram; 2. Importance of programming language; 3. High level languages; 4. Flow chart-the first step in problem solving by computers. 5. An illustrative FORTRAN-77 program to acquaint the students with various formats and types of instructions. 6. Process of running a programme - entering a program, compilation and execution. 7. Types of data; constants and variables; Real, Integer, Complex, Character, Logical and double precision. 8. Arithmetic operations on Data; 9. Priority of arithmetic operations. 10. IF statements - logical and arithmetic; 11. GOTO statements; illustrative examples; 12. Nested IF blocks. 13. DO structures; 14. Nested Do structures; 15. DO while and DO until structures; Illustrative examples. 16. Input/Output statements; 17. Formats; 18. Use of input/output data files. 19. Subscripted variables; 20. Use of subscripted variables for matrix operations. 21. Complex operations; Illustrative problems form Electrical Engineering. 22. Solution of algebraic equations using Gauss and Gauss-Seidel techniques; algorithm and programs. 23. Algorithm for solution of differential equations using step by step technique and 24. Runge-Kutta method. 25. Algorithm for solution of integration using Trapezoidal rule. Illustrative problems should be based on Electrical Engineering applications as far as practicable. TEXT BOOKS: 1. “Computer Programming in FORTRAN-77” by V. Raja Raman (Publisher: Prentice-hall of India Pvt Ltd) 2. “Introductory Methods of Numerical Analysis” by S S Sastry (Publisher: Prentice-hall Of India Pvt Ltd) 3. “Numerical Methods” by E Balagurusamy (Publisher: Tata Mgraw Hill) Evaluation: First Pre-final Exam 20 Marks Second Pre-final Exam 20 Marks Assignment/Quiz/Attendance/Lab 20 Marks Final Exam 40 Marks Total 100 Marks TEE-200: Computer Methods in Electrical Engineering (Lecture Notes): Abhishek Yadav 1. Computer block diagram At the heart of the computer is the microprocessor. This contains several REGISTERS to store data and an ARITHMETIC LOGIC UNIT (ALU) which manipulates data. It acts as the central processing unit (CPU) of the computer, carrying out a sequence of instructions, called a program. The program may be stored in memory, as software, or written into the memory from tape or disk. There are two types of memory. Read Only Memory (ROM) which stores software permanently. The software is not lost when the computer is switched off but the stored data cannot be changed. Random Access Memory (RAM) which can be written to and read from. The stored data is volatile. It is lost when the computer is switched off. The actual computer, its case and printed circuit boards etc are known as hardware. The computer needs to communicate with the outside world. It does this via interfaces which are usually a plug or socket of some type. 2. Importance of programming Language A programming language is a machine-readable artificial language designed to express computations that can be performed by a machine, particularly a computer. Programming languages can be used to create programs that specify the behavior of a machine, to express algorithms precisely, or as a mode of human communication. 3. High level languages A high-level programming language is a programming language with strong abstraction from the details of the computer. In comparison to low-level programming languages, it may use natural language elements, be easier to use, or be more portable across platforms. Such languages hide the details of CPU operations such as memory access models. This greater abstraction and hiding of details is generally intended to make the language user-friendly, as it includes concepts from the problem domain instead of those of the machine used. A high-level language isolates the execution semantics of computer architecture from the specification of the program, making the process of developing a program simpler and more understandable with respect to a low-level language. FORTRAN, PASCAL, C, C++, C#, JAVA are examples of high level languages. TEE-200: Computer Methods in Electrical Engineering (Lecture Notes): Abhishek Yadav 4. Flow chart-the first step in problem solving by computers A graphic representation of an algorithm used in the design phase of programming to work out the logical flow of a program. Connector or joining of two parts of program Annotation The following are some guidelines in flowcharting: a. In drawing a proper flowchart, all necessary requirements should be listed out in logical order. b. The flowchart should be clear, neat and easy to follow. There should not be any room for ambiguity in understanding the flowchart. c. The usual direction of the flow of a procedure or system is from left to right or top to bottom. d. Only one flow line should come out from a process symbol. e. Only one flow line should enter a decision symbol, but two or three flow lines, one for each possible answer, should leave the decision symbol. f. Only one flow line is used in conjunction with terminal symbol. g. Write within standard symbols briefly. As necessary, you can use the annotation symbol to describe data or computational steps more clearly. h. If the flowchart becomes complex, it is better to use connector symbols to reduce the number of flow lines. Avoid the intersection of flow lines if you want to make it more effective and better way of communication. i. Ensure that the flowchart has a logical start and finish. j. It is useful to test the validity of the flowchart by passing through it with a simple test data. TEE-200: Computer Methods in Electrical Engineering (Lecture Notes): Abhishek Yadav ADVANTAGES OF USING FLOWCHARTS The benefits of flowcharts are as follows: 1. Communication: Flowcharts are better way of communicating the logic of a system to all concerned. 2. Effective analysis: With the help of flowchart, problem can be analysed in more effective way. 3. Proper documentation: Program flowcharts serve as a good program documentation, which is needed for various purposes. 4. Efficient Coding: The flowcharts act as a guide or blueprint during the systems analysis and program development phase. 5. Proper Debugging: The flowchart helps in debugging process. 6. Efficient Program Maintenance: The maintenance of operating program becomes easy with the help of flowchart. It helps the programmer to put efforts more efficiently on that part Example 1 Draw a flowchart to find the sum of first 50 natural numbers. Example 2 Draw a flowchart to find the largest of three numbers A,B, and C. TEE-200: Computer Methods in Electrical Engineering (Lecture Notes): Abhishek Yadav Example 3 Draw a flowchart for computing factorial N (N!) . Where N! = 1 ´ 2 ´ 3 ´ …… N . 5. An illustrative FORTRAN-77 program: various formats and types of instructions Column 1: may contain a * or a C to indicate a comment. The rest of the line is ignored. A completely blank line is also treated as a comment. Comments are a good way of increasing the legibility of code, and should be used liberally. Columns 1-5: statement label or number. Occasionally you may need a label in order to refer to a particular line of text. Column 6: usually left blank. If a character is present, it implies that the line is being used to continue the statement on the previous line. The character itself is ignored. Columns 7-72: used for the Fortran statement itself. Columns 73-80: ignored. Problem: Write a program in FORTRAN-77 to compute the area of a circle. program circlearea real r, area, pi parameter (pi = 3.14159) c This program computes the area of a circle. print *, "What is the radius?" read *, r area = pi * r ** 2 print *, "The area is", area print *, "Bye!" end 6. Process of running a programme - entering a program, compilation and execution Click on “Start” Button; go to “Run” and type “cmd” to open a command window. Type “edit filename.for” and press enter key. Type your program, save it end exit. To compile the program, type “FL filename.for” and press enter key. To execute the program, type “filename” and press the enter key. 7. Types of data; constants and variables; Real, Integer, Complex, Character, Logical and double precision TEE-200: Computer Methods in Electrical Engineering (Lecture Notes): Abhishek Yadav Constants are quantities whose values do not change through program execution. Variables are used in FORTRAN to indicate the values that may change during the program execution. FORTRAN supports six different data types: 1. Integer 4. Complex 2. Real 5. Character 3. Double precision 6. Logical The first four types are used for processing various types of NUMERIC data. Character type is used for processing STRINGS of characters. Logical type is for processing LOGICAL data; such data may have either .TRUE. or .FALSE. as its value. Each of these data types has two representations: they can be CONSTANTS or VARIABLES. Constants: INTEGER constants are a sequence of digits without commas or decimal points. A negative integer constant must be preceded by a minus sign, a plus sign is optional for positive integers. VALID: 0 137 -2517 +3245 INVALID: 3,456 24.5 REAL constants are numbers that have a decimal point without a comma. The rules for negative and positive signs also apply here. VALID: 0. 12345.678 12.0 -.0001 INVALID: 0 1,234.99 22 REAL constants can also be represented in scientific notation where the exponent is raised to a power of ten and multiplied by the base. The following real number can be expressed in many ways: Eg; 3.374 X 102 = 3.374E2 .3374E3 337.4E0 33740E-2 Notice in the last example the decimal was omitted. Scientific Notation automatically implies a REAL CONSTANT. CHARACTER constants, or STRINGS, are sequences of characters or symbols enclosed within apostrophes (single quotes). VALID: 'This is a STRING' ' ' '&h"y^5)' INVALID: '' "OC3030" Variables: Variable names must be between 1 and 32 letters or numbers in length with the first character being a letter. VALID: A AB Z1 MASS LENGTH ABC123 NOTTOOLONG INVALID: 1VAR R2-D2 A*AB4 6FEET Variables must be one of the six types discussed earlier. The type of the variable determines what kind of data can be assigned to it. Variables must be declared at the top of the program before any executable statements. REAL ALPHA, BETA, KOUNT INTEGER Z, A1 CHARACTER*8 HGT, STRING*13, L9*2 TEE-200: Computer Methods in Electrical Engineering (Lecture Notes): Abhishek Yadav If a variable is not explicitly declared, it will be assigned a type based on the default naming convention. This convention states that the variable is: INTEGER if the first letter is I,J,K,L,M,N (I-N) REAL any other letter NOTE: There is no default convention for CHARACTER, LOGICAL, DOUBLE PRECISION or COMPLEX types. 8. Arithmetic operations on Data FORTRAN variables and constants can be processed using operations and functions appropriate to their types. The five arithmetic operators in FORTRAN are: 1. Addition + 2. Subtraction - 3. Multiplication * 4. Division / 5. Exponentiation ** Example: B 2 - 4 A C would be written B**2 - 4.0*A*C 9. Priority of arithmetic operations Expressions are always carried out based on a set of priority rules: 1. All exponentiations are performed first; consecutive exponentiation is performed right to left. 2. All multiplication and division is performed next, in order from left to right. 3. Addition and subtraction is performed last, in order from left to right. Example 1: 2 ** 3 ** 2 = 2 ** 9 = 512 Example 2: 5 * 11 - 5 ** 2 * 4 + 9 = 5 * 11 - 25 * 4 + 9 = 55 - 100 + 9 = -36 Parentheses can modify the order of evaluation. Expressions within parentheses will be evaluated first, using the priority rules. If parentheses are nested, the inner expression is evaluated first. Using parentheses is recommended (even if their use does not change the order of evaluation) because it makes the expression easier to read. 10. IF statements - logical and arithmetic The logical IF statement conditionally executes a single Fortran statement based on the current value of a logical expression within the logical IF statement. It takes the following form: IF (e) st e is a logical expression. st is any complete, executable Fortran statement-except any of the block IF statements, DO, END DO, another logical IF statement, a SELECT CASE, CASE, or CASE DEFAULT statement. The logical IF statement first evaluates the logical expression e and then acts as follows: If e is true, st executes. If e is false, control transfers to the next executable statement after the logical IF; st does not execute. The following examples show valid logical IF statements: IF (J.GT.4 .OR. J.LT.1) GO TO 250 IF (REF(J,K) .NE. HOLD) REF(J,K) = REF(J,K) * (-1.5D0) IF (ENDRUN) CALL EXIT TEE-200: Computer Methods in Electrical Engineering (Lecture Notes): Abhishek Yadav The arithmetic IF statement conditionally transfers control to one of three statements, based on the value of an arithmetic expression. The arithmetic IF statement takes the following form: IF (expr) label1, label2, label3 expr: a scalar numeric expression of type integer or real (enclosed in parentheses). label1, label2, label3: the labels of valid branch target statements that are in the same scoping unit as the arithmetic IF statement. All three labels are required, but they do not need to refer to three different statements. The same label can appear more than once in the same arithmetic IF statement. During execution, the expression is evaluated first. Depending on the value of the expression, control is then transferred as follows: If the Value of expr is: Control Transfers To: Less than 0 Statement label1 Equal to 0 Statement label2 Greater than 0 Statement label3 The following example transfers control to statement 50 if the real variable THETA is less than or equal to the real variable CHI. Control passes to statement 100 only if THETA is greater than CHI. IF (THETA-CHI) 50,50,100 The following example transfers control to statement 40 if the value of the integer variable NUMBER is even. It transfers control to statement 20 if the value is odd. IF (NUMBER / 2*2 - NUMBER) 20,40,20 11.GOTO statements Unconditional GOTO Statement GOTO label When an unconditional GOTO is encountered, control is immediately transferred to the executable statement labelled with the label. Computed GOTO Statement GOTO (label1, label2, ..., labeln), integer-expression This form of GOTO statement is obscure and its use in modern programs is strongly discouraged. It is equivalent to this block IF statement: IF (integer-expression .EQ. 1) THEN GOTO label1 ELSE IF (integer-expression .EQ. 2) THEN GOTO label2 ... ELSE IF (integer-expression .EQ. n) THEN TEE-200: Computer Methods in Electrical Engineering (Lecture Notes): Abhishek Yadav GOTO labeln ENDIF If the integer-expression is less than 1 or greater than n, then control passes on and no GOTO is executed. It is permissible for two or more of the labels to be the same. 12.Nested IF blocks if statements can be nested in several levels. To ensure readability, it is important to use proper indentation. Here is an example: if (x .GT. 0) then if (x .GE. y) then write(*,*) 'x is positive and x >= y' else write(*,*) 'x is positive but x < y' endif elseif (x .LT. 0) then write(*,*) 'x is negative' else write(*,*) 'x is zero' endif 13.DO structures The do-loop is used for simple counting. Here is a simple example that prints the cumulative sums of the integers from 1 through n: integer i, n, sum n = 10 sum = 0 do i = 1, n sum = sum + i write(*,*) 'i =', i write(*,*) 'sum =', sum enddo 14.Nested Do structures One DO loop may be completely contained within another DO loop. The two DO loops must use different control variables. : DO statement1 DO statement2 END DO statement3 END DO 15.DO while and DO until structures The most intuitive way to write a while-loop is while (logical expr) do statements enddo TEE-200: Computer Methods in Electrical Engineering (Lecture Notes): Abhishek Yadav or alternatively, do while (logical expr) statements enddo The statements in the body will be repeated as long as the condition in the while statement is true. If the termination criterium is at the end instead of the beginning, it is often called an until-loop. The pseudocode looks like this: do statements until (logical expr) 16.Input/Output statements Input statements provide the means of transferring data from external media to internal storage or from an internal file to internal storage. This process is called reading. Output statements provide the means of transferring data from internal storage to external media or from internal storage to an internal file. This process is called writing. Some input/output statements specify that editing of the data is to be performed. In addition to the statements that transfer data, there are auxiliary input/output statements to manipulate the external medium, or to inquire about or describe the properties of the connection to the external medium. There are nine input/output statements: 1. READ 2. WRITE 3. PRINT 4. OPEN 5. CLOSE 6. INQUIRE 7. BACKSPACE 8. ENDFILE 9. REWIND The READ, WRITE, and PRINT statements are data transfer input/output statements). The OPEN, CLOSE, INQUIRE, BACKSPACE, ENDFILE, and REWIND statements are auxiliary input/output statements. The BACKSPACE, ENDFILE, and REWIND statements are file positioning input/output statements. The simplest form of the I/O statement is the list-directed form which is represented by: READ(*,*) item1, item2, item3... WRITE(*,*) item1, item2, item3... where ITEMx = a variable, a constant or math expression Example: WRITE(*,*) 'ALPHA=', ALPHA The first asterisk (*) means the input comes from the keyboard in a READ statement and goes to the screen in a WRITE statement. The second asterisk (*) means the computer decides how the I/O elements should look based on the TYPE of data in the input/output list. This is sometimes called "FREE-FORMAT". NOTES ON LIST-DIRECTED I/O TEE-200: Computer Methods in Electrical Engineering (Lecture Notes): Abhishek Yadav SPACES may be inserted between components of the statement to add clarity. The READ statement causes the program to PAUSE and allow you to enter values. The program will not continue until all values have been entered. Separate values with SPACES when typing data into the program with a READ operation. Make sure you press the ENTER key. Insert a WRITE statement to PROMPT yourself for input just before a READ statement. This statement should tell the user WHAT to enter. GENERAL FORM OF I/O STATEMENTS The general form of the FORTRAN I/O statements allow data transfer to FILES, TAPE, PRINTER and other devices as well as the TERMINAL. The general form is: WRITE(unit#, format, options) item1, item 2,... READ(unit#, format, options) item1, item2,... NOTE: We will restrict our use to TERMINAL and FILE I/O. Typically, you transfer data files to/from tape, diskette, and printer by using UNIX (or DOS on PCs) commands, rather than a Fortran program reading or writing directly from/to the device. In this form, PARENTHESES are used to enclose information about the UNIT, the FORMAT (if any) and other options. Again, SPACES may be used to add clarity to the statement. The input and output lists (item1, item2,...) are composed of constants, variables or expressions, separated by COMMAS. 17.Formats The FORMAT IDENTIFIER as used in a WRITE or READ statement generally has one of two forms; 1. An asterisk (*), indicates LIST-DIRECTED or "free format". 2. A LABEL designates a FORMAT statement which specifies the format to use. For simple I/O, the * is used as we have already discussed. When you want to FORMAT your I/O, you will generally use method 2. - The FORMAT statement is a list of FORMAT DESCRIPTORS separated by commas which describe what each FIELD of the output should look like. Example: WRITE(*,10) 'USING', L2, 'AREA =', AREA 10 FORMAT(1X, A5, 2X, I3, 4X, A6, 2X, F6.2) In a format statement, there is a one-to-one correspondence with the items in the I/O list, with the exception of certain positional descriptors; X, T, /. The FORMAT statement is defined only once (for each label referenced) in the program but may be used by any number of I/O statements. FORMAT DESCRIPTORS - The elements in the I/O list must agree in TYPE with the FORMAT DESCRIPTOR being used. - There are a dozen or so descriptors, The most commonly used are: Descriptor Use rIw Integer data rFw.d Real data in decimal notation TEE-200: Computer Methods in Electrical Engineering (Lecture Notes): Abhishek Yadav rEw.d Real data in scientific notation Aw Character data 'x...x' Character strings nX Horizontal spacing (skip spaces) / Vertical spacing (skip lines) Tc Tab Where: w = positive integer specifying FIELD WIDTH r = positive integer REPEAT COUNT d = non-negative integer specifying number of digits to display to right of decimal. x = any character n = positive integer specifying number of columns c = positive integer representing column number NOTES ON DESCRIPTORS: - Values for I, F, and E descriptors will be right justified. - You must leave enough space for negative sign in the field width for I, E, F, and for decimal point for E, F and for exponent for E. - When using the E descriptor, the number will be NORMALIZED. This means the number is shifted over so that the first digit is in the tenths position. A leading zero and decimal are always present and a minus sign, if needed. The exponent always requires 4 spaces; E, the sign and 2 digits. - If you overflow the space allotted to a descriptor, it will print out asterisks (at RUN TIME) in the field space. - Character data (A) is also right-justified but if it overflows the allotted space it is left-most truncated. - The first character of each line of output is used to control the vertical spacing. Use 1X to give normal spacing. 18.Use of input/output data files The UNIT is a number which has an association with a particular device. The device can be the TERMINAL or a FILE (or something else too). The UNIT is an INTEGER or INTEGER EXPRESSION, generally between 1-30. Standard FORTRAN reserves two UNIT numbers for I/O to user. They are: UNIT = 5 for INPUT from the keyboard with the READ statement UNIT = 6 for OUTPUT to the screen with the WRITE statement Most versions of FORTRAN will also let you use the ASTERISK (*) for I/O to the TERMINAL. The asterisk can be used with both the READ and WRITE statements, thus there is no need to remember whether 5 or 6 is for input or output. When I/O is to a file you must ASSOCIATE a UNIT number (which you choose) with the FILENAME. Use any unit number other than 5 and 6. On some computers, some unit numbers are reserved for use by the computer operating system. The association of the unit number and filename occurs when the OPEN statement is executed. OPEN(UNIT=n, FILE='filename', options...) TEE-200: Computer Methods in Electrical Engineering (Lecture Notes): Abhishek Yadav This statement associates UNIT n with the file mentioned. All subsequent READs or WRITEs using unit n will be to or from this file. MISCELLANEOUS FILE I/O NOTES When doing I/O to a FILE, each READ statement inputs data from a NEW LINE and each WRITE statement outputs data on a NEW LINE. Most files are SEQUENTIALLY organized. You can't easily go back up in the file, but you can REWIND the file with: REWIND unit_number When you are finished with the file, you may CLOSE it with: close(10) or close(unit=10) (This is optional.) 19.Subscripted variables It is often useful to be able to store several values in one variable, for example, to store the three components of a vector, or one hundred data points from an experiment. Fortunately, in Fortran we are able to do this, by using arrays. Arrays are sometimes referred to as subscripted variables. They may be declared as follows: integer n (10) double precision e (25) In the above, the integer variable n refers to 10 separate variables: n(1), n(2), n(3), ... , n(10). Similarly, the double precision variable e refers to 25 separate variables: e(1), e(2), ... , e(25). The number in brackets is referred to as the index or subscript of the array. The useful feature about such subscripts is that they may also be variables or expressions. For example, consider the following fragment of code: integer i, j, k double precision a(10) c i = 1 j = 2 k = 3 c a(7) = 3.0d0 a(i) = 1.5d1 a(i*j+k) = 2.0d-3 In the above, three elements of the array a are initialised: a(7), a(1) and a(5). The range of subscript values for a particular array may be chosen when the array is declared, by using a more flexible form of declaration. For example: integer a(0:100) double precision b(-10:10) Here, a takes subscripts from 0 to 100 (101 elements) whereas b takes subscripts from -10 to +10 (21 elements). As well as simple one-dimensional arrays, multi-dimensional arrays can be used in Fortran. The most common example is a two-dimensional array, eg. a(n,m), which is useful for storing matrices. In this case, the element a(1,1) is the value for the first row, first column; a(1,2) is the value for the first row, second column; etc... 20.Use of subscripted variables for matrix operations TEE-200: Computer Methods in Electrical Engineering (Lecture Notes): Abhishek Yadav Matrices are useful for the problems related with linear algebra, networks, population dynamics, Markov processes and many similar issues. In this chapter we deal with multidimensional arrays, specifically focusing on commonly used two-dimensional arrays. A matrix is a two-dimensional array which may be used in a wide variety of representations. For example, a distance array representing the lengths of direct connections in a network is a matrix. We will deal mainly with square matrices in this chapter (i.e. matrices having the same number of rows as columns), although in principle a matrix can have any number of rows or columns. A matrix with only one column is also called a vector. A matrix is usually denoted by a bold capital letter, e.g. A. Each entry, or element, of the matrix is denoted by the small letter of the same name followed by two subscripts, the first indicating the row of the element, and the second indicating the column. So a general element of the matrix A is called aij, meaning its placement at row i and column j. If A has three rows and columns —(3 times 3) for short—it will look like this in general: a11 a 12 a 12 A= a 21 a 22 a 23 a 31 a 32 a 33 Various mathematical operations are defined on matrices. Addition and subtraction are obvious, and may be done with the intrinsic operators in Fortran 90. So the matrix addition [A = B + C] translates directly into A=B+C where the arrays must clearly all have the same shape, i.e. the same extent along corresponding dimensions. Matrix multiplication is used widely in network theory, solution of linear systems of equations, transformation of co-ordinate systems, and population modeling, for examples. Product C = A B of a p by q matrix A by a q by r matrix B is defined by q cij = k 1 aik bkj , for all cij , i=1 ,…, p , j=1, … , r As seen in this example, A.B is not always equal to B.A. A Fortran program would calculate the matrix multiplications in do loops. Following example makes a scalar matrix multiplication, a matrix by matrix multiplication, and inverse of a matrix, which may be useful in solving linear sets of equations. REAL A(4,8), B(4,4), C(4,4), Ainv(4,4) REAL TM(8), TC(4) ! temporary vectors REAL PE, TE INTEGER I,J,K, PR, TR A(1,1:2) = (/ 1, 0 /) A(2,1:2) = (/ 2, 1 /) DATA (( B(I,J), J = 1,2), I = 1,2) / 1, 2, 0, 1 / ! Multiply matrices C = A * B C=0 do I = 1,2 do J = 1,2 do K=1,2 C(I,J) = C(I,J) + A(I,K) * B(K,J) end do end do end do PRINT*, " A * B = "; do i=1,2; PRINT*, (C(I,J), J=1,2 ) ; end do ! Multiply matrices C= B*A C=0 do I = 1,2; do J = 1,2 do K=1,2 TEE-200: Computer Methods in Electrical Engineering (Lecture Notes): Abhishek Yadav C(I,J) = C(I,J) + B(I,K) * A(K,J) end do end do; end do PRINT*, " B * A = "; do i=1,2; PRINT*, (C(I,J), J=1,2 ) ; end do ! We will solve a linear set of equations !Ax=b ! to find the values of x. n=4 A(1,1:N)= (/ 1,2,5,6 /) ! assignment goes columnwise A(2,1:N)= (/ 3,2,2,6 /) A(3,1:N)= (/ 2,6,5,4 /) A(4,1:N)= (/ 4,3,2,5 /) ! a is now a 4 by 4 matrix print*,"A="; do i=1,4; PRINT*, A(i,1:4) ; end do B(1:N,1)= (/ 0, 1, 0, 2 /) ! b is a column vector ! now solve equation sum(a(i,j)*x(j), j=1,n)=b(i), i=1,n ! select a pivot, use pivot to eliminate the elements ! of the selected colon. ! Inverse of a matrix is obtained if a(1:n,n:2n)=I(nxn) do i=1,N; do j=1,N; a(i,j+N)=0; end do; a(i,i+N)=1; end do; DO PR = 1, n ! process every row PE = A( PR, PR ) ! choose pivot element IF (PE == 0) THEN ! check for zero pivot k = PR + 1 ! run down rows to find a non-zero pivot DO WHILE (PE == 0 .AND. k <= n) PE = A( k, PR ) ! try next row k = k + 1 ! K will be 1 too big end do if (PE == 0) then ! it's still zero print*, "Couldn't find a non-zero pivot" else ! non-zero pivot in row K, so swop rows PivRow and K: TM = A( PR, 1:2*n ) K = K - 1 ! adjust for overcount A( PR, 1:2*n ) = A( K, 1:2*n ) A( K, 1:2*N ) = TM end if end if A( PR, 1:2*n ) = A( PR, 1:2*n ) / PE ! scalar divide ! whole row ! now replace all other rows by target row minus pivot row ... ! ... times element in target row and pivot column: do TR = 1, n TE = A( TR, PR ) if (TR /= PR) then A( TR, 1:2*N ) = A( TR, 1:2*N ) & - A( PR, 1:2*N ) * TE end if end do end do TEE-200: Computer Methods in Electrical Engineering (Lecture Notes): Abhishek Yadav ! second half of a() is inv(a) AINV(1:n,1:n)=A(1:n,n+1:2*n) print*," AINV="; do i=1,n; print*, AINV(i,1:n); end do ! Multiply C = AINV*b to find x C=0 do I = 1,n; do J = 1,1 do K=1,n C(I,J) = C(I,J) + AINV(I,K) * B(K,J) end do end do; end do ! print calculated x to the screen print*," X="; print*, C(1:n,1) end 21.Complex operations A complex number has a real and an imaginary component. In Fortran, complex numbers are stored as a pair of real numbers (real part first, followed by the imaginary part), or a pair of double-precision numbers. The standard way to declare a complex number is simply using the COMPLEX statement: COMPLEX myVariable, anotherVariable, anArray(2,3) COMPLEX*8 myVariable, anotherVariable, anArray(2,3) If you want double-precision complex numbers, you're pretty much stuck with specifying a precision, and hoping that the compiler likes that format. Here's where portability comes in: If, for instance, the machine has floating point types of lengths 4, 8, and 16, and I specify a complex length of 8, 16, or 32, I'm likely to get a pair of floating point numbers at the size I expect. If, however, I specify a length of some other value, for instance, 10, then I'm likely to get a pair of 4s or 8s, depending on what the compiler likes to do. So, in general, to get the values as double-precision, I'd code: COMPLEX*16 myVariable, anotherVariable, anArray(2,3) Constants and Expressions Complex constants are specified as (, a real value, a comma, another real value, and ). Some examples are 1 as (1.0,0.0), i as (0.0,1.0), 2-3i as (2.0,-3.0), and 1000000000i as (0.0,1.0E+09). The same constants can be coded as double-precision complex constants by the simple expedient of using a D in the exponent. Thus, the same constants can be coded in double-precision using (1.0D+00,0.0D+00), (0.0D+00,1.0D+00), (2.0D+00,-3.0D+00), and (0.0D+00,1.0D+09), respectively. Any expression involving complex numbers and other numbers is promoted to complex. Operators and Functions All of the arithmetic operators can take a complex number on either side. Care must be taken in using function calls that might cause an error. For instance, taking the square root of the real number -1 (coded as -1.0) will result in an error, because -1 is outside of the domain of real square root. Taking the square root of the complex number -1 (coded as (-1.0,0.0)) is OK, because -1 is in the domain of complex square root. 22.Solution of algebraic equations using Gauss and Gauss-Seidel techniques In numerical linear algebra, the Gauss–Seidel method, also known as the Liebmann method or the method of successive displacement, is an iterative method used to solve a linear system of equations. It is named after the German mathematicians Carl Friedrich Gauss and Philipp Ludwig von Seidel. The Gauss TEE-200: Computer Methods in Electrical Engineering (Lecture Notes): Abhishek Yadav method or Jacobi method is easily derived by examining each of the n equations in the linear system Ax = b n in isolation. If in the i-th equation a x j 1 ij j bi we solve for the value of xi while assuming the other entries of x remain fixed, we obtain This suggests an iterative method defined by which is the Jacobi method or Gauss method. n Consider the linear equations in ( aij x j bi ). If we proceed as with the Gauss method, but now assume j 1 that the equations are examined one at a time in sequence, and that previously computed results are used as soon as they are available, we obtain the Gauss-Seidel method: 23.Algorithm for solution of differential equations using step by step technique What is Euler’s Method? Euler’s method is a numerical technique to solve ordinary differential equations of the form f x, y , y0 y0 dy (1) dx So only first order ordinary differential equations can be solved by using Euler’s method. In other sections, we will discuss how Euler’s method is used to solve higher order ordinary differential equations or coupled (simultaneous) differential equations. How does one write a first order differential equation in the above form? Example 1: 2 y 1.3e x , y0 5 dy dx is rewritten as 1.3e x 2 y, y0 5 dy dx In this case f x, y 1.3e x 2 y Another example is given as follows Example 2: x 2 y 2 2Sin(3x ), y 0 5 dy ey dx TEE-200: Computer Methods in Electrical Engineering (Lecture Notes): Abhishek Yadav is rewritten as dy 2 Sin(3x ) x 2 y 2 , y 0 5 dx ey In this case 2 Sin(3x ) x 2 y 2 f x, y ey At x 0 , we are given the value of y y 0 . Let us call x 0 as x 0 . Now since we know the slope of y with respect to x , that is, f x, y , then at x x0 , the slope is f x0 , y 0 . Both x 0 and y 0 are known as they from the initial condition y x0 y 0 . y True value y1, x0,y0 Φ Predicted value Step size, h x Figure 1. Graphical interpretation of the first step of Euler’s method So the slope at x=x0 as shown in Figure 1 is Rise Slope Run y y0 1 x1 x 0 f x0 , y 0 From here y1 y 0 f x0 , y 0 x1 x0 Calling x1 x0 as a step size h , we get y1 y 0 f x0 , y 0 h . (2) One can now use the value of y1 (an approximate value of y at x x1 ) to calculate y 2 , and that would be the predicted value at x 2 , y 2 y1 f x1 , y1 h x2 x1 h Based on the above equations, if we now know the value of y yi at x i , then yi 1 yi f xi , yi h (3) This formula is known as the Euler’s method and is illustrated graphically in Figure 2. In some books, it is also called the Euler-Cauchy method. TEE-200: Computer Methods in Electrical Engineering (Lecture Notes): Abhishek Yadav y True Value yi+1, Predicted value Φ yi h Step size x xi xi+1 Figure 2. General graphical interpretation of Euler’s method Example 3: A rectifier-based power supply requires a capacitor to temporarily store power when the rectified waveform from the AC source drops below the target voltage. To properly size this capacitor a first- order ordinary differential equation must be solved. For a particular power supply, with a capacitor of 150 μF, the ordinary differential equation to be solved is dv(t ) 1 18 cos(120 (t )) 2 v(t ) 6 0.1 max ,0 , v(0) 0 dt 150 10 0.04 Find voltage across the capacitor at t= 0.00004s. Use step size h=0.00002 Solution dv 1 18 cos(120 (t )) 2 v 0.1 max ,0 , dt 150 10 6 0.04 18 cos(120 (t )) 2 v f t , v 1 0.1 max ,0 150 10 6 0.04 Per equation 3, the Euler’s method reduces to vi 1 vi f t i , vi h For i 0 , t 0 0 , v0 0 v1 v0 f t 0 , v0 h 0 f 0,00.00002 1 18 cos(120 (0)) 2 (0) 0.1 max ,0 0.00002 150 10 6 0.04 0 2.6667 10 0.00002 6 53.322V v1 is the approximate voltage at t t1 t 0 h 0 0.00002 0.00002 v1 v0.00002 53 .322V For i 1 , t1 0.00002 , v1 53 .322 TEE-200: Computer Methods in Electrical Engineering (Lecture Notes): Abhishek Yadav v2 v1 f t1 , v1 h 53 .322 f 0.00002 ,53 .322 0.00002 1 18 cos(120 (0.00002 )) 2 (53 .322 ) 53 .322 0.1 max ,0 0.00002 150 10 6 0.04 53 .322 0.000015000 0.00002 53.322V v 2 is the approximate voltage at t t 2 t t 2 t1 h 0.00002 0.00002 0.00004s v2 v0.00004 53 .322V Figure 3 compares the exact solution with the numerical solution from Euler’s method for the step size of h=0.00004s. 60 h=0.00002 50 Voltage, v(V) 40 30 Exact Solution 20 10 0 0 0.00001 0.00002 0.00003 0.00004 Time, t(sec) Figure 3. Comparing exact and Euler’s method The problem was solved again using smaller step sizes. The results are given below in Table 1. Table 1. Voltage at 0.00004 seconds as a function of step size, h Step size, v0.00004 Et |t | % h 106.64 -90.667 567.59 0.00004 53.307 -37.333 233.71 0.00002 26.640 -10.666 66.771 0.00001 15.995 - 0.13146 0.000005 15.992 0.021000 0.11268 0.0000025 - 0.018000 Figure 4 shows how the voltage varies as a function of time for different step sizes TEE-200: Computer Methods in Electrical Engineering (Lecture Notes): Abhishek Yadav 100 h=0.00004 v(V) 80 Voltage, 60 h=0.00002 40 h=0.00001 20 Exact solution 0 0 0.00001 0.00002 0.00003 0.00004 Tim e, t (sec) Figure 4. Comparison of Euler’s method with exact solution for different step sizes while the values of the calculated voltage at t=0.00004s as a function of step size are plotted in Figure 5. Voltage, v(V) Step size, h (s) 0 0 0.00001 0.00002 0.00003 0.00004 Figure 5. Effect of step size in Euler’s method. The solution to this nonlinear equation is v(0.00004) 15.974V It can be seen that Euler’s method has large errors. This can be illustrated using Taylor series. 2 3 y i 1 y i dy xi 1 xi 1 d 2 xi 1 xi 2 1 d y xi 1 xi 3 ... (4) y dx xi , yi 2! dx x , y 3! dx 3 x , y i i i i yi 1 yi f ( xi , yi ) f ' ( xi , yi )xi 1 xi f ' ' ( xi , yi )xi 1 xi ... 1 2 1 3 (5) 2! 3! As you can see the first two terms of the Taylor series yi 1 yi f xi , yi h are the Euler’s method. The true error in the approximation is given by f xi , yi 2 f xi , yi 3 Et h h ... (6) 2! 3! TEE-200: Computer Methods in Electrical Engineering (Lecture Notes): Abhishek Yadav The true error hence is approximately proportional to the square of the step size, that is, as the step size is halved, the true error gets approximately quartered. However from Table 1, we see that as the step size gets halved, the true error only gets approximately halved. This is because the true error being proportioned to the square of the step size is the local truncation error, that is, error from one point to the next. The global truncation error is however proportional only to the step size as the error keeps propagating from one point to another. 24.Runge-Kutta method. What is the Runge-Kutta 2nd Order Method? Runge-Kutta 2nd order method is a numerical technique to solve ordinary differential equation of the form f x, y , y0 y0 dy dx Only first order ordinary differential equations can be solved by using Runge-Kutta 2nd order method. In other sections, we have discussed how Euler and Runge-Kutta methods are used to solve higher order ordinary differential equations or coupled (simultaneous) differential equations. How does one write a first order differential equation in the above form? Example 1: 2 y 1.3e x , y0 5 dy dx is rewritten as 1.3e x 2 y, y0 5 dy dx In this case f x, y 1.3e x 2 y Another example is given as follows Example 2: x 2 y 2 2Sin(3x ), y 0 5 dy ey dx is rewritten as dy 2Sin(3x) x 2 y 2 , y 0 5 dx ey In this case 2 Sin(3 x) x 2 y 2 f x, y ey Euler’s method is given by yi 1 yi f xi , yi h (1) where x0 0 y0 y0 h xi 1 xi To understand Runge-Kutta 2nd order method, we need to derive Euler’s method from Taylor series. TEE-200: Computer Methods in Electrical Engineering (Lecture Notes): Abhishek Yadav 2 3 y i 1 y i dy xi 1 xi 1 d y 2 xi 1 xi 2 1 d y 3 xi 1 xi 3 ... dx xi , yi 2! dx xi , yi 3! dx xi , yi f ' ( xi , yi )xi 1 xi f ' ' ( xi , yi )xi 1 xi ... (2) 1 1 yi 1 yi f ( xi , yi ) 2 3 2! 3! As you can see the first two terms of the Taylor series yi 1 yi f xi , yi h are the Euler’s method and hence can be considered to be Runge-Kutta 1st order method. The true error in the approximation is given by f xi , yi 2 f xi , yi 3 Et h h ... (3) 2! 3! So how would a 2nd order method formula look like. It would include one more term of the Taylor series as follows. yi 1 yi f xi , yi h f xi , yi h 2 1 (4) 2! Let us take a generic example of a first ordinary differential equation e 2 x 3 y, y 0 5 dy dx f x, y e2 x 3 y Now since y is a function of x, f x, y f x, y dy f x, y (5) x y dx e 2 x 3 y x 2 x y e 3 y e 2 x 3 y 2e ( 3)e 3 y 2 x 2 x 5e2 x 9 y The 2nd order formula for the above example would be yi 1 yi f xi , yi h f xi , yi h 2 1 2! yi 1 yi e 2 xi 3 yi h 5e 2 xi 9 yi h 2 1 2! However, we already see the difficulty of having to find f x, y in the above method. What Runge and Kutta did was write the 2nd order method as yi 1 yi a1k1 a 2 k 2 h (6) where k1 f xi , yi k2 f xi p1h, yi q11k1h (7) This form allows one to take advantage of the 2 order method without having to calculate f x, y . nd So how do we find the unknowns a1 , a2 , p1 and q11 . Without proof, equating Equation (4) and (6) , gives three equations. a1 a2 1 1 a2 p1 2 1 a2q11 2 Since we have 3 equations and 4 unknowns, we can assume the value of one of the unknowns. The other three will then be determined from the three equations. Generally the value of a2 is chosen to TEE-200: Computer Methods in Electrical Engineering (Lecture Notes): Abhishek Yadav 1 2 , 1 and , and are evaluate the other three constants. The three values generally used for a2 are 2 3 known as Heun’s Method, Midpoint method and Ralston’s method, respectively. Heun’s Method 1 Here a2 is chosen, giving 2 1 a1 2 p1 1 q11 1 resulting in 1 1 yi 1 yi k1 k2 h (8) 2 2 where k1 f x i , y i (9a) k 2 f x i h, y i k 1 h (9b) This method is graphically explained in Figure 1. y Slope f xi h, yi k1h yi+1, predicted Slope f xi , yi Average Slope 1 f xi h, yi k1h f xi , yi 2 yi x xi xi+1 Figure 1. Runge-Kutta 2nd order method (Heun’s method) Midpoint Method Here a2 1 is chosen, giving a1 0 1 p1 2 1 q11 2 resulting in yi 1 yi k2h (10) where TEE-200: Computer Methods in Electrical Engineering (Lecture Notes): Abhishek Yadav k1 f xi , yi (11a) 1 1 k 2 f xi h, y i k1h (11b) 2 2 Ralston’s method 2 Here a 2 is chosen, giving 3 1 a1 3 3 p1 4 3 q11 4 resulting in 1 2 yi 1 yi ( k1 k 2 )h (12) 3 3 where k1 f xi , yi (13a) 3 3 k 2 f x i h , y i k1 h (13b) 4 4 Example 3: A rectifier-based power supply requires a capacitor to temporarily store power when the rectified waveform from the AC source drops below the target voltage. To properly size this capacitor a first- order ordinary differential equation must be solved. For a particular power supply, with a capacitor of 150 μF, the ordinary differential equation to be solved is dv(t ) 1 18 cos(120 (t )) 2 v(t ) 0.1 max ,0 , v(0) 0 dt 150 10 6 0.04 Find voltage across the capacitor at t= 0.00004s. Use step size h=0.00002 Solution dv 1 18 cos(120 (t )) 2 v 0.1 max ,0 , dt 150 10 6 0.04 18 cos(120 (t )) 2 v f t , v 1 0.1 max ,0 150 10 6 0.04 Per Heun’s method given by Equations (8) and (9) 1 1 vi 1 vi k1 k 2 h 2 2 k1 f t i , vi k 2 f t i h, vi k1 h i 0, t 0 0, v0 v(0) 0 k1 f t 0 , vo f 0,0 TEE-200: Computer Methods in Electrical Engineering (Lecture Notes): Abhishek Yadav 1 18 cos( (0)) 2 (0) 120 0.1 max ,0 150 106 0.04 1 0.1 max399.90,0 150 106 1 0.1 399.90 150 106 2.6667 106 k 2 f t 0 h, v0 k1 h f 0 0.00002 ,0 2.6667 10 6 0.00002 f 0.00002 ,53 .32 1 18 cos( (0.00002)) 2 (53.32) 120 0.1 max ,0 150 106 0.04 1 0.1 max 933.013,0 150 106 1 0.1 0 150 106 666.67 1 1 v1 v0 k1 k 2 h 2 2 0 2.6667 10 6 666 .67 0.00002 1 1 2 2 0 1.3330 10 0.00002 6 26.660V i 1, t1 t 0 h 0 0.00002 0.00002 v1 26 .660 V k1 f t1 , v1 f 0.00002 ,26 .660 1 18 cos(120 (0.00002 )) 2 (26 .660 ) 0.1 max ,0 150 10 6 0.04 1 0.1 max 266.51,0 150 106 1 0.1 0 150 106 666.67 k 2 f t1 h, v1 k1 h f 0.00002 0.00002 ,26 .660 666 .67 0.00002 f 0.00004 ,26 .647 1 18 cos(120 (0.00004 )) 2 (26 .647 ) 0.1 max ,0 150 10 6 0.04 1 0.1 max 266.22,0 150 106 1 0.1 0 150 106 666.67 TEE-200: Computer Methods in Electrical Engineering (Lecture Notes): Abhishek Yadav 1 1 v 2 v1 k1 k 2 h 2 2 1 26 .660 666 .67 666 .67 0.00002 1 2 2 26 .660 666 .67 0.00002 26.647V v2 v0.00004 26 .647 V The results from Heun’s method are compared with exact results in Figure 2. The solution to this nonlinear equation is v(0.00004) 15.974V Voltage, v(V) h=0.00004 h=0.00002 Exact h=0.00001 0 0 0.00001 0.00002 0.00003 0.00004 Time, t(sec) Figure 2. Heun’s method results for different step sizes Using smaller step size would increases the accuracy of the result as given in Table 1 and Figure 3 below. Table 1. Effect of step size for Heun’s method Step size, h v0.00004 Et t % 0.00004 53.307 -37.333 233.71 0.00002 26.400 -10.426 65.269 0.00001 15.979 -0.0050000 0.031301 0.000005 15.917 0.057000 0.35683 0.0000025 15.968 0.0060000 0.037561 v(0.00004) Voltage, Step size, h 0 0 0.00001 0.00002 0.00003 0.00004 Figure 3. Effect of step size in Heun’s method TEE-200: Computer Methods in Electrical Engineering (Lecture Notes): Abhishek Yadav In Table 2, the Euler’s method and Runge-Kutta 2nd order method results are shown as a function of step size, Table 2. Comparison of Euler and the Runge-Kutta methods Step size, v(0.00004) h Euler Heun Midpoint Ralston 0.00004 106.64 53.307 -0.026667 35.529 0.00002 53.307 26.400 -0.026667 17.751 0.00001 26.640 15.979 11.642 15.363 0.000005 15.995 15.917 15.917 15.917 0.0000025 15.992 15.968 15.968 15.968 while in Figure 4, the comparison is shown over the range of time. 30 Euler 25 v(V) 20 Analytical Heun Voltage, 15 Ralston 10 Midpoint 5 0 0 0.00001 0.00002 0.00003 0.00004 Tim e, t (sec) Figure 4. Comparison of Euler and Runge Kutta methods with exact results over time. How do these three methods compare with results obtained if we found f x, y directly? Of course, we know that since we are including first three terms in the series, if the solution is a polynomial of order two or less (that is, quadratic, linear or constant), any of the three methods are exact. But for any other case the results will be different. Let us take the example of TEE-200: Computer Methods in Electrical Engineering (Lecture Notes): Abhishek Yadav e 2 x 3 y, y0 5 . dy dx If we directly find the f x, y , the first three terms of Taylor series gives yi 1 yi f xi , yi h f xi , yi h2 1 2! where f x, y e2 x 3 y f x, y 5e2 x 9 y For a step size of h 0.2 , using Heun’s method, we find y0.6 1.0930 The exact solution yx e2 x 4e3x gives y0.6 e20.6 4e30.6 0.96239 Then the absolute relative true error is 0.96239 1.0930 t 100 0.96239 13.571% For the same problem, the results from the Euler and the three Runge-Kutta method are given below Table 3. Comparison of Euler’s and Runge-Kutta 2nd order methods y(0.6) Exact Euler Direct 2nd Heun Midpoint Ralston Value 0.96239 0.4955 1.0930 1.1012 1.0974 1.0994 t % 48.514 13.571 14.423 14.029 14.236 What is the Runge-Kutta 4th Order Method? Runge-Kutta 4th order method is a numerical technique to solve ordinary differential equation of the form f x, y , y0 y0 dy dx So only first order ordinary differential equations can be solved by using Runge-Kutta 4th order method. In other sections, we have discussed how Euler and Runge-Kutta methods are used to solve higher order ordinary differential equations or coupled (simultaneous) differential equations. How does one write a first order differential equation in the above form? Example 1: 2 y 1.3e x , y0 5 dy dx is rewritten as 1.3e x 2 y, y0 5 dy dx In this case TEE-200: Computer Methods in Electrical Engineering (Lecture Notes): Abhishek Yadav f x, y 1.3e x 2 y Another example is given as follows Example 2: x 2 y 2 2Sin(3x ), y 0 5 dy ey dx is rewritten as dy 2 Sin(3x ) x 2 y 2 , y 0 5 dx ey In this case 2 Sin(3x ) x 2 y 2 f x, y ey Runge-Kutta 4th order method is based on the following yi 1 yi a1k1 a2 k2 a3k3 a4 k4 h (1) where knowing the value of y yi at x i , we can find the value of y=yi+1 at x i 1 , and h xi 1 xi Equation (1) is equated to the first five terms of Taylor series 2 3 dy xi 1 xi 1 d y xi , yi xi 1 xi 2 1 d y xi , yi xi 1 xi 3 yi 1 yi xi , y i 2 3 dx 2! dx 3! dx 4 (2) x , y xi 1 xi 1d y 4 4! dx 4 i i f x, y and xi 1 xi h dy Knowing that dx yi 1 yi f xi , yi h f ' xi , yi h 2 f '' xi , yi h3 f ''' xi , yi h 4 1 1 1 (3) 2! 3! 4! Based on equating Equation (2) and Equation (3), one of the popular solutions used is yi 1 yi k1 2k2 2k3 k4 h 1 (4) 6 k1 f xi , yi (5a) 1 1 k2 f xi h, yi k1h (5b) 2 2 1 1 k3 f xi h, yi k2 h (5c) 2 2 k4 f xi h, yi k3h (5d) Example 3: A rectifier-based power supply requires a capacitor to temporarily store power when the rectified waveform from the AC source drops below the target voltage. To properly size this capacitor a first- order ordinary differential equation must be solved. For a particular power supply, with a capacitor of 150 μF, the ordinary differential equation to be solved is dv(t ) 1 18 cos(120 (t )) 2 v(t ) 0.1 max ,0 , v(0) 0 dt 150 10 6 0.04 TEE-200: Computer Methods in Electrical Engineering (Lecture Notes): Abhishek Yadav Find voltage across the capacitor at t= 0.00004s. Use step size h=0.00002 Solution dv 1 18 cos(120 (t )) 2 v 0.1 max ,0 , dt 150 10 6 0.04 18 cos(120 (t )) 2 v f t , v 1 0.1 max ,0 150 10 6 0.04 vi 1 vi 1 k1 2k 2 2k 3 k 4 h 6 For i 0 , t 0 0 , v0 0V k1 f t 0 , v0 f 0,0 1 18 cos(120 (0)) 2 (0) 0.1 max ,0 150 10 6 0.04 1 0.1 max400,0 150 106 1 0.1 400 150 106 2.6660 106 1 1 k 2 f t 0 h, v0 k1 h 2 2 f 0 0.00002 ,0 2.6660 10 6 0.00002 1 1 2 2 f 0.00001 ,26 .660 1 18 cos( (0.00001)) 2 ( 26.660) 120 0.1 max ,0 150 106 0.04 1 0.1 max 266.50,0 150 106 1 0.1 266.50 150 106 666.67 1 1 k 3 f t 0 h, v 0 k 2 h 2 2 f 0 0.00002 ,0 666 .67 0.00002 1 1 2 2 f 0.00001 ,0.0066667 1 18 cos(120 (0.00001 )) 2 (0.0066667 ) 0.1 max ,0 150 10 6 0.04 1 0.1 max400.16,0 150 106 1 0.1 400.16 150 106 2.6671 106 TEE-200: Computer Methods in Electrical Engineering (Lecture Notes): Abhishek Yadav 1 1 k 4 f t 0 h, v 0 k 3 h 2 2 f 0 0.00002 ,0 2.6671 10 6 0.00002 1 1 2 2 f 0.00001 ,26 .671 1 18 cos( (0.00001)) 2 (26.671) 120 0.1 max ,0 150 106 0.04 1 0.1 max 10.671,0 150 106 1 0.1 0 150 106 666 .67 1 v1 v0 (k1 2k 2 2k 3 k 4 )h 6 0 2.6660 106 2 666.67 22.6671 106 666.670.00002 1 6 0 7.9922 106 0.00002 1 6 26.641 V v1 is the approximate voltage at t t1 t0 h 0 0.00002 0.00002 v1 v 0.00002 26 .641V For i 1, t1 0.00002 , v1 26 .641V k1 f t1 , v1 f 0.00002 ,26 .641 1 18 cos( (0.00002)) 2 ( 26.641) 120 0.1 max ,0 150 106 0.04 1 0.1 max 266.04,0 150 106 1 0.1 0 150 106 666.67 1 1 k 2 f t1 h, v1 k1 h 2 2 f 0.00002 0.00002 ,26 .641 666 .67 0.00002 1 1 2 2 f 0.00003 ,26 .634 1 18 cos( (0.00003)) 2 ( 26.634) 120 0.1 max ,0 150 106 0.04 1 0.1 max 265.88,0 150 106 1 0.1 0 150 106 666 .67 TEE-200: Computer Methods in Electrical Engineering (Lecture Notes): Abhishek Yadav 1 1 k 3 f t1 h, v1 k 2 h 2 2 f 0.00002 0.00002 ,26 .641 666 .67 0.00002 1 1 2 2 f 0.00003 ,26 .634 1 18 cos( (0.00003)) 2 (26.634) 120 0.1 max ,0 150 106 0.04 1 0.1 max 265.87,0 150 106 1 0.1 0 150 106 666 .67 k 4 f t1 h, v1 k 3 h f 0.00002 0.00002 ,26 .641 666 .67 0.00002 1 1 2 2 f 0.00003 ,26 .634 1 18 cos( (0.00003)) 2 (26.634) 120 0.1 max ,0 150 106 0.04 1 0.1 max 265.87,0 150 106 1 0.1 0 150 106 666 .67 1 v2 v1 (k1 2k 2 2k 3 k 4 )h 6 26.641 666.67 2 666.67 2 666.67 666.670.00002 1 6 26.641 4000.0 0.00002 1 6 26 .627V v 2 is the approximate voltage at t t2 t2 t1 h = 0.00002 0.00002 0.00004 v2 v.00004 26 .627V Figure 1 compares the exact solution with the numerical solution using Runge-Kutta 4th order method step size of h=0.00004. TEE-200: Computer Methods in Electrical Engineering (Lecture Notes): Abhishek Yadav v(V) h=0.00004 Voltage, h=0.00002 h=0.00001 Exact 0 0 0.00001 0.00002 0.00003 0.00004 Time,t(sec) Figure 1. Comparison of Runge-Kutta 4th order method with exact solution for different step sizes Table 1 and Figure 2 shows the effect of step size on the value of the calculated temperature at t=0.00004 seconds. Table 1. Value of voltage at time, t=0.00004s for different step sizes Step size, h v0.00004 Et |t | % 0.00004 53.335 -37.361 233.89 0.00002 26.627 -10.673 66.815 0.00001 15.985 -0.011000 0.068862 0.000005 15.973 0.0010000 0.0062600 0.0000025 15.974 0 0 v(0.00004) Voltage, 10 0 0.00001 0.00002 0.00003 0.00004 Step size, h Figure 2. Effect of step size in Runge-Kutta 4th order method In Figure 3, we are comparing the exact results with Euler’s method (Runge-Kutta 1st order method), Heun’s method (Runge-Kutta 2nd order method) and Runge-Kutta 4th order method. TEE-200: Computer Methods in Electrical Engineering (Lecture Notes): Abhishek Yadav 60 v(V) 50 Euler Voltage, 40 4th order Heun 30 20 10 Exact 0 0 0.00001 0.00002 0.00003 0.00004 Time, t(sec) Figure 3. Comparison of Runge-Kutta methods of 1st, 2nd, and 4th order. The formula described in this chapter was developed by Runge. This formula is same as Simpson’s rd rule, if f x, y is only a function of x . There are other versions of the fourth order method just 1 3 like there are several versions of the second order methods. Formula developed by Kutta is yi 1 yi k1 3k 2 3k 3 k 4 h 1 (6) 8 where k1 f xi , yi (7a) 1 1 k 2 f xi h, yi hk1 (7b) 3 3 2 1 k 3 f xi h, yi hk1 hk 2 (7c) 3 3 k 4 f xi h, yi hk1 hk 2 hk3 (7d) This formula is the same as the Simpson’s th rule, if f x, y is only a function of x . 3 8 25.Algorithm for solution of integration using Trapezoidal rule What is the Trapezoidal Rule? Trapezoidal rule is based on the Newton-Cotes formula that if one approximates the integrand by an nth order polynomial, then the integral of the function is approximated by the integral of that nth order polynomial. Integrating polynomials is simple and is based on the calculus formula. b b n1 a n1 x dx n 1 , n 0 n (1) a So if we want to approximate the integral b I f ( x)dx (2) a to find the value of the above integral, one assumes TEE-200: Computer Methods in Electrical Engineering (Lecture Notes): Abhishek Yadav f ( x) f n ( x) (3) where f n ( x) a 0 a1 x ... a n 1 x n 1 a n x n . (4) where f n (x ) is an n th order polynomial. Trapezoidal rule assumes n 1 , that is, the area under the linear polynomial (straight line), b b a f ( x)dx f1 ( x)dx a Derivation of the Trapezoidal Rule Method 1: Derived from Calculus Hence b b f ( x)dx f ( x)dx a a 1 b (a0 a1 x)dx a b2 a2 2 . a0 (b a) a1 (5) But what is a0 and a1? Now if one chooses, (a, f (a)) and (b, f (b)) as the two points to approximate f (x) by a straight line from a to b , f (a) f1 (a) a0 a1 a (6) f (b) f1 (b) a 0 a1b (7) Solving the above two equations for a and b , f (b) f (a) a1 ba f (a)b f (b)a a0 (8) ba Hence from Equation (5), f (a)b f (b)a f (b) f (a) b 2 a 2 b a f ( x)dx ba (b a) ba 2 f (a ) f (b) (b a ) (9) 2 Method 2: Also derived from Calculus f 1 ( x ) can also be approximated by using the Newton’s divided difference polynomial as f (b) f (a) f 1 ( x) f (a) ( x a) (10) ba Hence b b f ( x)dx f ( x)dx a a 1 TEE-200: Computer Methods in Electrical Engineering (Lecture Notes): Abhishek Yadav f (b) f (a) b f (a) ( x a) dx a ba b f (b) f (a ) x 2 f (a) x ba 2 ax a f (b) f (a) b 2 a2 f (a)b f (a)a ab 2 a2 ba 2 f (b) f (a) b a2 2 f (a)b f (a)a ab ba 2 2 f (b) f (a ) 1 f (a)b f (a)a b a 2 ba 2 f (a)b f (a)a f (b) f (a) b a 1 2 1 1 1 1 f (a)b f (a)a f (b)b f (b)a f (a)b f (a)a 2 2 2 2 1 1 1 1 f (a)b f (a)a f (b)b f (b)a 2 2 2 2 f (a ) f (b) (b a ) (11) 2 This gives the same result as Equation (10) because they are just different forms of writing the same polynomial. Method 3: Derived from Geometry The Trapezoidal rule can also be derived from geometry. Look at Figure 2. The area under the curve f1(x) is a trapezoid. The integral b f ( x)dx Area of a trapezoid 1 ( Sum of parallel sides)(height) 2 f (b) f (a) (b a) 1 2 f (a ) f (b) (b a ) (12) 2 TEE-200: Computer Methods in Electrical Engineering (Lecture Notes): Abhishek Yadav Figure 2: Geometric Representation of Trapezoidal Rule Method 4: Derived from Method of Coefficients Trapezoidal rule can also be derived by the method of coefficients. The formula ba ba b f ( x)dx 2 f (a) 2 f (b) a (13) 2 ci f ( xi ) i 1 where ba c1 2 ba c2 2 x1 a x2 b Figure 3: Area by Method of Coefficients The interpretation is that f (x) is evaluated at points a and b , and each function evaluation is given ba a weight of . Geometrically, Equation (12) is looked at as an area of a trapezoid, while 2 Equation (13) is viewed as the sum of the area of two rectangles, as shown in Figure 3. How can one derive trapezoidal rule by the method of coefficients? Assume TEE-200: Computer Methods in Electrical Engineering (Lecture Notes): Abhishek Yadav b f ( x)dx c a 1 f (a) c2 f (b) (14) b b Let the right hand side be an exact expression for integrals of 1dx and xdx , that is, the formula a a will then also be exact for linear combinations of f ( x) 1 and f ( x) x , that is, f ( x) a 0 (1) a1 ( x) . b 1dx b a c a 1 c2 (15) b2 a2 b xdx a 2 c1a c2 b (16) Solving the above two equations gives ba c1 2 ba c2 (17) 2 Hence ba ba b f ( x)dx 2 f (a) 2 f (b) . a (18) Method 5: Another approach on the Method of Coefficients Trapezoidal rule can also be derived by the method of coefficients by another approach ba ba b f ( x)dx 2 f (a) 2 f (b) a Assume b f ( x)dx c a 1 f (a) c2 f (b) (19) Let the right hand side be exact for integrals of the form b a a 0 a1 x dx So b b x2 a0 a1 x dx a a 0 x a1 2 a b2 a2 a0 b a a1 2 (20) But we want b a a 0 a1 x dx c1 f (a) c2 f (b) (21) to give the same result as Equation (20) for f ( x) a 0 a1 x . b a a 0 a1 x dx c1 a0 a1a c2 a0 a1b a 0 c1 c 2 a1 c1 a c 2 b (22) Hence from Equations (20) and (22), TEE-200: Computer Methods in Electrical Engineering (Lecture Notes): Abhishek Yadav b2 a2 2 a0 c1 c2 a1 c1a c2 b a0 b a a1 Since a 0 and a1 are arbitrary for a general straight line c1 c2 b a b2 a2 c1 a c 2 b (23) 2 Again, solving the above two equations (23) gives ba c1 2 ba c2 (24) 2 Therefore b f ( x)dx c a 1 f (a) c2 f (b) ba ba f (a) f (b) (25) 2 2 Example 1: The concentration of benzene at a critical location is given by c 1.75 erfc 0.6560 e 32.73erfc 5.758 where x erfcx e z dz 2 So in the above formula 0.6560 erfc0.6560 e z2 dz Since e z decays rapidly as z , we will approximate 2 0.6560 erfc0.6560 e z2 dz 5 a) Use single segment Trapezoidal rule to find the value of erfc(0.6560). b) Find the true error, Et for part (a). c) Find the absolute relative true error for part (a). Solution f (a ) f (b) a) I (b a ) , where 2 a5 b 0.6560 f ( z) e z 2 f (5) e 5 1.38881011 2 TEE-200: Computer Methods in Electrical Engineering (Lecture Notes): Abhishek Yadav f (0.6560) e 0.6560 0.65029 2 1.3888 10 11 0.65029 I (0.6560 5) 2 1.4124 b) The exact value of the above integral cannot be found. We assume the value obtained by adaptive numerical integration using Maple as the exact value for calculating the true error and relative true error. 0.6560 erfc0.6560 e dz z 2 5 0.31333 so the true error is Et True Value Approximate Value 0.31333 (1.4124) 1.0991 c) The absolute relative true error, t , would then be True Error t 100 True Value 0.31333 (1.4124 ) 100 0.31333 350.79% Multiple-segment Trapezoidal Rule: In Example 1, the true error using a single segment trapezoidal rule was large. We can divide the interval [5, 0.6560] into [5, 2.828] and [2.828, 0.6560] intervals and apply Trapezoidal rule over each segment. f ( z) e z 2 0.6560 2.828 0.6560 5 f ( z )dz 5 f ( z )dz f ( z)dz 2.828 f (5) f (2.828 ) f (2.828 ) f (0.6560 ) (2.828 5) (0.6560 2.828 ) 2 2 11 f (5) 1.3888 10 f (2.828) 0.000336 f (0.6560) 0.65029 Hence 0.6560 1.3888 10 11 0.000336 0.000336 0.65029 5 f ( z )dz (2.828 5) 2 (0.6560 2.828) 2 0.70695 The true error, E t TEE-200: Computer Methods in Electrical Engineering (Lecture Notes): Abhishek Yadav Et 0.31333 (0.70695 ) 0.39362 The true error now is reduced from 1.0991 to 0.39362. Extending this procedure to dividing [a, b] into n equal segments and apply the Trapezoidal rule over each segment, the sum of the results obtained for each segment is the approximate value of the integral. Example 2: The concentration of benzene at a critical location is given by c 1.75 erfc 0.6560 e 32.73erfc 5.758 where x erfcx e z dz 2 So in the above formula 0.6560 erfc0.6560 e z2 dz z2 Since e decays rapidly as z , we will approximate 0.6560 erfc0.6560 e z2 dz 5 d) Use two segment Trapezoidal rule to find the value of erfc(0.6560). e) Find the true error, Et for part (a). f) Find the absolute relative true error for part (a). Solution: a) The solution using 2-segment Trapezoidal rule is ba n1 I f (a) 2 f (a ih) f (b) 2n i 1 n2 a5 b 0.6560 ba h n 0.6560 5 2 2.172 0.6560 5 21 I f (5) 2 f (a ih) f (0.6560) 2(2) i 1 4.344 f (5) 2 f (2.828) f (0.6560) 4 4.344 4 1.3888 10 11 2(0.000336) 0.65029 TEE-200: Computer Methods in Electrical Engineering (Lecture Notes): Abhishek Yadav 0.70695 b) The exact value of the above integral cannot be found. We assume the value obtained by adaptive numerical integration using Maple as the exact value for calculating the true error and relative true error. 0.6560 erfc0.6560 e dz z 2 5 0.31333 so the true error is Et True Value Approximate Value 0.31333 (0.70695) 0.39362 c) The absolute relative true error, t , would then be True Error t 100 True Value 0.31333 (0.70695 ) 100 0.31333 125.63% 0.6560 Table 1: Values obtained using multiple-segment Trapezoidal rule for erfc0.6560 e z2 dz 5 n Value Et t % a % 1 -1.4124 1.0991 350.79 --- 2 -0.70695 0.39362 125.63 99.793 3 -0.48812 0.17479 55.787 44.829 4 -0.40571 0.092379 29.483 20.314 5 -0.37028 0.056957 18.178 9.5662 6 -0.35212 0.038791 12.380 5.1591 7 -0.34151 0.028182 8.9946 3.1063 8 -0.33475 0.021426 6.8383 2.0183 Example 3: Use Multiple-segment Trapezoidal Rule to find the area under the curve 300x f ( x) 1 ex from x 0 to x 10 . Solution Using two segments, we get TEE-200: Computer Methods in Electrical Engineering (Lecture Notes): Abhishek Yadav 10 0 h 5 2 300(0) f (0) 0 1 e0 300(5) f (5) 10.039 1 e5 300(10) f (10) 0.136 1 e10 ba n1 I f (a) 2 f (a ih) f (b) 2n i 1 10 0 21 f (0) 2 f (0 5) f (10) 2(2) i 1 f (0) 2 f (5) f (10) 10 4 0 2(10.039) 0.136 50.535 10 4 So what is the true value of this integral? 10 300x 1 e x dx 246.59 0 Making the absolute relative true error 246 .59 50 .535 t 100 % 246 .59 79.506% Why is the true value so far away from the approximate values? Just take a look at Figure 5. As you can see, the area under the “trapezoids” (yeah, they really look like triangles now) covers a small the area under the curve. As we add more segments, the approximated value quickly approaches the true value. Figure 5: 2-Segment Approx. of an Exp. Function 10 300 x Table 2: Values obtained using Multiple-segment Trapezoidal Rule for 1 e 0 x dx TEE-200: Computer Methods in Electrical Engineering (Lecture Notes): Abhishek Yadav n Approximate Value Et t 1 0.681 245.91 99.724% 2 50.535 196.05 79.505% 4 170.61 75.978 30.812% 8 227.04 19.546 7.927% 16 241.70 4.887 1.982% 32 245.37 1.222 0.495% 64 246.28 0.305 0.124% Example 4: Use multiple-segment Trapezoidal Rule to find 2 1 I dx 0 x Solution We cannot use Trapezoidal Rule for this integral, as the value of the integrand at x 0 is infinite. However, it is known that a discontinuity in a curve will not change the area under it. We can assume any value for the function at x 0 . The algorithm to define the function so that we can use multiple-segment Trapezoidal Rule is given below. Function f(x) If x = 0 Then f = 0 If x <> 0 Then f = x ^ (-0.5) End Function Basically, we are just assigning the function a value of zero at x 0 . Everywhere else, the function is continuous. This means the true value of our integral will be just that—true. Let’s see what happens using multiple-segment Trapezoidal Rule: Using two segments, we get 20 h 1 2 f (0) 0 1 f (1) 1 1 1 f ( 2) 0.70711 2 ba n1 I f (a) 2 f (a ih) f (b) 2n i 1 TEE-200: Computer Methods in Electrical Engineering (Lecture Notes): Abhishek Yadav 20 21 f (0) 2 f (0 1) f (2) 2(2) i 1 f (0) 2 f (1) f (2) 2 4 0 2(1) 0.70711 2 4 1.3536 So what is the true value of this integral? 2 1 x dx 2.8284 0 Thus making the absolute relative true error 2.8284 1.3536 t 100 % 2.8284 52.143% 2 1 Table 3: Values obtained using Multiple-segment Trapezoidal Rule for 0 x dx Number of Approximate Value Et t Segments 2 1.354 1.474 52.14% 4 1.792 1.036 36.64% 8 2.097 0.731 25.85% 16 2.312 0.516 18.26% 32 2.463 0.365 12.91% 64 2.570 0.258 9.128% 128 2.646 0.182 6.454% 256 2.699 0.129 4.564% 512 2.737 0.091 3.227% 1024 2.764 0.064 2.282% 2048 2.783 0.045 1.613% 4096 2.796 0.032 1.141% Error in Multiple-segment Trapezoidal Rule The true error for a single segment Trapezoidal rule is given by (b a) 3 Et f " ( ), a b 12 where is some point in a, b . TEE-200: Computer Methods in Electrical Engineering (Lecture Notes): Abhishek Yadav What is the error, then in the multiple-segment Trapezoidal rule? It will be simply the sum of the errors from each segment, where the error in each segment is that of the single segment Trapezoidal rule. The error in each segment is E1 (a h) a3 f " ( ), a a h 1 1 12 h3 f " ( 1 ) 12 E2 (a 2h) (a h)3 f " ( ), a h a 2h 2 2 12 h3 f " ( 2 ) 12 . . . Ei (a ih) (a (i 1)h)3 f " ( ), a (i 1)h a ih i i 12 h3 f " ( i ) 12 . . . E n 1 a (n 1)h a (n 2)h3 f " ( ), a (n 2)h a (n 1)h n 1 n 1 12 h3 f " ( n 1 ) 12 En b a (n 1)h3 f " ( ), a (n 1)h b n n 12 h3 f " ( n ) 12 Hence the total error in multiple-segment Trapezoidal rule is n Et Ei i 1 h3 n 12 f " ( i 1 i ) (b a ) 3 n 12 n 3 f " ( i 1 i ) n (b a) 3 f " ( i ) i 1 12n 2 n n f " ( i 1 i ) The term is an approximate average value of the second derivative f " ( x), a x b . n TEE-200: Computer Methods in Electrical Engineering (Lecture Notes): Abhishek Yadav Hence n f " ( i ) (b a) 3 i 1 Et 12n 2 n Below is given the table for the integral 30 140000 2000 ln 140000 2100t 9.8t dt 8 as a function of the number of segments. You can visualize that as the number of segments are doubled, the true error gets approximately quartered Table 4: Values obtained using Multiple-segment Trapezoidal Rule for 30 140000 x 2000 ln 9.8t dt 8 140000 2100t n Value Et t % a % 2 11266 -205 1.854 5.343 4 11113 -51.5 0.4655 0.3594 8 11074 -12.9 0.1165 0.03560 16 11065 -3.22 0.02913 0.00401 For example, for 2-segment Trapezoidal rule, the true error is 205, and a quarter of that error is 51.25. That is close to the true error of 51.5 for the 4-segment Trapezoidal rule. Can you answer the question why is the true error not exactly 51.25? How does this information help us in numerical integration? You will find out that this forms the basis of Romberg integration based on Trapezoidal Rule, where we use the argument that true error gets approximately quartered when the number of segments is doubled. Romberg integration based on Trapezoidal rule is computationally more efficient than using Trapezoidal rule by itself in developing an adaptive integration scheme. TEE-200: Computer Methods in Electrical Engineering (Lecture Notes): Abhishek Yadav