# Matlab by elafanto

VIEWS: 53 PAGES: 48

• pg 1
```									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

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?"
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:
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

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:

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:

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,...

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 , y0  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 , y0  5
dy
dx
is rewritten as
 1.3e  x  2 y, y0  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,00.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  v0.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  v0.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,        v0.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 , y0  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 , y0  5
dy
dx
is rewritten as
 1.3e  x  2 y, y0  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   e2 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

 5e2 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  106   
                     0.04              

1
 0.1  max399.90,0
150  106

1
 0.1  399.90
150  106
 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  106   
                          0.04                   

1
 0.1  max 933.013,0
150  106

1
 0.1  0
150  106
 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  106

1
 0.1  0
150  106
 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  106

1
 0.1  0
150  106
 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  v0.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 v0.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, y0  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   e2 x  3 y
f x, y   5e2 x  9 y
For a step size of h  0.2 , using Heun’s method, we find
y0.6  1.0930
The exact solution
yx  e2 x  4e3x
gives
y0.6  e20.6  4e30.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 , y0  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 , y0  5
dy
dx
is rewritten as
 1.3e  x  2 y, y0  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  max400,0
150  106

1
 0.1  400
150  106
 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  106     
                          0.04                     

1
 0.1  max 266.50,0
150  106

1
 0.1  266.50
150  106
 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  max400.16,0
150  106

1
 0.1  400.16
150  106
 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  106   
                          0.04                    

1
 0.1  max 10.671,0
150  106

1
 0.1  0
150  106
 666 .67
1
v1  v0  (k1  2k 2  2k 3  k 4 )h
6
 0  2.6660  106  2 666.67  22.6671  106    666.670.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  106   
                          0.04                     

1
 0.1  max 266.04,0
150  106

1
 0.1  0
150  106
 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  106   
                          0.04                     

1
 0.1  max 265.88,0
150  106

1
 0.1  0
150  106

 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  106   
                          0.04                    

1
 0.1  max 265.87,0
150  106

1
 0.1  0
150  106
 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  106   
                          0.04                    

1
 0.1  max 265.87,0
150  106

1
 0.1  0
150  106
 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.670.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      v0.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 n1  a n1 
 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 
ba
f (a)b  f (b)a
a0                                                                                (8)
ba
Hence from Equation (5),
f (a)b  f (b)a           f (b)  f (a) b 2  a 2
b


a
f ( x)dx 
ba
(b  a) 
ba           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)
ba
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
ba               
b
           f (b)  f (a )  x 2     
  f (a) x 
ba

 2  ax 

                                    a
 f (b)  f (a)  b              
2
a2
 f (a)b  f (a)a                     ab 
 2          a2 

     ba               2       
 f (b)  f (a)  b          a2 
2
 f (a)b  f (a)a                         ab  

      ba        2         2 
 f (b)  f (a )  1
 f (a)b  f (a)a                       b  a 
2

      ba       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
ba          ba
b

 f ( x)dx  2 f (a)  2 f (b)
a
(13)
2
  ci f ( xi )
i 1

where

ba
c1 
2
ba
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
ba
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
ba
c1 
2
ba
c2                                                                               (17)
2
Hence
ba         ba
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
ba          ba
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
ba
c1 
2
ba
c2                                                                                (24)
2
Therefore
b

 f ( x)dx  c
a
1   f (a)  c2 f (b)

ba         ba
            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
erfcx    e  z dz
2



So in the above formula
0.6560
erfc0.6560                   e
z2
dz


Since e  z decays rapidly as z   , we will approximate
2

0.6560
erfc0.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       
a5
b  0.6560
f ( z)  e  z
2

f (5)  e 5  1.38881011
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
erfc0.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
erfcx    e  z dz
2



So in the above formula
0.6560
erfc0.6560              e
z2
dz

 z2
Since e          decays rapidly as z   , we will approximate
0.6560
erfc0.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
ba             n1                
I        f (a)  2 f (a  ih)  f (b)
2n             i 1               
n2
a5
b  0.6560
ba
h
n
0.6560  5

2
 2.172
0.6560  5              21                     
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
erfc0.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                erfc0.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
ba               n1                
I          f (a)  2 f (a  ih)  f (b)
2n               i 1               
10  0            21                
          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
20
h        1
2
f (0)  0
1
f (1)      1
1
1
f ( 2)        0.70711
2
ba              n1                
I         f (a)  2 f (a  ih)  f (b)
2n              i 1               

TEE-200: Computer Methods in Electrical Engineering (Lecture Notes): Abhishek Yadav
20              21               
        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)  a3 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)h3 f " ( ), a  (n  2)h    a  (n  1)h
n 1                   n 1
12
h3
      f " ( n 1 )
12

En 
b  a  (n  1)h3 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

```
To top