OTDM PRACTICAL

ASSIGNMENT USING
SCILAB
SEMESTER- 4th
B.TECH ECE-AI-1

SUBMITTED TO: NAME: JULIE GAYARI

DR. JYOTI SHOKEEN ROLL NO.- 04801182023
Contents
1 Introduction to SCILAB
1.1 Overview
1.2 Become Familiar With Scilab
1.3 The Menu Bar
1.4 Pre-defined mathematical variables and operators
1.5 Booleans
1.6 Complex Numbers
1.7 Strings
1.8 Matrices
1.9 Looping and Branching
1.10 Programming in Scilab
1.10.1 Editor
1.11 Functions
1.12 The Graphic Window

2 Lab Objective and Requirements

2.1 Lab Objectives
2.2 Course Outcome:
2.3 Lab Requirements
2.4 List of Experiments
 Matrix Operations
2.5 Matrix Addition
2.6 Matrix Multiplication
2.7 Matrix Transpose
2.8 Inverse of 3X3 Matrix using Gauss Jordan Method
2.9 Eigen Values & Eigen Vectors of 2 X 2 Matrix

3 Measures of Central Tendency

3.1 Measures of Central Tendency (Ungrouped Data)
3.2 Measures of Central Tendency (Grouped Data)

4 Plotting of 2D Graphs Using Scilab

4.1 A simple graph plotting raw data
4.2 Generation of square wave
4.3 Unit Step function I
4.4 Unit Step function II

5 Numerical Methods for Finding Roots of Algebraic and

Transcendental Equations 6.1 Bisection Method
6.2 Newton-Raphson Method

6 Numerical Integration for Solving Definite Integrals

6.1 The Trapezoidal Rule
6.2 Simpson’s 1/3rd Rule
6.3 Simpson’s 3/8th Rule

7 Numerical Solution of Ordinary Differential equations using

Euler’s Method
7.1 Euler’s Method

8 Numerical Solution of Ordinary Differential Equations

Using Runge- Kutta Method
8.1 Runge-Kutta method of 4th Order 2
Chapter 1:Introduction to Scilab

Scilab is a programming language associated with a rich collection of numerical

algorithms covering many aspects of scientific computing problems.
From the software point of view, Scilab is an interpreted language. This generally
allows to get faster development processes, because the user directly accesses to a
high level language, with a rich set of features provided by the library. The Scilab
language is meant to be extended so that user-defined data types can be defined with
possibly overloaded operations. Scilab users can develop their own module so that
they can solve their particular problems. The Scilab language allows to dynamically
compile and link other languages such as Fortran and C and in this way, external
libraries can be used as if they were a part of Scilab built-in features.
Scilab is also a numerical computation software that anybody can freely download.
Available under Windows, Linux and Mac OS X, Scilab can be downloaded at the
following address http://www.scilab.org/.
An online help is provided in many local languages.
From a scientific point of view, Scilab comes with many features. At the very
beginning of Scilab, features were focused on linear algebra. But rapidly, the number
of features extended to cover many areas of scientific computing such as matrices,
statistics, Ordinary differential equations, signal processing, interpolation,
approximation, linear, quadratic and non linear optimization and many more.

1.2 Become Familiar With Scilab

The Scilab workspace consists of several windows:
• The console for making calculations
• The editor for writing programs
• The graphics windows for displaying graphics
• The embedded help

1
1.3 The Menu Bar
The following options in menu are particularly useful:
Applications
• The command history allows you to find all the commands from previous
sessions to the current session.

• The variables browser allows you to find all variables previously used during
the current session.
Editing

• Preferences (in Scilab menu under Mac OS X) allow you to set and customize
colors, fonts and font size in the console and in the editor, which is very useful
for screen projection.

• Clicking on Clear Console clears the entire content of the console. In this case,
the command history is still available and calculations made during the session
remain in memory. Commands that have been erased are still available through
the keyboards arrow keys.

1.4 Pre-defined mathematical variables and operators

In Scilab, several mathematical variables are pre-defined variables, which name
begins with a percent % character. The variables which have a mathematical meaning
are summarized in the given table:

Also apart from usual operators for summation, subtraction, multiplication and
division, comparison operators are as given:

2
1.5 Booleans
Boolean variables can store true or false values. In Scilab, true is written with %t
or %T and false is written with %f or %F.

1.6 Complex Numbers

Scilab provides complex numbers, which are stored as pairs of floating point
numbers. The predefined variable %i represents the mathematical imaginary number i
which satisfies i2 = —1.

1.7 Strings
Strings can be stored in variables, provided that they are delimited by double
quotes. The concatenation operation is available from the + operator. They are many
functions which allow to process strings, including regular expressions.

1.8 Matrices
There is a simple and efficient syntax to create a matrix with given values.
The following is the list of symbols used to define a matrix:
• square brackets ”[” and ”]” mark the beginning and the end of the matrix

• commas ”,” separate the values on different columns

• semicolons ”;” separate the values of different rows

The following syntax can be used to define an m × n matrix, where blank spaces are
optional (but make the line easier to read) and ”...” are designing intermediate values:
A = [a 11, a 12, ..., a 1n; a 21, a 22, ..., a 2n; ...; am 1, am 2, ..., amn]

3
Use of colon operator in matrices
Following table depicts use of colon operator in matrices:

Matrix and usual operators

Following table shows various operator used in matrix operations:

1.9 Looping and Branching

In this section, we describe how to make conditional statements
The if statement
The if uses a boolean variable to perform its choice: if the boolean is true, then the
statement is executed. A condition is closed when the end keyword is met. In the
following script, we display the string ”Hello!” if the condition %t, which is always
true, is satisfied.
if ( %t ) then disp
(” Hello ! ”) end
The previous script produces:
Hello !
If the condition is not satisfied, the else statement allows to perform an alternative
statement, as in the following script:
if ( %f ) then disp
(” Hello ! ”)
else

4
disp (” Goodbye !” ) end
The previous script produces:
Goodbye !
In order to get a boolean, any comparison operator can be used, e.g. ==, >, etc or any
function which returns a boolean. In the following session, we use the == operator to
display the message ”Hello !”.
i = 2 if ( i
== 2
) then
disp (”
Hello ! ”) else disp (”
Goodbye !” ) end
The select statement
The select statement allows to combine several branches in a clear and simple way.
Depending on the value of a variable, it allows to perform the statement
corresponding to the case keyword. There can be as many branches as required. In the
following script, we want to display a string which corresponds to the given integer i.
i=2
select i case 1
disp (” One ”)
case 2 disp (”
Two ”) case 3
disp (” Three ”)
else

disp (” Other ”) end

The previous script prints out ”‘Two”, as expected. The else branch is used if all the
previous case conditions are false.
The for statement
The for statement allows to perform loops, i.e. to perform a given action several
times. Most of the time, a loop is performed over an integer value, which goes from a
starting to an ending index value. We will see, at the end of this section, that the for

5
statement is in fact much more general, as it allows to loop through the values of a
matrix. In the following Scilab script, we display the value of i, from 1 to 3.
for i = 1 : 3
disp (i) end
The previous script produces the following output.
1.
2.
3.
In the previous example, the loop is performed over a matrix of floating point
numbers containing integer values. Indeed, we used the colon ”:” operator in order to
produce the vector of index values [1, 2, 3]. The following session shows that the
statement 1:3 produces all the required integer values into a row vector.
—— > i = 1 : 3 i
= 1.2.3.
The while statement
The while statement allows to perform a loop while a boolean expression is true. At
the beginning of the loop, if the expression is true, the statements in the body of the
loop are executed. When the expression becomes false (an event which must occur at
certain time), the loop is ended. In the following script, we compute the sum of the
numbers i from 1 to 10 with a while statement.

s=0i
=1
while ( i <= 10 )
s=s+ii
=i
+ 1 end
At the end of the algorithm, the values of the variables i and s are: s =
55 and i = 11.
The break and continue statement
The break statement allows to interrupt a loop. Usually, we use this statement in loops
where, if some condition is satisfied, the loops should not be continued further. In the

6
following example, we use the break statement in order to compute the sum of the
integers from 1 to 10. When the variable i is greater than 10, the loop is interrupted.
s=0i
=1
while ( %t ) if ( i
> 10 ) then break
end
s=s+ii
=i
+ 1 end
At the end of the algorithm, the values of the variables i and s are:
s = 55 and i = 11
The continue statement allows to go on to the next loop, so that the statements in the
body of the loop are not executed this time. When the continue statement is executed,
Scilab skips the other statements and goes directly to the while or for statement and
evaluates the next loop. In the following example, we compute the sum s = 1 + 3 + 5
+ 7 + 9 = 25. The modulo(i, 2) function returns 0 if the
number i is even. In this situation, the script increments the value of i and use the
continue statement to go on to the next loop.
s=0i
=1
while ( i <= 10) if ( modulo (i,
2) == 0 ) then
i=i+
1
continue
else
s=s+ii
=i

7
Chapter 2 Lab Objective and Requirement

2.1 Lab Objectives

The purpose of this course is to introduce the students with Scilab applications in
Mathematical operations. The students will be able to enhance their analyzing and
problem solving skills and use the same for writing programs in Scilab.

2.2 Course Outcome:

At the end of the course, the students will be able to:

C154.1: write, compile and debug programs in Scilab.

C154.2: use conditional expressions and looping statement to solve problems associated
with conditions and repetitions.

C154.3: implement the programs using arithmetic and relational operators.

C154.4: understand the concept of inbuilt and user defined functions.

C154.5: understand and solve matrices operations effectively.

C154.6: able to choose programming components that efficiently solve computing

problems in Mathematics .

2.3 Lab Requirements

Hardware Details :
Intel i3/C2D Processor/2 GB RAM/500GB HDD/MB/LAN Card/ Key Board/
Mouse/CD Drive/15” Color Monitor/ UPS Software
Detail :
Operating System: BSDs (e.g., FreeBSD), Linux, macOS,
Windows Open source software: Scilab 5.5.2

8
Chapter 3 Matrix Operations

3.1 Matrix Addition

• (A + B)i,j = Ai,j + Bi,j
• Order of the matrices must be the same
• Matrix addition is commutative
• Matrix addition is associative Matrix Addition Algorithm:

1. Start

2. Declare variables and initialize necessary variables

3. Enter the elements of matrices row wise using loops

4. Add the corresponding elements of the matrices using nested loops

5. Print the resultant matrix as console output

6. Stop
Flow Chart for Matrix Addition

9
 Scilab Code for Matrix Addition
clc
m=input(”enter number of rows in the Matrices: ”);
n=input(”enter number of columns in the Matrices: ”);
disp(’enter the first Matrix’) for i=1:m for j=1:n
A(i,j)=input(’\’); end end
disp(’enter the second Matrix’)
for i=1:m for j=1:n
B(i,j)=input(’\’); end end for i=1:m for
j=1:n C(i,j)=A(i,j)+B(i,j); end end
disp(’The first matrix is’) disp(A)
disp(’The Second matrix is’) disp(B)
disp(’The sum of the two matrices is’)
disp(C)

Matrix Addition using functions

// Save file as addition.sce clc
function [ ]=addition(m, n, A, B)
C=zeros(m,n);
C=A+B;
disp(’The first matrix is’) disp (A)
disp(’The Second matrix is’) disp
(B) disp(’The sum of two matrices
is’) disp (C) endfunction

10
3.2 Matrix Multiplication

• Am×n × Bn×p = ABm×p

Where (AB)ij = Ai 1 B 1 j + Ai2 B 2 j + . . . + AinBnj

• The number of columns in the first matrix must be equal to the number of rows in
the second matrix

• Matrix multiplication is not commutative Matrix Multiplication Algorithm:

1. Start

2. Declare variables and initialize necessary variables

3. Enter the elements of matrices row wise using loops

4. Check the number of rows and column of first and second matrices

5. If number of rows of first matrix is equal to the number of columns of second matrix, go
to step 6. Otherwise, print ’Matrices are not conformable for multiplication’ and go to
step 3

6. Multiply the matrices using nested loops

7. Print the product in matrix form as console output

8. Stop
Flow Chart for Matrix Multiplication

Scilab Code for Matrix Multiplication clc

11
m=input(”Enter number of rows in the first Matrix: ”);
n=input(”Enter number of columns in the first Matrix: ”);
p=input(”Enter number of rows in the second Matrix: ”);
q=input(”Enter number of columns in the second Matrix: ”); if
n==p disp(’Matrices are conformable for multiplication’) else
disp(’Matrices are not conformable for multiplication’) break;
end disp(’enter the first Matrix’) for i=1:m for j=1:n
A(i,j)=input(’\’); end end
disp(’enter the second Matrix’)
for i=1:p for j=1:q
B(i,j)=input(’\’);
end end
C=zeros(m,q); for
i=1:m for j=1:q
for k=1:n
C(i,j)=C(i,j)+A(i,k)*B(k,j);
end end
end disp(’The first matrix is’) disp(A)
disp(’The Second matrix is’) disp(B)
disp(’The product of the two matrices is’)
disp(C)

Matrix Multiplication using functions

// Save file as multiplication.sce clc
function [ ] = multiplication(m, n, p, q, A, B)
C=zeros(m,n); if n==p
disp(’Matrices are conformable for multiplication’)
else disp(’Matrices are not conformable for
multiplication’) break; end C=A*B
disp(’The first matrix is’) disp (A) disp(’The
Second matrix is’) disp (B) disp(’The
multiplication of two matrices is’) disp (C)
endfunction

12
3.3 Matrix Transpose
T
• The transpose of an m × n matrix A is n × m matrix A

• Formed by interchanging rows into columns and vice versa

• (AT )i,j = Aj,i

Matrix Transpose Algorithm:

1. Start

2. Declare variables and initialize necessary variables

3. Enter the elements of matrix by row wise using loop

4. Interchange rows to columns using nested loops

5. Print the transposed matrix as console output

6. Stop
Flow Chart for Matrix Transpose

13
 Scilab Code for Matrix Transpose
clc
m=input(”Enter number of rows in the Matrix: ”);
n=input(”Enter number of columns in the Matrix: ”);
disp(’Enter the Matrix’) for i=1:m for j=1:n
A(i,j)=input(’\’); end end
B=zeros(n,m); for i=1:n for
j=1:m B(i,j)=A(j,i) end end
disp(’Entered matrix is’)
disp(A) disp(’Transposed
matrix is’) disp(B)
Matrix Transpose using functions
// Save file as transpose.sce function
[ ]=transpose(m, n, A)
B=zeros(m,n);
B=A’
disp(’The matrix is’) disp (A)
disp(’Transposed matrix is’)
disp (B) endfunction

3.4 Inverse of a 3×3 Matrix Using Gauss Jordan Method

Matrix Inverse

• The inverse of a square matrix A denoted by A −1 is such that AA −1 = I

• Inverse of a square matrix exists if and only if matrix is non singular Matrix Inverse

Algorithm:

1. Start

2. Enter the elements of the 3 × 3 matrix row wise using loop

3. Check Det(A). If Det(A) = 0, Print ’Inverse does not exist’ else go to step 4.

4. Make an augmented matrix B = [A, I], where I is unit matrix.

14
 5. Make element B(1, 1) = 1 and using this make B(2, 1) = 0 and B(3, 1) = 0 using
row/column transformations.

6. Make element B(2, 2) = 1 and using this make B(1, 2) = 0 and B(3, 2) = 0.

7. Make element B(3, 3) = 1 and using this make B(1, 3) = 0 and B(2, 3) = 0.

8. Thus final augmented matrix is in the form B[I, A−1]. Print A−1

9. Stop
Flow Chart for Matrix Inverse

Scilab Code for Matrix Inverse using Gauss Jordan Method clc
disp(’Enter a 3 by 3 matrix row-wise, make sure that diagonal elements are non
-zeros’)
for i=1:3
for j=1:3
A(i,j)=input(’\’); end end disp(’Entered Matrix
is’) disp(A) if det(A)==0 disp(’Matrix is singular,
Inverse does not exist’) break; end
//Taking the augmented matrix
B=[A eye(3,3)]
disp(’Augumented matrix is:’)
disp(B)
// Making B(1,1)=1
B(1,:) = B(1,:)/B(1,1);
//Making B(2,1) and
15
B(3,1)=0 B(2,:) = B(2,:) -
B(2,1)*B(1,:);
B(3,:) = B(3,:) - B(3,1)*B(1,:);
//Making B(2,2)=1 and B(1,2),
B(3,2)=0 B(2,:) = B(2,:)/B(2,2);
B(1,:) = B(1,:) - B(1,2)*B(2,:);
B(3,:) = B(3,:) - B(3,2)*B(2,:);
//Making B(3,3)=1 and B(1,3),
B(2,3)=0 B(3,:) = B(3,:)/B(3,3);
B(1,:) = B(1,:) - B(1,3)*B(3,:);
B(2,:) = B(2,:) - B(2,3)*B(3,:); disp(’Augumented
matrix after row operations is:’) disp(B) B(:,1:3)=[ ]
disp(’Inverse of the Matrix is’) disp(B)

Matrix Inverse using functions

// Save file as transpose.sce

function [ ]=inverse(m, A)
C=zeros(m,m);
B=det(A) if
B==0
disp(’Matrix is singular, Inverse does not exist’)
break; end C=inv(A) disp(’The matrix is’) disp
(A) disp(’Inverse of given matrix is:’) disp (C)
endfunction

3.5 Eigen Values and Eigen Vectors of a 2×2 Matrix

Eigen Values and Eigen Vectors

• Eigenvalues are a special set of scalars associated with a n × n matrix.

• These are roots of the characteristic equation |A—λI = 0|, hence also known as
characteristic roots or latent roots.

• X is a column vector associated with an eigen value λ such that AX = λX

Algorithm to find eigen values and eigen vectors of a 2 × 2 matrix

16
1. Start

2. Enter the elements of the 2 × 2 matrix A,row wise using loop

3. Find Trace(A) and Det(A) as the 2 eigen values λ1, λ2 are roots of the char-
-acteristic equation λ2 — trace(A) + det(A) = 0

4. Find λ1 and λ2

5. Find corresponding eigen vectors X 1 and X 2 as follows:

If a 1,2 /= 0, calculate X 1 and X 2 from the first row
If a2,1 /= 0, calculate X1 and X2 from the second row
If a 1 , 2 = 0 and a 2 , 1 = 0,
1
X1 =
0

0
X2 =
1

6. Stop
Flow Chart for Eigen Values and Vectors

Scilab Code for Eigen Values of a 2 × 2 Matrix

17
disp(’Enter the 2 by 2 Matrix row-wise’)
for i=1:2 for j=1:2
A(i,j)=input(’\’); end end
b=A(1,1)+A(2,2);
c=A(1,1)*A(2,2)-A(1,2)*A(2,1);
// characteristic equation is λ2 — trace(A) + det(A) = 0, here λ1 ≡ e1, λ2 ≡ e2 e1
= (b + sqrt(b∧2 — 4 ∗ c))/2; e2 =
(b — sqrt(b∧2 — 4 ∗ c))/2; if
A(1, 2)=˜ 0
X1 = [A(1,2); e1-A(1,1)]; X2
= [A(1,2); e2-A(1,1)]; elseif
A(2, 1)=˜ 0
X1 = [e1-A(2,2); A(2,1)];
X2 = [e2-A(2,2); A(2,1)];
else
X1 = [1; 0];
X2 = [0; 1];
end disp(’First Eigen value is:’);
disp(e1) disp(’First Eigen vector
is:’); disp (X1) disp(’Second
Eigen value is:’); disp(e2)
disp(’Second Eigen vector is:’);
disp (X2)

Chapter 4 Measures of Central Tendency

18
4.1 Measures of Central Tendency (Ungrouped Data)

Algorithm to find mean, median, mode and moments of ungrouped data

1. Start

2. Arrange the data in ascending order

3. Find number of observations (n)

4. Evaluate Mean = sum of observationsn

5. Compute Median = Value of (n+21) th observation, if n is odd, = Average value of nth and (n +
1)th observations, if n is even 2 2

6. For finding mode of given data:

I. Arrange the data in ascending order
II. Find the difference between the adjacent elements
III. Find the indices at which the value is non zero
IV. Locate the position where there is largest gap between the non zero indices
V. Term at next position gives mode
7. Find rth moment about mean x 8. Evaluate S.D.
µ
Flow Chart for Central Tendencies (Ungrouped Data)

//Scilab Code to find mean, mode, median, moments, skewness and kurtosis of
linear data clc

19
function [ ]= moments(A) B=gsort(A); n
= length(B); meanA = sum(B)/n; if
pmodulo(n,2)==0 medianA
=((B(n/2)+B(n/2 +1)))/2; else medianA
= B((n+1)/2); end
C = diff(B)
//C= diff(B) calculates differences between adjacent elements of B along the first array
dimension whose size exceeds 1:
//If B is a vector of length n, then C = diff(B) returns a vector of length n-1.
The elements of C are the differences between adjacent elements of B.
//C = [B(2)-B(1) B(3)-B(2)........B(m)-B(m-1)]
D = find(C)
//D = find(C) finds the idices(positions), where value is non zero E =
diff(D)
[m k] = max(E) // maximum ’m’ at kth position modeA =
B(D(k)+1)
printf(’Mean of the given data is : %f \n \n’, meanA);
printf(’Median of the given data is : %f \n \n’, medianA);
printf(’Mode of the given data is : %f \n \n’, modeA);
printf(’First moment about the mean(M1)= %f \n \n’, 0); for
i=1:n X(i)=A(i)-meanA; end
M2 = sum(X.*X)/n;
M3 = sum(X.*X.*X)/n;
M4 = sum(X.*X.*X.*X)/n; printf(’Second moment about the
mean(M2)= %f \n \n’, M2); printf(’Third moment about the
mean(M3)= %f \n \n’, M3); printf(’Fourth moment about the
mean(M4)= %f \n \n’, M4); sd= sqrt (M2);
printf(’Standard deviation: %f \n \n’, sd); Csk=
(meanA - modeA)/sd;
printf(’Coefficient of skewness: %f \n \n’, Csk); Sk=
(M 3)∧2/(M 2)∧3;
printf(’Skewness: %f \n \n’, Sk);
Kur= M 4/(M 2)∧2;
printf(’Kurtosis: %f \n \n’, Kur); endfunction

20
4.2 Measures of Central Tendency (Grouped Data)
Algorithm to find mean, median, mode and moments of grouped data

1. Start

2. Enter the number of observations (n)

3. Input the observations using loop

4. Input frequency of each observation using loop

Σ
5. Compute sum of all frequencies fi
Σ
6. Compute fixi
Σ
f i xi

7. Find mean = Σ
fi

8. Enter how many moments to be calculated (r)

9. Compute r moments about mean in a loop

10. Evaluate S.

11. Print Mean, Moments and Standard Deviation

12. Stop
Flow Chart for Central Tendencies(Grouped Data)

21
 Scilab code for Mean and Moments of Grouped Data clc
n=input(’Enter the no. of observations:’); disp(’Enter the values of
xi’); for i=1:n
x(i)=input(’\’); end; disp(’Enter the
corresponding frequencies fi:’) sum=0; for
i=1:n
f(i)=input(’\’); sum=sum+f(i); end; r=input(’How
many moments to be calculated:’); sum1=0 for
i=1:n sum1=sum1+f(i)*x(i); end
A=sum1/sum; //Calculate the average
printf(’Average=%f \n’,A); for j=1:r
sum2=0; for i=1:n y(i)=f (i) ∗ (x(i) —
A)∧ j; sum2=sum2+y(i); end
M(j)=(sum2/sum); //Calculate the moments
printf(’Moment about mean M(%d)=%f \n’,j,M(j));
end sd=sqrt(M(2)); //Calculate the standard deviation
printf(’Standard deviation=%f \n’,sd);

Chapter 5

Plotting of 2D Graphs Using Scilab

Two Dimensional Graphs
The generic 2D multiple plot is plot2di(x,y,<options>)
index of plot2d : i = none,2, 3, 4 For the different
values of i we have: i =none : piecewise
linear/logarithmic plotting i = 2 : piecewise constant
drawing style i = 3 : vertical bars i = 4 : arrows style
//Specifier Color
//r Red
//g Green
//b Blue
//c Cyan
//m Magenta
//y Yellow
//k Black

22
//w White
//Specifier Marker Type

//+ + + Plus sign

//◦ ◦ ◦ Circle
//∗ ∗ ∗ Asterisk
//· · · Point
//× × × Cross
//’square’ or ’s’ Square
//’diamond’ or ’d’ Diamond
//∧Upward-pointing triangle
//∨ Downward-pointing triangle
//’pentagram’ or ’p’ Five-pointed star (pentagram)
5.1 A simple graph plotting raw data
clc
x = [1 -1 2 3 4 -2 ]; y =
[2 0 3 -2 5 -3 ];
plot2d(x,y) xlabel(’x’);
ylabel(’y’)

5.2 Generation of square wave

clc
x = [1 2 3 4 5 6 7 8 9 10 ]; y =
[5 0 5 0 5 0 5 0 5 0 ];
plot2d2(x,y) xlabel(’x’);
ylabel(’y’); title(’Square Wave
Function’);

23
5.3 Unit Step function I
clc
x = [-1 -2 -3 -4 -5 0 1 2 3 4 5 ]; y =
[0 0 0 0 0 1 1 1 1 1 1 ];
plot(x,y, ’ro’) xlabel(’x’);
ylabel(’y’); title(’Unit Step
Function’);

5.4 Unit Step function II

function y=unitstep2(x)
y(find (x < 0)) = 0;
y(find (x >= 0)) =
1; endfunction
clc
// define your independent values x
= [-4 : 1 : 4]’;
// call your previously defined function y
= unitstep2(x); // plot plot(x, y, ’m*’)

24
Chapter 6

Numerical Methods for Finding Roots of Algebraic

and Transcendental Equations
6.1 Bisection Method
Bisection method is used to find an approximate root in an interval by repeatedly
bisecting into sub-intervals. It is a very simple and robust method but it is also
relatively slow. Because of this it is often used to obtain a rough approximation to a
solution which is then used as a starting point for more rapidly converging methods.
The method is also called the interval halving method, binary search method or the
dichotomy method. This scheme is based on the intermediate value theorem for
continuous functions.
Let f (x) be a function which is continuous in the interval (a, b). Let f(a) be
positive and f (b) be negative. The initial approximation is x 1 = a + b
2

Then to find the next approximation to the root, one of the three conditions arises:
1. f (x1) = 0, then we have a root at x1.
2. f (x1) < 0,then since f (x1 )f (a) < 0,the root lies between x1 and a.
3. f (x1) > 0, then since f (x 1 )f (b) < 0,the root lies between x1 and b. We can further divide
this sub-interval into two halves to get a new sub-interval which contains the root. We can
continue the process till we get desired accuracy.

Bisection Method Algorithm:

1. Start

2. Read y = f (x) whose root is to be computed.

3. Input a and b where a, b are end points of interval (a, b) in which the root lies.

4. Compute f (a) and f (b).

5. If f (a) and f (b) have same signs, then display function must have different signs at a
and b, exit. Otherwise go to next step.

6. Input e and set iteration number counter to zero. Here e is the absolute error i.e. the
desired degree of accuracy.

25
 7. root = (a + b)/2

8. If f (root) ∗ f (a) > 0, then a = root else b = root

9. Print root and iteration number
10. If |a — b| > 2e, print the root and exit otherwise continue in the loop
11. Stop
Flow Chart for Bisection Method

//Scilab Code for Bisection method clc

deff('y = f (x)','y = x3 + x2 — 3 ∗ x — 3') a
=input(”enter initial interval value: ”); b
=input(”enter final interval value: ”);
//compute initial values of f (a) and f (b) fa
= f (a); fb = f (b); if sign(fa) == sign(fb)
// sanity check: f (a) and f (b) must have different signs
disp(’f must have different signs at the endpoints a and b’)
error end e=input(” answer correct upto : ”); iter=0;
printf(’Iteration \t a \t\t b \t \t root \t \t f(root)\n’) while
abs(a — b) > 2 ∗ e root
= (a + b)/2
printf(’ %i\t\t %f \t %f \t %f \t %f \n’ ,iter,a,b,root,f(root)) iff
(root) ∗ f (a) > 0 a
= root
else

26
b = root end
iter=iter+1 end printf(’\n \n The solution of given equation is %f after %i
Iterations’,root,iter —1)

6.2 Newton-Raphson Method

Newton -Raphson method named after Isaac Newton and Joseph Raphson, is a
method for finding successively better approximations to the roots of a real-valued
function. The Newton-Raphson method in one variable is implemented as follows:

Let x0 be an approximate root of the equation f (x) = 0. Next approximation x1 is given

by x1 = x 0 —'f(0x0)
f (x )

and n approximation x 1 is given by xn +1 = xn — f(xn)

th

f‘(xn)

Newton Raphson Method Algorithm:

1. Start

2. Enter the function f (x) and its first derivative f (x)

3. Take an initial guess root say x1 and error precision e.

4. Use Newtons iteration formula to get new better approximate of the root, say x2. Repeat
the process for x3, x4 .....till the actual root of the function is obtained, fulfilling the
tolerance of error.

clc
deff(’y = f (x)’,’y = x3 + x2 — 3 ∗ x — 3’)
deff(’y = df(x)’,’y = 3 ∗ x2 + 2 ∗ x — 3’)
x(1)=input(’Enter Initial Guess:’); e= input(”
answer correct upto : ”); for i = 1 : 100

27
x(i + 1) = x(i) — f (x(i))/df (x(i)); err(i)
= abs((x(i + 1) — x(i))/x(i)); if err(i) <
e

28
 Chapter 7

Numerical Integration for Solving Definite

Integrals
Numerical integration is the approximate computation of an integral using numerical
techniques. The numerical computation of an integral is sometimes called quadrature. There
are a wide range of methods available for numerical integration.

7.1 The Trapezoidal Rule

Numerical Integration by Trapezoidal Rule Algorithm:

1. Start
1. Define and Declare function y = f(x) whose integral is to be computed.
∫b
2. Input a, b and n, where a, b are lower and upper limits of integral f (x)dx
a
and n is number of trapezoids in which area is to be divided.

3. Initialize two counters sum1 and sum2 to zero.

4. Compute x(i) = a + i ∗ h and y(i) = f (x(i)) in a loop for all n + 1 points dividing n
trapezoids.
5. sum1 = y(1) + y(n) and sum2 = y(2) + . . . + y(n — 1)
6. val = 2h(sum1 + 2sum2)

7. Print value of integral.

8. Stop

Flow Chart for Trapezoidal Rule

29
 // Scilab code for Trapezoidal Rule
clc
deff(’y = f (x)',' y = x/(x 2 + 5)’); a=input(”Enter
Lower Limit: ”) b=input(”Enter Upper Limit:
”) n=input(”Enter number of sum intervals: ”)
h = (b — a)/n
sum1 = 0 sum2
= 0 for i = 0 : n
x=a+i∗hy
= f (x)
disp([xy])
if (i == 0)|(i == n) sum1 =
sum1 + y
else
sum2 = sum2 + y
end
end
30
val = (h/2) ∗ (sum1 + 2 ∗ sum2)
disp(val,”Value of integral by Trapezoidal Rule is:”)

7.2 Simpson’s 1/3rd Rule

The Simpson’s 1/3rd rule is similar to the trapezoidal rule, though it approximates
the area using a series of quadratic functions instead of straight lines. It is used if the
number of segments is even.

Numerical Integration by Simpson’s 1/3rd Rule Algorithm:

1. Start

2. Define and Declare function y = f(x) whose integral is to be computed.

3. Input a, b and n, where a, b are lower and upper limits of integral and n is number of
intervals in which area is to be divided. Note that nmust be even.

4. Put x 1 = a and initialize sum = f (a)

5. Compute x(i) = x(i — 1) + h
6. sum = sum + 4[f (x2 ) + f (x4 ) +............+ f (xn))
7. sum = sum + 2[f (x3 ) + f (x5 ) +............+ f (xn−1))
8. sum = sum + f (b)
9. val = h(sum)
3

10. Print value of integral.

11. Stop

Flow Chart for Simpson’s 1/3rd Rule

31
 // Scilab Code for Simpsons 1/3rd Rule
clc
deff(’y = f (x)',' y = x/(x 2 + 5)’); a=input(”Enter
Lower Limit: ”) b=input(”Enter Upper Limit:
”) n=input(”Enter number of sum intervals: ”)
h = (b — a)/n x(1) = a;
sum = f (a); for i = 2 : n
x(i) = x(i — 1) + h;
end
for j = 2:2:n
sum = sum + 4 ∗ f (x(j)); end for k
= 3 : 2 : n sum = sum + 2 ∗ f
(x(k)); end sum = sum + f (b); val =
sum ∗ h/3;
disp(val,”Value of integral by Simpsons 1/3rd Rule is:”)

7.3 Simpson’s 3/8th Rule

The Simpson’s 3/8th Rule rule is similar to the 1/3 rule. It is used when it is
required to take 3 segments at a time. Thus number of intervals must be a multiple of
3.

32
∫b f
a

f (xn)]

Numerical Integration by Simpson’s 3/8th Rule Algorithm:

1. Start
2. Define and Declare function y = f(x) whose integral is to be computed.

3. Input a, b and n, where a, b are lower and upper limits of integral and n is number of
intervals in which area is to be divided. Note that nmust be a multiple of 3.

4. Put x 1 = a and initialize sum = f (a)

5. Compute x(i) = x(i — 1) + h
6. sum = sum + 3[f (x 2 ) + f (x 5 ) +...........]
7. sum = sum + 3[f (x3 ) + f (x6 ) +...........])

8. sum = sum + 2[f (x 4 ) + f (x 7 ) +...........]

9. sum = sum + f (b)
10. val = 3h(sum)
8

11. Print value of integral.

12. Stop
Flow Chart for Simpson’s 3/8th Rule

33
 Scilab Code for Simpson’s 3/8th Rule
clc
deff(’y = f (x)',' y = x/(x 2 + 5)’); a=input(”Enter
Lower Limit: ”) b=input(”Enter Upper Limit:
”) n=input(”Enter number of sum intervals: ”)
h = (b — a)/n x(1) =
a; sum = f (a); for i =
2 : n x(i) = x(i — 1) +
h;
end
for j = 2 : 3 : n
sum = sum + 3 ∗ f (x(j)); end
for k = 3 : 3 : n sum = sum +
3 ∗ f (x(k)); end for l = 4 : 3 :
n sum = sum + 2 ∗ f (x(l));
end
sum = sum + f (b); val =
sum ∗ 3 ∗ h/8;
disp(val,”Value of integral by Simpson’s 3/8th Rule is:”)
34
 Chapter 8

Numerical Solution of Ordinary Differential

equations using Euler’s Method
The problem of finding a function y of x when we know its derivative and its value y0 at
a particular point x0 is called an initial value problem.

8.1 Euler’s Method

Euler’s Method provides us with an approximation for the solution of a differential
equation of the form dy = f (x, y) , y(x0) = y0. The idea behind Euler’s
dx

Method is to use the concept of local linearity to join multiple small line segments
so that they make up an approximation of the actual curve, as shown below.

The upper curve shows the actual graph of a function. The curve 0˜AA5 give
approximations using Euler’s method. Generally, the approximation gets less accurate
the further we go away from the initial value. Better accuracy is achieved when the
points in the approximation are chosen in small steps . Each next approximation is
estimated from previous by using the formula yn+1 = yn + hf (xn, yn)

Algorithm to solve an initial value problem using Euler’s method:

1. Start

2. Define and Declare function ydot representing ddyx = f (x, y)

3. Input x0, y0 , h and xn for the given initial value condition y(x0) = y0. h is step size and xn is
final value of x.

4. For x = x0 to xn solve the D.E. in steps of h using inbuilt function ode.

35
 5. Print values of x 0 to xn and y 0 to yn

6. Plot graph in (x, y)

7. Stop

Flow Chart for Euler’s Method

Scilab Code for Solving Initial Value Problem using Euler’s Method // Solution of
Initial value problem dy = 2 — 2y — e−4x, y(0) = 1, h = .1, xn = 1
dx
// If f is a Scilab function, its calling sequence must be ydot = f (x, y) // ydot
is used for first order derivative dy
dx

// ode is scilab inbuilt function to evaluate the D.E. in the given format clc
function ydot = euler(x, y)
ydot= 2 — 2 ∗ y — %e∧(—4 ∗
x) endfunction x0=input(”Enter initial
value x0: ”) y0=input(”Enter initial
value y0: ”) h=input(”Enter step size h:
”) xn=input(”Enter final value xn: ”) x
= x0 : h : xn;
y=ode(y0, x0, x, euler)
disp(x,” x value:”)
disp(y,” y value:”) plot
(x, y)

36
Chapter 9

Numerical Solution of Ordinary Differential

Equations Using RungeKutta Method

The problem of finding a function y of x when we know its derivative and its value y0 at
a particular point x0 is called an initial value problem.

9.1 Runge-Kutta method of 4th Order

Eulers method is simple but not an appropriate method for integrating an ODE as
the derivative at the starting point of each interval is extrapolated to find the next
function value. The method has first-order accuracy.
In Fourth-order Runge-Kutta method, at each step the derivative is evaluated four
times; once at the initial point; twice at trial midpoints; and once at a trial endpoint.
The final function value is calculated from these derivatives as given below.

Algorithm for solving initial value problem using RUNGE

KUTTA METHOD :

1. Start

2. Define and Declare function ydot representing ddyx = f (x, y)

3. Initialize values of x and y and input h, the step size.

4. Calculate 2 sets (i= 1,2) of values for k1, k2, k3 and k4 and subsequently value of k = 1 [k1 +
2k2 + 2k3 + k4].
37
 5. Finally yi +1 6= yi + k

6. Print values of xi and yi.

7. Stop

Flow Chart for Runge Kutta Method

// Scilab Code for Runge Kutta Method clc

function ydot = f(x, y) ydot =
x + y∧2
endfunction
x1 = 0;
y1 = 1;
h=input(”Enter step size h: ”)
x(1) = x1;
y(1) = y1;
for i = 1 : 2
k1 = h ∗ f (x(i), y(i)); k2 = h ∗ f (x(i) + 0.5 ∗ h, y(i) + 0.5 ∗
k1); k3 = h ∗ f((x(i) + 0.5 ∗ h), (y(i) + 0.5 ∗ k2)); k4 = h ∗ f
((x(i) + h), (y(i) + k3)); k = (1/6) ∗ (k1 + 2 ∗ k2 + 2 ∗ k3 +
k4); y(i + 1) = y(i) + k; printf(’\n The value of y at x=%f is %f
’,i ∗ h, y(i + 1))
x(i + 1) = x(1) + i ∗ h;
end

38