GUDLAVALLERU ENGINEERING COLLEGE

(An Autonomous Institute with Permanent Affiliation to JNTUK, Kakinada)

Seshadri Rao Knowledge Village, Gudlavalleru – 521 356.

Department of Electrical and Electronics Engineering

HANDOUT

On

NUMERICAL METHODS WITH COMPUTER

APPLICATIONS
Vision

To be a pioneer in electrical and electronics engineering education and

research, preparing students for higher levels of intellectual attainment, and
making significant contributions to profession and society.

Mission

• To impart quality education in electrical and electronics engineering

in dynamic learning environment and strive continuously for the
interest of stake holders, industry and society.
• To create an environment conducive to student-centered learning and
collaborative research.
• To provide students with knowledge, technical skills, and values to
excel as engineers and leaders in their profession.

Program Educational Objectives

1. Graduates will have technical knowledge, skills and competence to

identify, comprehend and solve problems of industry and society.
2. Graduates learn and adapt themselves to the constantly evolving
technology to pursue higher studies and undertake research.
3. Graduates will engage in lifelong learning and work successfully in teams
with professional, ethical and administrative acumen to handle critical
situations.
 HANDOUT ON NUMERICAL METHODS WITH COMPUTER APPLICATIONS

Class & Sem. : II B.Tech – I Semester Year : 2018-19

Branch : EEE Credits : 4
===============================================================

1. Brief History and Scope of the Subject

As the long history of numerical techniques indicates, numerical analysis

does not require any particular computer resources. On the other hand, the
scale and complexity of the problems that can be solved in a particular
manner are strongly influenced by the availability of high-speed computers
and efficient software implementation of numerical algorithms.
Understanding the basic numerical methods and the issues involved in
designing and analyzing them is a first step in the use of these techniques
for real-world problems. Implementing the appropriate method, or using
existing software, or a combination of the two approaches, is necessary to
achieve the final goal.

MATLAB has been chosen as the programming environment for the

presentation of numerical techniques for two reasons. First, MATLAB
provides outstanding graphing and programming capabilities, together with
the ability to solve many types of problems symbolically as well. Second,
MATLAB’s underlying matrix structure makes the software especially useful
for focusing on the aspects of various numerical techniques that can be
described conveniently in vector form. Because, vectorization is an
important approach to parallel computing.

2. Pre-Requisites

Basic Knowledge of linear algebra, and calculus

3. Course Objectives:

 To introduce the various numerical techniques using MATLAB

 To be aware of different methods to solve first order differential
equations
 To construct a curve for the given data.
4. Course Outcomes:

Upon successful completion of the course ,the students will be able to

CO1: Demonstrate various commands in MATLAB programming.
CO2: Analyse a mathematical problem and select a suitable
numerical technique to implement it in MATLAB programming.
CO3: Construct an interpolating polynomial for the given data using
MATLAB.
CO4: Find derivatives and integrals by using numerical techniques
using MATLAB.
CO5: Utilize method of least squares to fit a curve for the given data
using MATLAB.

5. Program Outcomes:

The graduates of electronics and communication engineering program

will be able to
a) apply knowledge of mathematics, science, and engineering for
solving intricate engineering problems.
b) identify, formulate and analyze multifaceted engineering problems.
c) design a system, component, or process to meet desired needs
within realistic constraints such as economic, environmental,
social, political, ethical, health and safety, manufacturability, and
sustainability.
d) design and conduct experiments based on complex engineering
problems, as well as to analyze and interpret data.
e) use the techniques, skills, and modern engineering tools necessary
for engineering practice.
f) understand the impact of engineering solutions in a global,
economic and societal context.
g) design and develop eco-friendly systems, making optimal utilization
of available natural resources.
h) understand professional ethics and responsibilities.
i) work as a member and leader in a team in multidisciplinary
environment.
j) communicate effectively.
k) manage the projects keeping in view the economical and societal
considerations.
l) recognize the need for adapting to technological changes and
engage in life-long learning.
6. Mapping of Course Outcomes with Program Outcomes:

a B c d e f g h i j k l
CO1 M H H H H
CO2 H H H H
CO3 H H H H
CO4 H H H H
CO5 H H H H

7. Prescribed Text Books

1. Laurene V. Fausett ,Applied numerical analysis using MATLAB:

2nd edition, Pearson publications, 2012, NewDelhi
2. B.S.Grewal, Higher Engineering Mathematics : 42nd edition,
Khanna Publishers, 2012 , New Delhi.
3. B.V Ramana, Higher Engineering Mathematics, Tata-Mc Graw Hill
Company Ltd.

8. Reference Text Books

1. Erwin Kreyszig, Advanced Engineering Mathematics : 8th edition,

Maitrey Printech Pvt. Ltd, 2009, Noida.
2. Robert J.Schilling, SandrabL .Harries, Applied Numerical methods for
engineers using MATLAB & C, Thomson books.
3. John.H.Mathews, kurtis D.Fink, Numerical methods using MATLAB,
4th edition-PHI.
4. Steven C.Chapra, Raymond P.Canale, Numerical methods for
Engineers, 3rd Edition, TATA McGrawhill, 2000, NewDelhi.

9. URLs and Other E-Learning Resources.

Sonet CDs & IIT CDs on some of the topics are available in the digital
library.

10. Digital Learning Materials:

https://youtu.be/R5eSBMP3XAM?list=PLEJxKK7AcSEFtiPAlDTJywUIAIF8FGYXA
https://youtu.be/6F89wOGgRf0?list=PLTwPa5Tfu7AULpYLzEGS5c2SFyuD5sOFt
https://www.youtube.com/watch?v=sZ_nCZjokQs
https://www.youtube.com/watch?v=fCKUOWiM-6s
https://www.youtube.com/watch?v=-QoZcEoGDEQ
https://www.youtube.com/watch?v=NZfd-EuBYyo
 https://www.youtube.com/watch?v=ElEqbKICvEs

11. Lecture Schedule / Lesson Plan

S.No Topics Covered Periods Tutorial

UNIT-I: Introduction to MATLAB
1 Introduction and MATLAB Environment 1
2 Basic Commands 1
3 Variables 1
1
4 Arithmetic operations 1
One dimensional array-vectors and operations
5 1
on vectors
Two dimensional array-matrices and
6 1
operations on matrices
1
7 Scripts and Functions 1
8 2D-Plotting 1
UNIT-II: Algebraic and Transcendental Equations
9 Introduction and interval calculation 1
10 Bisection Method 2 1
11 Method of False Position 2
12 Newton- Raphson Method. 2 1
UNIT-III : Interpolation
13 Introduction 1
14 Finite differences 1 1
15 Construction of difference tables & problems 2
Newton's Forward Difference formula for
16 2
interpolation 1
17 Newton's Backward difference formula for 2
 interpolation
Gauss Forward Difference formula for
18 2
interpolation
Gauss Backward difference formula for 1
19 2
interpolation
20 Lagrange's Interpolation formula 2
UNIT-IV: Numerical differentiation and integration
Introduction to Numerical differentiation and
21 1
formulae
Numerical differentiation by Newton’s Forward
22 1 1
differences
Numerical differentiation by Newton’s
23 1
Backward differences
24 Numerical Integration by Trapezoidal rule 1
25 Numerical Integration by Simpson’s 1/3 rulerd 1 1
26 Numerical Integration by Simpson’s 3/8th rule 1
UNIT-V: Numerical solution of Ordinary Differential equations
27 Introduction 1
28 Taylor’s series Method 2 1
29 Euler’s Method 1
30 Modified Euler’s Method 1
1
31 Runge-Kutta Methods of 4th order 2
UNIT-VI: Curve Fitting
32 Introduction and Method of least squares 1
33 Fitting of a straight line 1 1
34 Fitting of a Parabolic curve 1
35 Fitting of an exponential curve 1
36 Fitting of an exponential curve 1 1
37 Fitting of a power curve 1
TOTAL 48 13

12. Seminar Topics:

1. Newton – Raphson method
2. Lagranges Interpolation formula
3. Trapezoidal, Simpson 1/3 rule.
4. Derive Normal equations to fit a parabola by Method of least squares.
13. List of Lab Experiments
 S.No Title of Lab Experiments Periods

1 MATLAB Programs - function programs – Expressions - Array

2
Operations.
2 To find out the root of the Algebraic and Transcendental
2
equations using (i).Bisection, (ii). Regula – falsi method
3 To find out the root of the Algebraic and Transcendental
2
equations using Newton – Raphson method.
4 To implement Newton’s Forward and Backward Interpolation
2
formula.
5 To implement Gauss Forward and Backward interpolation
2
formula.
6
To implement Lagranges Interpolation formula. 2
7 To implement Numerical Differentiations Newton’s interpolation
2
formulae.
8 To implement Numerical Integration using Trapezoidal, Simpson
2
1/3 rule.
9 To find out the Numerical solution of ordinary differential
2
equations using Euler’s Method
10 To find out the Numerical solution of ordinary differential
2
equations using R-K Method of fourth order.
11 To implement Least Square Method to fit a Straight line and
2
Parabolic curve.
12 To implement Least Square Method to fit an exponential curve
2
and power curve.
TOTAL 24

NUMERICAL METHODS WITH COMPUTER APPLICATIONS

UNIT– I: Introduction to MATLAB

Pre-requisite : Basics of Linear algebra

Syllabus: Variables - Arrays: Vectors & Matrices - Array Operations – functions

- plots in MATLAB

Course Objectives:
 To familiarize with the environment of MATLAB and its components.
 To introduce various data types and variables in MATLAB.
 To distinguish various aspects of MATLAB scalars, arrays-vectors,
matrices.
 To get acquainted with the built-in functions on scalars, arrays-vectors,
matrices.
 To learn how to create and execute script and function files in MATLAB.
  To learn about plot functions using MATLAB.

Course Outcomes:

The students are able to

 Launch a new MATLAB session, work on desktop environment and
terminate the session
 Use various commands in MATLAB
 Create variables, scalars, arrays-vectors and matrices
 Use various built-in functions on different types of data
 Plot graphs using MATLAB commands with different style options
 Create and execute script and function files

Introduction:

 The name MATLAB stands for MATrix LABoratory

 Developed primarily by Cleve Moler, Chairman of
the Computer Science Department at the University
of New Mexica in the 1970's
 Gained its popularity through word of mouth,
because it was not officially distributed.
 The MathWorks Inc., Natick, Massachusetts, USA
was created (1984) to market and continue
development of MATLAB.

Features of MATLAB

 MATLAB is a high level interpreted programming language with

interactive environment provides vast library of mathematical functions
for numerical computation, visualization and application development.
 It is a case sensitive language.

Schematic diagram of MATLAB’s main features

MATLAB Work Environment:
Four Important Windows
 The Command Window is where we type in Matlab commands.
 The Command History Window shows the commands we have entered
in the past. We can repeat any of these commands by double-clicking on
them, or by dragging them from the Command History Window into the
Command Window. We can also scroll back to previous commands by
using the up arrow in the Command Window.
 The Workspace shows the list of variables that are currently defined,
and what type of variable each is. (i.e. a simple scalar, a vector, or a
matrix, and the size of all arrays. ) Depending on the size ( i.e. type ) of
the variable, its value may also be shown.
 The Current Directory Window shows the contents of the current
working directory.

Creating MATLAB variables

 Variable is a name made of a letter or a combination of several letters

that is assigned a value or an expression.
  MATLAB variables are created with an assignment statement.

 The syntax of variable assignment is

variable_name = a value (or an expression)
For example,
>> x =2
Where expression can involve:
 manual entry
 built-in functions
 user-defined functions
Rules for variable names
 A valid variable name starts with a letter, followed by letters, digits, or
underscores.
 MATLAB is case sensitive, so, A and a are not the same variable.
 We cannot define variables with the same names as MATLAB keywords,
such as if or end.
 For a complete list, run the iskeyword command.

Examples of valid names: Invalid names:

x6 6x
lastValue end
n_factorial n!

MATLAB System Variables

MATLAB has certain variables that are recognized by MATLAB itself and not
defined by the users.
Variable Description
ans This variable is automatically generated by MATLAB when there is
no variable assigned to store the result of an expression
inf It represents infinity, generated usually when a number is divided
by zero
eps This is a constant value representing the floating point relative
accuracy uses in its calculations
NaN It represents Not a Number. Resulting from operations like 0/0 and
inf / inf
pi This represents the constant value of π = 3.14159…
i As the basic imaginary unit sqrt(-1), i is used to enter complex
 numbers
Ex.: Z = 2+3i
j Use the character j in place of the character i, if desired, as the
imaginary unit.

MATLAB Data types

MATLAB, as a computing language, recognizes different types of data.

Scalars :
 Any number that represents magnitude (quantity or measure) is known
as scalar.
 This includes integers, complex number and floating point numbers
 For example : 5, -7, 4 + 5i, -45.18 etc.
Characters :
 The character constant is an alphanumeric symbol enclosed in a single
quote.
 The character that is not represented in single quote is numeric or sign.
 For example : ‘H’, ‘4’, ‘*’, ‘+’ etc.
Arrays :
 An array is a list of homogenous data placed in rows and/or columns
form
 The elements of an array can be numeric or character or strings, but not
mixed up.
 An array is written using square brackets enclosing its elements
separated by commas or spaces
 For Example : [4, 8, 3, -5], [cleve, moler] etc.
Strings :
 Any two or more alphanumeric symbols enclosed in a single quote is
known as string data
 For example : ‘INDIA’, ‘MATLAB’ etc.
  A string is an array of characters, i.e., ‘INDIA’ is equivalent to [‘I’, ‘N’, ‘D’,
‘I’, ‘A’]
Cell arrays :
 A cell array is a special type of arrays, where the elements can be of
different types
 The elements of a cell array are enclosed using braces
 For example : {‘south’, 544, ‘+’, ‘book’}
COMMAND HANDLING

Common System Commands

Command Description Purpose / Action
>> clc Clear command Clears the command window.
window All the variables still appear in
the workspace
>> delete filename or Delete any The file is deleted from the
delete (‘file name’) undesired files current directory
>> cd pathname or Change the Used to the change the
cd (‘pathname’) directory working directory
>> copyfile (‘source file’, Copying file Used to copy a file from source
‘destination file’) to destination
>> dir directory_name Directory Used to list the name of
dir (‘directory_name’ ) subdirectories and files under
a directory
>> date Date Displays current date
>> save Saves workspace variables in a
file
>> load Loads workspace variables
from file
>> type Displays contents of a file
Workspace Commands
Command Description Purpose / Action
>> clear Clears all variables It deletes all variables from
the current directory
>> clear x y z Clears specific It deletes specific variables
variables from the current directory
>> who Used to list out variables
used in current workspace
window.
>> whos Used to list out variables
along with their size
>> help Used to get help for MATLAB
related function.
>> help sqrt
SQRT(X) is the square root of
the elements of X. Complex
results are produced if X is
not positive.
>> quit/exit Stops MATLAB
>> exist Checks for existence of a file
or variable

SCALARS
 In MATLAB a scalar is a variable with one row and one column.
 Scalars are the simple variables that we use and manipulate in simple
algebraic equations.

Creating scalars
To create a scalar you simply introduce it on the left hand side of an equal
sign.

>> x = 1;
>> y = 2;
>> z = x + y;
Scalar operations
MATLAB supports the standard scalar operations like addition, subtraction,
multiplication and division.
>> u = 5;
>> v = 3;

Operation Description Example

Performs addition on two >> w=u+v
Addition
scalar values w= 8
Performs subtraction on >> w=u-v
Subtraction
two scalar values w= 2
Performs multiplication >> w=u*v
Multiplication
on two scalar values w = 15
Performs Division on two >> w=u/v
Division
scalar values w = 1.6667
Performs Exponentiation >> w=u^v
Exponentiation
on two scalar values w =125
Some Built-in Scalar Functions:

Certain MATLAB functions are essentially used on scalars

Function Syntax Description Example

>> sin(pi/2)
Sin sin(x) trigonometric sine
ans = 1
>> cos(pi/2)
Cos cos(x) trigonometric cosine
ans = 6.1230e-017
>> tan(pi/4)
Tan tan(x) trigonometric tangent
ans = 1
trigonometric inverse sine >>asin(pi/4)
Asin asin(x)
(arcsine) ans = 0.9033
trigonometric inverse cosine >> acos(pi/4)
Acos acos(x)
(arccosine) ans = 0.6675
trigonometric inverse tangent >> atan(pi/4)
Atan atan(x)
(arctangent) ans = 0.6658
>>exp(2)
Exp exp(x) Exponential (ex)
ans = 7.3891
>> log(2)
Log log(x) natural logarithm
ans = 0.6931
>> abs(-2)
Abs abs(x) absolute value
ans = 2
>>sqrt(16)
Sqrt sqrt(x) square root
ans = 4
>> rem(12,5)
Rem rem(x,y) remainder
ans = 2
 >> round(5.45)
ans = 5
Round round(x) round towards nearest integer
>> round(5.75)
ans = 6
>> floor(5.45)
ans = 5
Floor floor(x) round towards negative infinity
>> floor(5.75)
ans = 5
>> ceil(5.45)
ans = 6
Ceil ceil(x) round towards positive infinity
>> floor(5.75)
ans = 6

VECTORS
 In MATLAB a vector is a matrix with either one row or one column.

Creating vectors

1. using the built-in functions ones, zeros, linspace, and logspace

2. assigning a mathematical expressions involving vectors
3. appending elements to a scalar
4. using colon operator

Creating vectors with ones, zeros, linspace, and logspace

Function Syntax Description Example

ones ones(1,n) creates a row vector of >> x=ones(1,5)
length n, filled with ones x= 1 1 1 1 1

ones(n,1) creates a column vector of >>x=ones(3,1)

length n, filled with ones x=1
1
1
 zeros zeros(1,n) creates a row vector of >> x=zeros(1,5)
length n, filled with zeros x= 0 0 0 0 0

creates a column vector of >>x=zeros(4,1)

zeros(n,1) length n, filled with zeros x=0
0
0
0
linspace linspace creates a vector with >> x = linspace(1,5,5)
(begin, end, linearly spaced elements x = 1 2 3 4
no. of starting from begin to end. 5
elements)
>> x = linspace(1,5,2)
x= 1 5
logspace logspace creates a vector with >> y = logspace(1,4,4)
(begin, end, logarithmically spaced y = 10 100 1000
no. of elements starting from 10000
elements) begin to end.
>> y = logspace(1,4,2)
y = 10 10000

Creating vectors with Colon (:) operator:

MATLAB colon (:) operator is often used in creation of vectors.

>> x = initial_value : increment : final_value

 If no increment is specified, MATLAB uses the default increment of 1

Examples:
>> x = 0:10:100
x = [0 10 20 30 40 50 60 70 80 90 100]

>> a = 0:pi/50:2*pi
a = [0 pi/50 2*pi/50 . . . . . .2*pi]

>> u = 3:10
u = [3 4 5 6 7 8 9 10]
  To create a column vector, append the transpose operator to the end of
the vector-creating expression

>> y = (1:5)'

y= 1
2
3
4
5

 Using colon operator to create a vector requires you to specify the

increment, whereas using the linspace command requires you to
specify the total number of elements.

Assigning vector expressions to a vector

 Once a vector is created, it may be assigned to another vector.
 The following statements create a row vector, x, and then copies the
third through seventh elements of x into y.
>> x = linspace(31,40,10);
>> y = x(3:7)
y = 33 34 35 36 37
>> y(3)
ans = 35
Addressing vector elements
 Individual elements of a vector can be addressed with a subscript.
Example :
>> x = linspace(11,15,5);
>> x(2)
ans = 12

Increasing the size of a vector (or scalar)

We can increase the size of a vector simply by assigning a value to an element
that has not been previously used.

>> x = linspace(21,25,5)
x = 21 22 23 24 25

>> x(7) = -9
 x = 21 22 23 24 25 0 -9

Vector Operations:

>> x = [1 2 3] >> y = [3 4 6]

Operation Description Example

+ Addition of two vectors >> Z=X+Y
Z= 4 6 9
- Subtraction of two vectors >>Z=X-Y
Z = -2 -2 -3
* Multiplication of two vectors >> Z=X * Y'
Z = 29
.* Element-by-element >> x.*y
multiplication ans = 3 8 18
./ Element-by-element left division >> x./y
ans = 0.3333 0.5000
0.5000
.\ Element-by-element right division >> x.\y
ans = 3 2 2
.^ Element-by-element >> x.^y
exponentiation ans = 1 16

Some Built-in vector functions:

The following MATLAB functions operate on vectors and return a scalar value.

Let z be the following row vector.

>> z = [34, 5, 11, 90]

Function Description Result

Largest component >>max(z)
Max
ans = 90
Smallest component >>min(z)
Min
ans = 5
length of a vector >>length(z)
Length
ans = 4
 sort in ascending
>>sort(z)
Sort
order ans = [5,11,34,90]

sum of elements >>sum(z)

Sum
ans = 140
product of elements >>prod(z)
Prod
ans = 168300
median value >>median(z)
median
ans = 22.5000
mean value >>mean(z)
Mean
ans = 35
Standard deviation >>std(z)
Std
ans = 38.7384

MATRICES

Creating Matrices

Function Description Example Result

eye Identity matrix » eye(3) 100
010
001

» eye(3,4) 1 0 0 0
0 1 0 0
0 0 1 0
zeros matrix of zeros » zeros(2) 0 0
0 0

» zeros(2,3) 0 00
0 00
Ones matrix of ones » ones(2) 1 1
1 1

» ones(2,3) 1 1 1
1 1 1
diag extract diagonal of a x= 1 3 4
matrix or create 4 5 7
diagonal matrices 5 9 0 1
>> diag(x) 5
0
triu upper triangular part of >> triu(x) 1 3 4
a matrix 0 5 7
 0 0 0
tril lower triangular part of >> tril(x) 1 0 0
a matrix 4 5 0
5 9 0
rand randomly generated >> rand(3) 0.9501 0.4860
matrix 0.4565
0.2311 0.8913
0.0185
0.6068 0.7621
0.8214
 The above (vector) commands can also be applied to a matrix.
 In this case, they act in a column-by-column fashion to produce a row vector
containing the results of their application to each column.
Matrix operations
Example:
A = [ 1 2 3 ; 4 5 6; 7 8 9];
B = [ 7 5 6 ; 2 0 8; 5 7 1];
Operation Oper Description Example
ator
Addition + cij = aij + bij >> C=A+B
C= 8 7 9
6 5 14
12 15 10
Subtraction - cij = aij - bij >> C=A-B
C = -6 -3 -3
2 5 -2
2 1 8
Multiplication * n >> C=A*B
C = 26 26 25
cij = ∑ aik * bkj 68 62 70
k=1 110 98 115
Right division / C = A/B >> C=A/B
=> A*inv(B) C = -0.5254 0.6864 0.6610
-0.4237 0.9407 1.0169
-0.3220 1.1949 1.3729
Left division \ C = A\B A = [4,5,3; 2,1,4; 3,2,6]
=> inv(A)*B B = [1,2,3; 4,5,6; 7,8,9]
>>C = A\B
C = -4.2000 -2.4000 -0.6000
2.0000 1.0000 -0.0000
2.6000 2.2000 1.8000
Exponentiation ^ n >> C=A^2
C = 30 36 42
cij = ∑ aik * akj
 k=1
66 81 96
102 126 150
Element-by- .* cij = aij * bij >>C = A .* B
element
multiplication C = [10, 40, 90 ; 160, 250, 360; 490,
640, 810]
Element-by- ./ cij = aij / bij >>C = A ./ B
element left
division C = [ 10, 10, 10; 10, 10, 10; 10, 10, 10
]
Element-by- .^ cij = aij ^ x >>C = A .^ 2
element
exponentiation C = [ 1, 4, 9; 16, 25, 36; 49, 64, 81 ]

Some Built-in Matrix functions:

A = [1 3; 2 4]

Function Description Example Result

Size size of a matrix >>Size(A) ans = 2 2
det Determinant of a square >>det(A) ans = -2
matrix
inv inverse of a matrix >>inv(A) ans = -2.0000 1.5000
1.0000 -0.5000
rank rank of a matrix >>rank(A) ans = 2
eig Eigen values and eigen >>eig(A) ans = -0.3723
vectors 5.3723
Produces a diagonal matrix >>[x d]=eig(A) x = -0.9094 -0.5658
‘d’ with the eigen values on 0.4160 -0.8246
the main diagonal, and a full d = -0.3723 0
matrix ‘x’ whose columns 0 5.3723
are the corresponding
eigenvectors.
max Produces a row vector with >>max(A) ans = 2 4
maximum element in each
column
>>max(max(A)) ans = 4
Maximum of all elements of
matrix
PLOTTING

 MATLAB has an excellent set of graphic tools.

  Plotting a given data set or the results of computation is possible with
very few commands.

The MATLAB command to plot a graph is plot(x,y).

Example:
The vectors
x = (1; 2; 3; 4; 5; 6)
y = (3;-1; 2; 4; 5; 1)
>> x = [1 2 3 4 5 6];
>> y = [3 -1 2 4 5 1];
>> plot(x,y)

Adding titles, axis labels, and annotations:

Multiple data sets in one plot

Multiple (x; y) pairs arguments create multiple graphs with a single call to plot.
For example,

>> x = 0:pi/100:2*pi;
>> y1 = 2*cos(x);
>> y2 = cos(x);
>> y3 = 0.5*cos(x);
>> plot(x,y1,'--',x,y2,'-',x,y3,':')
>> xlabel('0 \leq x \leq 2\pi')
>> ylabel('Cosine functions')
>>
legend('2*cos(x)','cos(x)','0.5*cos(x)')
>> title('Typical example of multiple
plots')
>> axis([0 2*pi -3 3])
Specifying line styles and colors

Using plot command, we can specify line styles, colors, and markers
(e.g., circles, plus signs, . . )

plot(x, y, 'style_color_marker')

where style_color_marker is a triplet of values

Example:
Time and Temp
55

% A simple plot of just one point!

% Create coordinate variables and plot a red ‘*’
x = [11];
50

y = [48];
Temperature

plot(x, y, 'r*') 45

% Change the axes and label them

axis([9 20 35 55]) 40

xlabel('Time')
ylabel('Temperature')
% Put a title on the plot
35
9 10 11 12 13 14 15 16 17 18 19 20
Time

title('Time and Temp')

Example: Time and Temp

55

x = [11,15,18];
y = [48,50,51]; 50
plot(x, y, ':gx')
% Change the axes and label them
Temperature

axis([9 20 35 55]) 45

40

35
9 10 11 12 13 14 15 16 17 18 19 20
Time
xlabel('Time')
ylabel('Temperature')
% Put a title on the plot
title('Time and Temp')
Introduction to programming in MATLAB

The commands entered in the Command Window cannot be saved and

executed again for several times. To execute the commands repeatedly with
MATLAB is:
1. to create a file with a list of commands,
2. save the file, and
3. run the file.
 These files are called script files or scripts
 They must have the file extension “.m”
 Corrections or changes can be made to the commands in the file
 There are two types of m-files: script files and function files.
 M-files can be scripts that simply execute a series of MATLAB
statements, or they can be functions that can accept arguments and can
produce one or more outputs.

Script files: A script file is an external file that contains a sequence of Matlab
statements. By typing the filename, subsequent Matlab input is obtained from
the file. Script files have a filename extension of .m and are often called M-files.
Example:

Consider the system of equations:

Find the solution x to the system of equations.

Solution:
 Use the MATLAB editor to create a file: File -> New -> M-file.
 Enter the following statements in the file:
A = [1 2 3; 3 3 4; 2 3 3];
b = [1; 1; 2];
x = A\b
 Save the file, for example, example1.m.
  Run the file, at the command prompt , by typing:
>> example1
x=
-0.5000
1.5000
-0.5000
Example:

Function files: Function file is a script file (M-file) that adds a function
definition to Matlab’s list of functions.
Syntax:
function [output variables] = function_name ( input variables) ;

Example 1:
function y = myfunc (x)
y = 2*x^2 - 3*x + 1;
end
Save this file as: myfunc.m in your working directory.
This file can now be used in the command window just like any predefined
Matlab function. In the command window enter:
x = -2:.1:2; % Produces a vector of x values
y = myfunc(x); % Produces a vector of y values
plot (x,y)
Example 2:
Functions can have multiple inputs, which are separated by commas.
For example:

function y = myfunc2d (x,p)

y = 2*x.^p - 3*x + 1;
end

 Functions can have multiple outputs, which are collected into a vector

function [x2 x3 x4] = mypowers (x)

x2 = x .^2;
x3 = x .^3;
x4 = x .^4;
end

We can use the results of the program to make graphs:

x = -1:.1:1
[x2 x3 x4] = mypowers (x);
plot (x,x,'black',x,x2,'blue',x,x3 ,'green',x,x4,'red')

Difference between Script files and Function files:

Script files Function files

Do not accept input arguments or Can accept input arguments and
return output arguments return output arguments
Store variables in a workspace that is Store variables in a workspace
shared with other scripts internal to the function
Are useful for automating a series of Are useful for extending the
commands MATLAB language for your
application
 UNIT-II : Algebraic and Transcendental Equations

Pre-requisite : Commands of MATLAB

Syllabus : Solution of Algebraic and Transcendental Equations-Introduction -

Bisection Method - Method of False Position - Newton-Raphson Method.

Course Objectives:

 To recongnize the algebraic and Transcendental Equations.

 To Understand the Bisection method, method of False Position and
Newton Raphson Method.
 To Implement Bisection, False Position and Newton Raphson Methods in
MATLAB

Course Outcomes:

Students will be able to

 Solve an Algebraic and Transcendental equation using Numerical

Methods
 Find the roots of non-linear equations using MATLAB programs.

Solutions of Algebraic and Transcendental equations

Introduction : A problem of great importance in science and engineering is

that of determining the roots/ zeros of an equation of the form f(x) = 0

 Polynomial function: A function f(x) is said to be a polynomial function

if f(x) is a polynomial in x.
i.e. f(x)  a0 x  a1 xn1  ............  an1 x  an where a0  0 ,
n

the coefficients a0 , a1...........an are real constants and

n is a non-negative integer.

 Algebraic function: A function which is a sum (or) difference (or) product of

two polynomials is called an algebraic function. Otherwise, the function is
called a transcendental (or) non-algebraic function.
Eg: f  x   c1e x  c2e x
x3
f  x   e5 x  3
2
 Algebraic Equation: If f(x) is an algebraic function, then the equation f(x) =0
is called an algebraic equation.
 Transcendental Equation: An equation which contains polynomials,
 exponential functions, logarithmic functions and Trigonometric functions
etc. is called a Transcendental equation.
Ex:- xe2x– 1 = 0, cos x – x ex= 0, tan x = x are transcendental equations.

 Root of an equation: A number α is called a root of an equation f(x) = 0

if f(α) =0. We also say that α is a zero of the function.

Note: (1) The roots of an equation are the abscissas of the points where the
graph y = f(x) cuts the x-axis.
(2) A polynomial equation of degree n will have exactly n roots, real or
complex, simple or multiple. A transcendental equation may have
one root or infinite number of roots depending on the form of f(x).

Methods for solving the equation

Direct method:
We know the solution of the polynomial equations such as linear equation
ax  b  0 and quadratic equation ax 2  bx  c  0 , will be obtained using direct
methods or analytical methods. Analytical methods for the solution of cubic
and quadratic equations are also well known to us.

There are no direct methods for solving higher degree algebraic equations or
equations involving transcendental functions. Such equations are solved by
numerical methods. In these methods we find an interval in which the root lies.

We use the following theorem of calculus to determine an initial approximation.

It is also called the Intermediate value theorem.

Intermediate value theorem : If f(x)is continuous on some interval [a, b]and

f(a)f(b)< 0, then the equation f(x) = 0 has at least one real root in the interval
(a, b).

In this unit we will study some important methods of solving algebraic and
transcendental equations.

Bisection method:
Bisection method is a simple iteration method to solve an equation. This
method is also known as "Bolzano method of successive bisection". Sometimes
it is referred to as "Half-interval method". Suppose we know an equation of the
form f  x   0 has exactly one real root between two real numbers x0 , x1 . The
number is chosen such that f  x0  and f  x1  will have opposite sign. Let us
bisect the interval  x0 , x1  into two half intervals and find the midpoint
x  x1
x2  0 . If f  x2   0 then x2 is a root.
2
If f  x1  and f  x2  have same sign then the root lies between x0 and x2. The
interval is taken as  x0 , x2  Otherwise the root lies in the interval  x2 , x1  .
Repeating the process of bisection , we obtain successive subintervals which
are smaller. At each iteration, we get the mid-point as a better approximation of
the root. This process is terminated when interval is smaller than the desired
accuracy.

Problem:-Find a root of the equation x 3  5 x  1  0 using the bisection method

in 5 – stages
Sol: Let f  x   x3  5 x  1
we note that f(0) 0 and f(1)
 Root lies between 0 and 1

Consider x0  0 and x1  1
By bisection method the next approximation is
x x 1
x2  0 1   0  1  0.5
2 2
 f  x2   f  0 : 5   1.375  0 and f  0   0
We have the root lies between 0 and 0.5
0  0.5
Now x3   0.25
2
We find f  x3   0.234375  0 and f  0   0
Since f  0   0 , we conclude that root lies between x0 and x3
The third approximation of the root is
x x 1
x4  0 3   0  0.25 
4 2
 0.125
We have f  x4   0.37495  0
Since f  x4   0 and f  x3   0 , the root lies between
x4  0.125 and x3  0.25
Considering the 4th approximation of the roots
x  x4 1
x5  3   0.125  0.25   0.1875
2 2
f  x5   0.06910  0 ,
 since f  x5   0 and f  x3   0 the root must lie between x5  0.18758 and
x3  0.25
Here the fifth approximation of the root is
1
x6   x5  x3 
2
1
  0.1875  0.25 
2
 0.21875
We are asked to do up to 5 stages. We stop here and 0.21875 is taken as an
approximate value of the root and it lies between 0 and 1

MATLAB Program for Bisection method

function c = bisection(f,a,b)

if f(a)*f(b)>0
disp('Interval has no root')
else
c = (a + b)/2;
while abs(f(c)) > 1e-7
if f(a)*f(c)> 0
a = c;
else
b = c;
end
c = (a + b)/2;
end
end

Output

>> f=@(x)x^3-5*x+1;
>> bisect(f,0,1)

ans =
0.2016
False Position Method ( Regula – Falsi Method)

In the false position method we will find the root of the equation f  x   0 .
Consider two initial approximate values x0 and x1 near the required root so
that f  x0  and f  x1  have different signs. This implies that a root lies between
x0 and x1 . The curve f  x  crosses x- axis only once at the point x2 lying
between the points x0 and x1 , Consider the point A   x0 , f  x0   and B   x1 , f  x1  
on the graph and suppose they are connected by a straight line, Suppose this
line cuts x-axis at x2 , We calculate the values of f  x2  at the point. If
f  x0  and f  x2  are of opposite sign, then the root lies between x0 and x2 and
value x1 is replaced by x2
Otherwise the root lies between x2 and x1 and the value of x0 is replaced by x2 .
Another line is drawn by connecting the newly obtained pair of values. Again
the point here the line cuts the x-axis is a closer approximation to the root.
This process is repeated as many times as required to obtain the desired
accuracy. It can be observed that the points x2 , x3 , x4 .....obtained converge to
the expected root of the equation y  f  x 

To obtain the equation to find the next approximation to the root

Let A   x0 , f  x0   and B   x1 , f  x1   be the points on the curve y  f  x  Then the

y  f  x0  f  x1   f  x0 
equation to the chord AB is   1
x  x0 x  x0
At the point C where the line AB crosses the x – axis, we have f(x) =0 i.e. y =0
x1  x0
From (1), we get x  x0  f  x0    2 
f  x1   f  x0 
x is given by (2) serves as an approximated value of the root, when the interval
in which it lies is small. If the new values of x is taken as x2 then (2) becomes
 x1  x0  f x
x2  x0   0
f  x1   f  x0 
x f  x1   x1 f  x0 
 0   3 (3)
f  x1   f  x0 

Now we decide whether the root lies between 0

x and x
2  or  x
2 and x
1

We name that interval as  x1 , x2  . The line joining  x1 , y1  x2 , y2  meets x – axis at

x1 f  x2   x2 f  x1 
x3 is given by x3 
f  x2   f  x1 
This will in general, be nearest to the exact root we continue this procedure till
the root is found to the desired accuracy.
The iteration process based on (3) is known as the method of false
position. The successive intervals where the root lies, in the above procedure
are named as  x0 , x1  ,  x1 , x2  ,  x2 , x3  etc
Where xi  xi  1 and f  x0  , f  xi  1 are of opposite signs
xi 1 f  xi   xi f  xi 1 
Also xi 1 
f  xi   f  xi 1 
Problem:-
Find out the roots of the equation x3  x  4  0 using false position method
Sol: Let f  x   x3  x  4  0
f  0   4, f 1  4, f  2   2
Since f 1 and f  2  have opposite signs the root lies between 1 and 2
x0 f  x1   x1 f  x0 
By false position method x2 
f  x1   f  x0 

x2 
1 2   2  4 
2   4 
2  8 10
   1.666
6 6
f 1.666   1.666   1.666  4
3

 1.042
Now, the root lies between 1.666 and 2
1.666  2  2   1.042 
x3   1.780
2   1.042 
f 1.780   1.780   1.780  4
3

 0.1402
Now, the root lies between 1.780 and 2
1.780  2  2   0.1402 
x4   1.794
2   0.1402 
f 1.794   1.794   1.794  4
3

 0.0201
Now, the root lies between 1.794 and 2
1.794  2  2   0.0201
x5   1.796
2   0.0201
f 1.796   1.796   1.796  4  0.0027
3

Now, the root lies between 1.796 and 2

 1.796  2  2   0.0027 
x6   1.796
2   0.0027 
The root is 1.796

MATLAB Program for False Position method

function c=falsePos(f, a, b)
%False Position method for nonlinear equation
if f(a)*f(b) > 0
disp('Interval has no root')
else
c = (a*f(b)-b*f(a))/(f(b)-f(a));
while abs(f(c)) > 1e-7
if f(a)*f(c) > 0
a=c;
else
b=c;
end
c = (a*f(b)-b*f(a))/(f(b)-f(a));
end

Output:

>> f=@(x)x^3-x-4;
>> falsePos(f,1,2)

ans =

1.7963

Newton- Raphson Method:

The Newton-Raphson method is a powerful and elegant method to find the root
of an equation. This method is generally used to improve the results obtained
by the previous methods.
Let x0 be an approximate root of f  x   0 and let x1  x0  h be the correct root
which implies that f  x1   0 .
By Taylor’s theorem neglecting second and higher order terms
f  x1   f  x0  h   0
 f  x0   hf 1  x0   0
f  x0 
h
f 1  x0 
Substituting this in x1 we get
x1  x0  h
f  x0 
 x0 
f 1  x0 
 x1 is a better approximation than x0
Successive approximations are given by

f  xn 
xn 1  xn 
f 1  xn 

Problem:-
Find by Newton’s method, the real root of the equation xe x  2  0 Correct to
three decimal places.
Sol: Let f  x   xe  2  1
x

Then f  0   2 and f 1  e  2  0.7183

So root of f  x  lies between 0 and 1
It is near to 1. so we take x0  1 and f  x   xe  e and f 1  e  e  5.4366
1 x x 1

 By Newton’s Rule
f  x0 
First approximation x1  x0 
f 1  x0 
0.7183
 1  0.8679
5.4366
 f  x1   0.0672 f 1  x1   4.4491
f  x1 
The second approximation x2  x1 
f 1  x1 
 0.0672
 0.8679 
4.4491
 0.8528
 Required root is 0.853 correct to 3 decimal places.

MATLAB Program for Newton-Raphson method

function c=newtonR(f,df,a)

if(df(a)==0) disp('differentiation of a is zero');

else
c=a-f(a)/df(a);
while abs(f(c)) > 1e-7
a=c;
c=a-f(a)/df(a);
end
end

Output:

>> f=@(x)x^3-x-4;
>> df=@(x)3*x^2-1;
>> newtonR(f,df,1)

ans =

1.7963
 UNIT-III

INTERPOLATION

Objectives:
 Develop an understanding of the use of numerical methods in modern
scientific computing.
 To gain the knowledge of Interpolation

Syllabus:
Pre-requisite : Commands of MATLab
Interpolation – Introduction - Finite differences - Forward Differences -
Backward differences - Central differences - Newton Interpolation formulae -
Gauss Interpolation formulae - Lagrange’s interpolation.

Learning Outcomes:
Student should be able to
 Know about the Interpolation, and Finite Differences.
 Utilize the Newton’s formulae for interpolation.
 Utilize the Gauss formulae for interpolation.
 Operate Lagrange’s Interpolation formula.

Introduction:
Consider the statement y =f(x) , x0  x  xn we understand that we can find
the value of y, corresponding to every value of x in the range x 0  x  xn. If the
function f(x) is single valued and continuous and is known explicitly then the
values of f(x) for certain values of x like x0, x1,…….xn can be calculated.
Now the problem is, if we are given the set of tabular values

Satisfying the relation y = f(x) and the explicit definition of f(x) is not
known, is it possible to find a simple function say  (x) such that f(x) and  (x)
agree at the set of tabulated points. This process of finding  (x) is called
interpolation. If  (x) is a polynomial then the process is called polynomial
interpolation and  (x) is called interpolating polynomial. In our study we are
concerned with polynomial interpolation.
Interpolation is the process of deriving a simple function from a set of
discrete data points so that the function passes through all the given data
points (i.e. reproduces the data points exactly) and can be used to estimate
data points in-between the given ones.

Finite Differences:
 Here we introduce forward, backward and central differences of a
function y = f(x). These differences play a fundamental role in the study of
differential calculus, which is an essential part of numerical applied
mathematics.

1. Forward Differences:
Consider a function y = f(x) of an independent variable x. let
y0,y1,y2…….yn be the values of y corresponding to the values x0,x1,x2…….xn of x
respectively. Then the differences y1-y0, y2 – y1.............. are called the first
forward differences of y, and we denote them by y0 , y1 ,........
that is y0  y1  y0 , y1  y2  y1 , y2  y3  y2 ,...........
In general yr  yr 1  yr where r  0,1, 2,3,......n  1
Here the symbol  is called the forward difference operator.
The second forward differences and are denoted by 2 y0 , 2 y1 , 2 y2 ,...........
that is

In general 2 yr  yr 1  yr where r  0,1, 2,3,......n  2

similarly, the kth forward differences are defined by the formula.
k yr  k 1 yr 1  k 1 yr , where r  0,1, 2,.....n  k and k  1, 2,......n 1
The symbol k is referred as the kth forward difference operator.
If f(x) is a function of x and h be the increment in x then forward
difference of f(x) is defined as f ( x)  f ( x  h)  f ( x)
Forward Difference Table:
The forward differences are usually arranged in tabular columns
as shown in the following table called a forward difference table

Values Values First order Second order Third order

of x of y differences differences differences

2. Backward Differences:
Consider a function y = f(x) of an independent variable x. let
y0,y1,y2…….yn be the values of y corresponding to the values x0,x1,x2…….xn of x
respectively. Then the differences y1-y0, y2 – y1.............. are called the first
backward differences of y, and we denote them by y1 , y2 ,........
that is y1  y1  y0 , y2  y2  y1 , y3  y3  y2 ,...........
In general yr  yr  yr 1 where r  1, 2,3,......n
Here the symbol  is called the backward difference operator.
The second backward differences and are denoted by 2 y2 , 2 y3 , 2 y4 ,...........
That is

In general 2 yr  yr  yr 1 where r  2,3,......n

similarly, the kth backward differences are defined by the formula as
k yr  k 1 yr  k 1 yr 1 , where r  k , k  1,.....n and k  1, 2,......n 1
The symbol  k is referred as the kth backward difference operator.
If f(x) is a function of x and h be the increment in x then bckward
difference of f(x) is defined as f ( x)  f ( x)  f ( x  h)

Backward Difference Table:-

X Y

3.Central Differences:
Consider a function y = f(x) of an independent variable x. let
y0,y1,y2…….yn be the values of y corresponding to the values x0,x1,x2…….xn of x
respectively. We define the first central differences
as follows

The symbol is called the central differences operator. This operator is a linear
operator
Comparing expressions (1) above with expressions earlier used on forward and
backward differences we get
 In general
The first central differences of the first central differences are called the
second central differences and are denoted by
Thus

Higher order central differences are similarly defined. In general the nth
central differences are given by
for odd
for even
while employing for formula (4) for , we use the notation
If y is a constant function, that is if a constant, then
If f(x) is a function of x and h be the increment in x then forward
 h  h
difference of f(x) is defined as  f ( x)  f  x    f  x  
 2  2

Central Difference Table:

Note: The forward difference operator  , backward difference operator  and

central difference operator  are linear operators.

E-operator: The shift operator E is defined by the equation Eyr  yr 1 . This

shows that the effect of E is to shift the functional value y r to the next higher
value yr 1 . This is also called forward shift operator.
A second operation with E gives E 2 yr  E( Eyr )  Eyr 1  yr 2
Generalizing E k yr  yk r
Inverse operator E-1: Inverse operator E 1 is defined as E 1 yr  yr 1. This is also
called backward operator.
In general E  k yr  yr k
Relationship Between  and E
We have

Some more relations

Note: We can easily establish the following relations

Theorem: If f(x) is a polynomial of degree n and the values of x are equally

spaced then  n f ( x ) is constant.
Note: As  n f ( x ) is a constant, it follows that  n 1 f ( x)  0,  n  2 f ( x)  0,........ The
converse of above result is also true that is, if  n f ( x ) is tabulated at equal
spaced intervals and is a constant, then the function f(x) is a polynomial of
degree n.
Q.
Q.Find the missing term in the following data

X 0 1 2 3 4
Y 1 3 9 - 81
Why this value is not equal to . Explain
Sol. Consider

Substitute given values we get

From the given data we can conclude that the given function is . To
find , we have to assume that y is a polynomial function, which is not so.
Thus we are not getting

Q.Evaluate

Sol. Let h be the interval of differencing

Proceeding on, we get

Newton’s Forward Interpolation Formula:

Let y=f(x) be a polynomial of degree n and taken in the following form

This polynomial passes through all the points (for i = 0 to n. Therefore, we can
obtain the yi ' s by substituting the corresponding xi ' s as

Let ‘h’ be the length of interval such that represent

This implies
From (1) and (2), we get
Solving the above equations for , we get

Similarly, we can see that

If we use the relationship

Then

Equation (3) becomes

Q. Find the melting point of the alloy containing 54% of lead, using appropriate
interpolation formula
Percentage of
50 60 70 80
lead(p)
Temperature 205 225 248 274

Sol. The difference table is

x Y
50 205
20
60 225 3
23 0
70 248 3
26
80 274

Let temperature =

By Newton’s forward interpolation formula

Melting point = 212.6

Q. Consider the following table of values

Sol.

MatLab Code for Newton’s Forward Interpolation Formula:

function fp = newtonint( x,y,p )

%To find a function value y at certain value of x when we have %x values and
their corresponding y values
% By using Newton’s Forward Interpolation Method
% x is data for x
% y is data for y
% p is a point where we have to calculate y value
n=length(x);
for i= 1:n
 diff(i,1)=y(i);
end
for j=2:n
for i=1:n-j+1
diff(i,j)=diff(i+1,j-1)-diff(i,j-1);
end
end
fp=y(1);
h=x(2)-x(1);
u=(p-x(1))/h;
for i=1:n-1
factor=1;
for j=0:i-1
factor=factor*(u-j);
end
fp=fp+factor*diff(1,i+1)/factorial(i);
end
end

Newton’s Backward Interpolation Formula:

If we consider

and impose the condition that y and should agree at the tabulated points

We obtain

Where
This uses tabular values of the left of yn. Thus this formula is useful
formula is for interpolation near the end of the tabular values.
Q. The sales for the last five years is given in the table below. Estimate the
sales for the year
1979

Sol.
Q. Consider the following table of values

Use Newton’s Backward Difference Formula to estimate the value of f(1.45).

Sol.
MatLab Code for Newton’s Backward Interpolation Formula:

function yval=nbdi(xd,yd,xval)
n=length(xd);
bdt(:,1)=xd';
bdt(:,2)=yd';
for j=3:n+1
for i=n:-1:j-1
bdt(i,j)=bdt(i,j-1)-bdt(i-1,j-1);
end
end
bdt
h=xd(2)-xd(1)
p=(xval-xd(n))/h
c=ones(1,n);
for r=1:n-1
c(r+1)=prod(p:1:p+r-1)/factorial(r);
end
terms=bdt(n,2:n+1).*c;
yval=sum(terms);
fprintf('\n The approximated y value at given x=%f is: %f',xval,yval)
end
Gauss Forward Interpolation Formula:
Gauss forward interpolation formula is given by

p( p  1) 2 ( p  1) p( p  1) 3 ( p  1) p( p  1)( p  2) 4
y p  y0  py0   y1   y1   y2
2! 3! 4!
( p  2)( p  1) p( p  1)( p  2) 5
  y2  .....
4!
Where

 The value p is measured forwardly from the origin and 0<p<1

 The above formula involves odd differences below the central horizontal
line and even differences on the line. This is explained in the following
figure.

Q. Find f(30) from the following table values using Gauss forward difference
formula:

x: 21 25 29 33 37

F(x): 18.4708 17.8144 17.1070 16.3432 15.5154

Sol:

The difference table is

x f f  2f  3f  4f
21 18.4708
-0.6564
25 17.8144 -0.0510
-0.7074 -0.0054
29 17.1070 0.0564 -0.0022
-0.7638 -0.0076
33 16.3432 -0.0640
-0.8278
37 15.5154

Let given x=30 and h=4 then

Gauss forward interpolation formula is given by

p ( p  1) 2 ( p  1) p ( p  1) 3 ( p  1) p( p  1)( p  2) 4
y p  y0  py0   y1   y1   y2  .....
2! 3! 4!
Therefore f (30) = 16.9217

MatLab Code for Gauss Forward Interpolation Formula:

function [ yval ] = gauss_p( xd,yd,xp )

n=length(xd);
if(length(yd)==n)
tbl=yd';
for j=2:n
for i=1:n-j+1
tbl(i,j)=tbl(i+1,j-1)-tbl(i,j-1);
end
end
tbl
h=xd(2)-xd(1);
if rem(n,2)==0
k=n/2+1;
else
k=n/2+0.5;
end

p=(xp-xd(k))/h;

pt=cumprod([1,p-(0:n-3)]);

dt=[tbl(k,1),tbl(k,2),tbl(k-1,3:n-1)+tbl(k-1,4:n)]./factorial(0:n-2);

yval=sum(pt.*dt);

end

Gauss Backward Interpolation Formula:

Gauss backward interpolation formula is given by
 ( p  1) p 2 ( p  1) p( p  1) 3 ( p  2)( p  1) p ( p  1) 4
y p  y0  py0   y1   y1   y2
2! 3! 4!
( p  2)( p  1) p( p  1)( p  2) 5
  y2  .....
4!
Where

 The value p is measured forwardly from the origin and -1<p<0.

 The above formula involves odd differences above the central horizontal
line and even differences on the line.

Q. From the following data find y when x=38 by using Gauss backward
interpolation formula
x 30 35 40 45 50
y 15.9 14.9 14.1 13.3 12.5

Sol: The difference table is

x y y 2 y 3 y 4 y
30 15.9
-1
35 14.9 0.2
-0.8 -0.2
40 14.1 0 0.2
-0.8 0
45 13.3 0
-0.8
50 12.5

Let given x=38 and h=5 then

Gauss forward interpolation formula is given by

( p  1) p 2 ( p  1) p( p  1) 3 ( p  2)( p  1) p ( p  1) 4
y p  y0  py0   y1   y1   y2
2! 3! 4!
( p  2)( p  1) p( p  1)( p  2) 5
  y2  .....
4!
Therefore y(38)=14.4245

Lagrange’s Interpolation Formula:

Let be the values of x which are not necessarily
equally spaced. Let be the corresponding values of y  f ( x) let the

polynomial of degree n for the function y  f ( x) passing through the

points be in
the following form

Where an are constants

Since the polynomial passes through , .
The constants can be determined by substituting one of the values of
in the above equation
Putting in (1) we get,

Putting in (1) we get,

Similarly substituting in (1), we get

Continuing in this manner and putting in (1) we get

Substituting the values of , we get

Q . Using Lagrange's formula calculate from the following table

x 0 1 2 4 5 6
 1 14 15 5 6 19

Sol. Given

From lagrange’s interpolation formula

Here then

MatLab Code for Lagrange’s Interpolation Formula:

function [ yval ] = lagrange(xd,yd,x)

%To find a function value y at certain value of x when we have x values and
%their corresponding y values by using Lagranges Interpolation Method
% xd is data for x
% yd is data for y
% x is a point where we have to calculate y value
n=length(xd);
if(length(yd)==n)
p=zeros(1,n);
for i=1:n
temp=xd;
temp(i)=[];
p(i)=prod((x-temp)./(xd(i)-temp));
end
yval=sum(p.*yd);
fprintf('\n The value of y at x=%f is %f',x,yval)
else
error('xd and yd must be of same size');
end

Assignment-Cum-Tutorial Questions

Section-A
Objective Questions:

1. A linear version of the Lagrange’s interpolation formula for f(x) is [ ]

a) b)

c) d)

2. The following is used for unequal interval of x values [ ]

a) Lagrange’s formula
b) Newton’s forward interpolation formula
c) Newton’s backward interpolation formula
d) Gauss forward interpolation formula

3. The order difference a polynomial of nth degree is [ ]

a) polynomial of nth degree b) zero
c) polynomial on first degree d) constant

X 1 2 3 4
4. f(x) 1 4 27 64

If x=2.5 then p= [ ]
a) 1.5 b) 1 c) 2.5 d) 2
5. X 0.1 0.2 0.3 0.4
f(x) 1.005 1.02 1.045 1.081
When p=0.6, x= [ ]
a) 0.16 b) 0.26 c) 0.1 d) 3.0 2.

6. Relation between Backward and Shifting operator is ________________.

7. When do we apply Lagrange’s interpolation?

8. If then [ ]
a) 1 b) 2 c) 0 d) 3

9. [ ]

a) b) c) d)

10. [ ]
a) b) c) d)

11.
X 0 1 2
F(x) 7 10 13

By Newton’s forward formula f(2.5)= [ ]

a) 15.25 b) 16.75 c) 16.25 d)16.108.

12. Using Lagrange’s interpolation, find the polynomial through (0,0),(1,1) and
(2,2).

13.   cosx  = _______________________.

Section-B
Subjective Questions
1. Prepare a MATLAB code to construct a forward difference table for the
given data.
2. Prepare a MATLAB code to construct a backward difference table for the
given data.
3. Certain corresponding values of x and log x are (300, 2.4771), (304,
2.4829), (305, 2.4843) and (307, 2.4871). Find log 301.
4. If the interval of differencing unity prove that

5. Find a cubic polynomial in x which takes on the values -3, 3, 11, 27, 57
and 107, when x=0, 1, 2, 3, 4 and 5 respectively.
6. Using Newton’s forward interpolation formula, for the given table of values
X 1.1 1.3 1.5 1.7 1.9
0.21 0.69 1.25 1.89 2.61
Obtain the value of f(x) when x= 1.4.

7. Develop a MATLAB code to find y value corresponding to some x value by

using Newton’s forward difference interpolation formula.
8. The population of a town in the decimal census was given below. Estimate
the population for the 1895
Year x 1891 1901 1911 1921 1931
Population of y 46 66 81 93 101

9. Find the cubic polynomial which takes the values

y(0)=1,y(1)=0,y(2)=1,y(3)=10
10. Using Newton's backward formula find the value of sin38 ?
x: 0 10 20 30 40
sin x: 0 .17365 .34202 .50000 .64279
11. Develop a MATLAB code to find y value corresponding to some x value by
using Newton’s backward difference interpolation formula.
12. Fit a polynomial of degree three which takes the following values
x: 3 4 5 6
y: 6 24 60 120
13. Using Newton’s forward formula, find the value of f(1.6) if
X 1 1.4 1.8 2.2 2.6
y 3.49 4.82 5.96 6.5 8.4

14. Find log 58.75 from the following data:

X 40 45 50 55 60 65
Log x 1.60206 1.65321 1.69897 1.74036 1.77815 1.81291
Using Newton’s Backward Interpolation formula.

15. Write a MATLAB program to implementation of Gauss backward difference

interpolation.
16. Using Gauss forward formula find y(3.3) from the following data
X 1 2 3 4 5
Y 15.30 15.10 15.0 14.5 14.0
17. Develop a MATLAB code to find y value corresponding to some x value by
using Gauss backward difference interpolation formula.
18. Write a MATLAB program to implementation of Lagrange’s interpolation.
19. Find the Lagrange’s interpolating polynomial and using it find y when x =
10, if the values of x and y are given as follows:
x 5 6 9 11
y 12 13 14 16

20. Find the number of students who got marks between 40 and 45
Marks : 30-40 40-50 50-60 60-70 70-80
No. of students : 31 42 51 35 31

21. The area A of a circle of diameter d is given below:

d: 80 85 90 95 100
A: 5026 5674 6362 7088 7854
Find approximately the areas of the circles of diameters 82 and 91.
Section-C
GATE/IES/Placement Tests/Other competitive examinations
1. Evaluate ∆10(1-x)(1-2x)(1-3x)..........(1-10x) taking h=1
 UNIT-IV
NUMERICAL DIFFERENTIATION AND INTEGRATION

Objectives:
 To understand the concepts of numerical differentiation and integration.

Syllabus:
Pre-requisite : Commands of MATLab
Approximation of derivative using Newton’s forward and backward formulae -
Integration using Trapezoidal and Simpson’s rules.

Learning Outcomes:
At the end of the unit, Students will be able to
 Utilize numerical techniques to evaluate derivatives at a point.
 Utilize numerical techniques to evaluate definite integrals.
 Calculate the area and slope of a given curve.

Introduction:

Suppose a function y = f(x) is given by a table of values (xi, yi). The

dy
process of computing the derivative for some particular value of x is called
dx
Numerical differentiation.

Derivatives using Newton’s forward difference formula:

 dy d2y
Q. Find and at x  51 from the following data.
dx dx 2
x 50 60 70 80 90
y 19.96 36.65 58.81 77.21 94.61

Solution: Here h = 10. To find the derivatives of y at x = 51 we use Newton’s

Forward difference formula taking the origin at a = 50.
x  x0 51  50
We have p   0.1
h 10

 dy 

 dy 

1
  
 2 p  1 2
 
3 p2  6 p  2 3
 
 
4 p3  18 p 2  22 p  6 4
 




    y0 y0 y0 y0 ...
 dx  x 51  dx  p 0.1 h  2! 3! 4! 


The difference table is given by

x  50
x p y y 2 y 3 y 4 y
10
50 0 19.96
16.69
60 1 36.65 5.47
22.16 -9.23
70 2 58.81 -3.76 11.99
18.40 2.76
80 3 77.21 -1.00
17.40
90 4 94.61

  4  0.13  18  0.12  22  0.1  6  

 dy   0.2  1  3  0.12  6  0.1  2     11.99  ...
 16.69 
1
   
5.47     9.23  
 dx  p 0.1 10  2

6

24 
 

1
 16.69  2.188  2.1998  1.9863  1.0316
10
 d2 y  1  2 
6 p 2  18 p  11 4  
 2    y0   p  1  3
y0   y0  ...
 dx  p 0.1 h  
2
12

 6 .12  18 .1  11 
5.47   0.1  1 9  23      11.99
1 

100  12 
 
1
 5.47  8.307  9.2523
100
=0.2303.
Q. The population of a certain town is shown in the following table
Year x 1931 1941 1951 1961 1971
Population y 40.62 60.80 79.95 103.56 132.65

Find the rate of growth of the population in 1961.

dy
Solution. Here h = 10 Since the rate of growth of population is we have to
dx
dy
find at x = 1961, which lies nearer to the end value of the table. Hence we
dx
choose the origin at x = 1971 and we use Newton’s backward interpolation
formula for derivative.

dy 1 
 y4 
 2 p  1 2
 y4 
3 p2  6 p  2 3
 y4 

2 p3  9 p 2  11 p  3 4  
 y4  ...

dx h  2 6 12 
 
x  x0 1961  1971
Where p    1
10 10
The backward difference table
x Year y Population y 2 y 3 y 4 y
1931 40.62
20.18
1941 60.80 -1.03
19.15 5.49 -4.47
1951 79.95 4.46
23.61 1.02
1961 103.56 5.48
29.09
1971 132.65

 3  12  6  1  2   2  13  9  12  11 1  3 

 dy  1  1     4.47 
   29.09      5.48    1.02   
 dx  p 1 10  2 6 12 
 
1
  29.09  2.74  0.17  0.3725
10
1
  26.5525  2.6553
10

The rate of growth of the population in the year 1961 is 2.6553.

MATLAB Code to numerical differentiation by using Newton’s forward
difference Interpolation:
function[d1val,d2val]=derbynfdi(xd,yd,xval)
n=length(xd);
syms p;
fdt=zeros(n,n+1);
fdt(:,1)=xd';
fdt(:,2)=yd';
for j=3:n+1
for i=1:n-j+2
fdt(i,j)=fdt(i+1,j-1)-fdt(i,j-1);
end
end
fdt
temp=[1,cumprod(p-(0:n-2))];
a=diff(temp);
b=diff(a);
h=xd(2)-xd(1)
p=(xval-xd(1))/h;
ap=eval(a);
bp=eval(b);
c=fdt(1,2:n+1)./factorial(0:n-1);
d1val=sum(ap.*c)/h;
d2val=sum(bp.*c)/(h*h);
fprintf('\n the approximated value of first derivative of y at given x value %f is
%f ',xval,d1val)
fprintf('\n the approximated value of second derivative of y at given x value %f
is %f ', xval,d2val)

MATLAB Code to numerical differentiation by using Newton’s backward

difference Interpolation:
function[d1val,d2val]=derbynbdi(xd,yd,xval)
n=length(xd);
syms p;
bdt=zeros(n,n+1);
bdt(:,1)=xd';
bdt(:,2)=yd';
for j=3:n+1
for i=n:-1:j-1
bdt(i,j)=bdt(i,j-1)-bdt(i-1,j-1);
end
end
bdt;
temp=[1,cumprod(p+(0:n-2))];
a=diff(temp);
b=diff(a);
h=xd(2)-xd(1);
p=(xval-xd(n))/h;
ap=eval(a);
bp=eval(b);
c=bdt(n,2:n+1)./factorial(0:n-1);
d1val=sum(ap.*c)/h;
d2val=sum(bp.*c)/(h*h);
fprintf('\n the approximated value of first derivative of y at given x value %f is
%f ',xval,d1val)
fprintf('\n the approximated value of second derivative of y at given x value %f
is %f ', xval,d2val)

Numerical Integration :
Simpson’s 1/3rd Rule:
Memorise:

Trapezoidal Rule:
b
h
 f ( x)dx  2  y
a
0  yn  2( y1  y2  y3  .......  yn 1 ) 

h
 ( sum of first and last ordinates  2( sum of remaining ordinates))
2

Simpson’s 1/3rd Rule:

b
h
 f ( x)dx  3  y
a
0  yn  2( y2  y4  y6  .......  yn 2 )  4( y1  y3  y5  .......  yn 1 ) 

h
 ( sum of first and last ordinates  2( sum of even ordinates)  4( sum of odd ordinates ))
2

Note: For Simpsons 1/3rd Rule number of subintervals should be even.

Simpson’s 3/8th Rule:
b
h
 f ( x)dx  3  y
a
0  yn  2( y3  y6  y9  .......  yn 3 )  3( y1  y2  y4  y5  .......  yn 2  yn 1 ) 

h
 ( sum of first and last ordinates  2(sum of multiples of 3 ordinates)  3( sum of remaining ordinates ))
2

Note: For Simpsons 3/8th Rule number of subintervals should be multiple of 3.

1 dx
Q .Evaluate  using (i) Trapezoidal rule (ii) Simpson’s one third rule (iii)
0 1 x

1
Simpson’s three eight rule. Take h  for all cases.
6
1 1
Solutions: Here h  , Let y  f ( x)  . The values of f(x) for the points of
6 1 x
subdivisions are as follows:
1 2 3 4 5
x 0 1
6 6 6 6 6
1
y 1 0.8571 0.75 0.6667 0.6 0.5455 0.5
1 x

(i) Tapezoidal rule

1 dx h
0 1  x  2  y0  y6   2  y1  y2  y3  y4  y5 
1
; 1  0..5   2  0.8571  0.755  0.6667  0.6  0.5455  
12 
=0.6949.

(ii) Simpson’s one third rule

dx h
 1  x  3  y  y6   2  y2  y4 )  4( y1  y3  y5  
1

0 0

1
; 1  0.5   2  0.75  0.6)  4(0.8571  0.6667  0.5455  
18 

= 0.6932.

(iii) Simpson’s three eight rule

dx 3h
 y0  y6   3  y1  y2  y4  y5   2 y3 
1
 1 x 
0 8 
 1  0.5   3  0.8571  0.75  0.6  0.5455  2  0.6667   
1
;
16  

= 0.6932.

MATLAB Code to numerical integration by using Trapezoidal Rule:

function intval=trapezoidal(f,a,b,n)
h=(b-a)/n;
x=linspace(a,b,n+1);
y=f(x)
intval=(h/2)*(y(1)+y(n+1)+2*sum(y(2:n)));
fprintf('\nthe approximate value of given integral is %f',intval)

1
MATLAB Code to numerical integration by using Simpson’s rd Rule:
3
function intval=simpsons1(f,a,b,n)
if rem(n,2)==0
h=(b-a)/n;
x=linspace(a,b,n+1);
y=f(x);
intval=(h/3)*(y(1)+y(n+1)+2*sum(y(3:2:n-1))+4*(sum(y(2:2:n))));
fprintf('\nthe approximate value of given integral is %f',intval)
else
error('number of subintervals n must be even')
end
3
MATLAB Code to numerical integration by using Simpson’s rd Rule:
8
function intval=simpsons2(f,a,b,n)
if rem(n,3)==0
h=(b-a)/n;
x=linspace(a,b,n+1);
y=f(x);
intval=(3*h/8)*(y(1)+y(n+1)+2*sum(y(4:3:n-2))+3*(sum(y)-y(1)-y(n+1)-
sum(y(4:3:n-2))));
fprintf('\n the approximate value of given integral is %f',intval)
else
error('number of subintervals n must be even')
end
 Assignment-Cum-Tutorial Questions

Section-A
Objective Questions:

1. By Newton’s forward interpolation formula

dy
 _________________________________________
dx
d2y
 _________________________________________
dx 2
2. By Newton’s backward interpolation formula
dy
 _________________________________________
dx
d2y
 _________________________________________
dx 2
3. Trapezoidal rule to find definite integral is_____________________
4. Simpson’s 1/3rd rule to find definite integral is_____________________
5. Simpson’s 3/8th rule to find definite integral is_____________________
6. If we put n = 2 in a general quadrature formula, we get [ ]
(a) Trapezoidal rule (b) Simpson’s 1/3rd rule
(c) Simpson’s 3/8th rule (d) Boole’s rule
7. In Simpson’s 1/3 rule the number of subintervals should be
rd [ ]
(a) Even (b) odd
(c) multiples of 3’s (d) more than ‘n’ interval
8. If the distance d(t) is traversed by a particle in the ‘t’ sec and d(0) = 0, d(2)
= 8, d(4) = 20 and d(6) = 28, then its velocity in cm after 6 secs is [ ]
(a) 1.67 (b) 16.67 (c) 2 (d) 2.003
1 1 1 
9. The formula y0  2 y0  3 y0  .......  is used only when the point x is
h 2 3 
at [ ]
(a) end of the tabulated set (b) middle of the tabulated set
(c) Beginning of tabulated set (d) none of these
10. To increase the accuracy in evaluating a definite integral by Trapezoidal
rule, we should take ______________________
11. Values of y = f(x) are known as x = x0, x1 and x2. Using Newton’s forward
 dy 
integration formula, the approximate value of   is
 dx  x x0
________________________________
12. Numerical differentiation gives [ ]
(a) exact value (b) approximate value
 (c) no result (d) negative value

13. The general quadrature formula is [ ]

(a) always same (b)depends upon interpolation formula
(c) not easy to derive (d) is also given approximate result
x1

14. For n =1in quadrature formula,  f ( x) dx

x0
equals to [ ]

 f 0  f1  (b)  f 0  f1  (c)  f 0  f1  (d)  f 0  f1 

h h h
(a)
2 2 4
15. To apply, Simpson’s 1/3 rd rule, always divide the given range of
integration into ‘n’ subintervals, where n is [ ]
(a) even (b) odd (c) 1,2,3,4 (d) 1,3,5,7
16. The process of calculating derivative of a function at some particular
value of the independent variable by means of a set of given values of
that function is [ ]
(a) Numerical value (b) Numerical differentiation
(c) Numerical integration (d) quadrature
17. While evaluating definite integral by Trapezoidal rule, the accuracy can
be increased by [ ]
(a) h =4 (b) even number of sub-intervals
(c) multiples of 3 (d) large number of sub-intervals

Section-B
Subjective Questions

1. A curve is expressed by the following values of x and y. Find the slope at

the point x =1.5

x 0.0 0.5 1.0 1.5 2.0

y 0.4 0.35 0.24 0.13 0.05
2. The population of a certain town is given below. Find the rate of growth of
the population in 1961:

Year 1931 1941 1951 1961 1971

Population 40.62 60.80 71.95 103.56 132.65
3. In a machine a slider moves along a fixed straight rod. Its distance x cms
along the rod is given below for various values of time ‘t’ seconds. Find the
velocity and acceleration of the slider when t = 0.3

t(sec) 0 0.1 0.2 0.3 0.4 0.5 0.6

x(cms) 30.13 31.62 32.87 33.64 33.95 33.81 33.24
4. The velocity of a train which starts from rest is given by the following
table being reckoned in minutes from the start and speed in miles per
hour
 Minutes 2 4 6 8 10 12 14 16 18
Miles per 10 18 25 29 32 20 11 5 2
hour

Estimate approximately the total distance travelled in 20 minutes.

5. The distance covered by an athlete for the 50 meter is given in the
following table

Time(sec) 0 1 2 3 4 5 6
Distance(meter) 0 2.5 8.5 15.5 24.5 36.5 50
Determine the speed of the athlete at t = 5 sec. correct to two decimals.
6. A curve is drawn to pass through the points given by following table:
x 1 1.5 2.0 2.5 3 3.5 4.0
y 2 2.4 2.7 2.8 3 2.6 2.1
Find the slope of the curve at x=1.25.
7. Write a MATLAB code to find derivatives at a point by using Newtons
forward interpolation formula.
8. Prepare a MATLAB code to find derivatives at a point by using Newtons
backward interpolation formula.
2

9. Evaluate  e
x 3
dx using Simpson’s rule taking h=0.25
0
10. A river is 80 meters wide. The depth 'd’ in meters at a distance x from the
bank is given in the following table. Calculate the cross section of the
river using Trapizoidal rule.
x 10 20 30 40 50 60 70 80

d(x) 4 7 9 12 15 14 8 3
5.2 5.2
11. Compute the value of the definite integral  log xdx or  ln xdx
4 4
using

i.Trapezoidal Rule ii. Simpson’s 1/3rd Rule and iii. Simpson’s 3/8th Rule.
12. Create a MATLAB function file to find a definite integral by using
Trapezoidal rule.
13. Design a MATLAB code to find a definite integral by using Simpson’s
1/3rd rule.
14. Write a MATLAB function file to find a definite integral by using
Simpson’s 3/8th rule.

Section-C
GATE/IES/Placement Tests/Other competitive examinations

1. If f(2) = 5, f(4) = 8, f(6) = 10, and f(8) = 16 then

 1 .4
2. Using Simpson’s 1/3rd rule, find the value of the integr  (sin x  log x  e x ) dx
0 .2

by taking 6 sub-intervals.
2
1
3. Minimum number of subintervals required to evaluate the integral  x dx
1

by using Simpson’s 1/3rd

rule so that the value is corrected up to 4
decimal places.
4. The following table gives the velocity v of a particle at time ‘t’

t (seconds) 0 2 4 6 8 10 12
v meters per 4 6 16 34 60 94 136
second

Find (i) the distance moved by the particle in 12 seconds and also (ii) the
acceleration at t = 2 sec.

5. A rocket is launched from the ground. Its acceleration is registered

during the first 80 seconds and is given in the table below. Using
Simpson’s 1/3rd rule, find the velocity of the rocket at t= 80 seconds.

t sec 0 10 20 30 40 50 60 70 80
f(cm/sec2) 30 31.63 33.34 35.47 37.75 40.33 43.25 46.69 50.67

-----------------@-@-@-@-@----------------
 UNIT–V: Numerical Solutions of Ordinary Differential Equations
Pre-requisite : Commands of MATLab

Objectives:
• To the numerical solutions of a first ordered Ordinary Differential Equation together with initial
condition.

Syllabus:
Taylor’s Series Method - Euler Method - Modified Euler Method - Runge – Kutta Fourth order Method.

Subject Outcomes:
At the end of the unit, Students will be able to

• solve Ordinary Differential equations using Numerical methods.

• implement the numerical methods in Matlab in order to solve Ordinary Differential equations.

The important methods of solving ordinary differential equations of first order numerically are as
follows

 Taylors series method

 Euler’s method

 Modified Euler’s method of successive approximations

 Runge- kutta method

To describe various numerical methods for the solution of ordinary differential equations we consider
the general 1st order differential equation

dy/dx = f(x,y) ------- (1)

with the initial condition y(x0) = y0

The methods will yield the solution in one of the two forms:

i) A series for y in terms of powers of x, ,from which the value of y can be obtained by direct
substitution.
ii ) A set of tabulated values of y corresponding to different values of x.

TAYLOR’S SERIES METHOD

To find the numerical solution of the differential equation

dy
 f ( x, y ) (1)
dx

With the initial condition y ( x0 )  y0  (2)

y ( x ) can be expanded about the point x0 in a Taylor’s series in powers of ( x  x0 ) as

( x  x0 ) ( x  x0 ) 2 ( x  x0 ) n n
y ( x)  y ( x0 )  y( x0 )  y( x0 )  ............  y ( x0 ) (3)
1 2! n!

In equation3, y ( x0 ) is known from Initial Condition . The remaining coefficients

y( x0 ), y( x0 ),......... y n ( x0 ) etc are obtained by successively differentiating equation1 and evaluating at
x0 . Substituting these values in equation3, y ( x ) at any point can be calculated from equation3.
Provided h  x  x0 is small.

When x0  0 , then Taylor’s series equ3 can be written as

x2 xn n
y ( x)  y (0)  x. y(0)  y(0)  ......  y (0)  ........ (4)
2! n!

Note: We know that the Taylor’s expansion of y(x) about the point x0 in a power of (x – x0)is.

( x  x0 ) I ( x  x0 ) 2 II ( x  x0 )3 III
y(x) = y(x0) + y (x0) + y (x0) + y (x0) + … (1)
1! 2! 3!

Or

( x  x0 ) I ( x  x0 ) 2 II ( x  x0 )3 III
y(x) = y0 + y0 + y0 + y0 + …..
1! 2! 3!

If we let x – x0 = h. (i.e. x = x0 + h = x1) we can write the Taylor’s series as

h I h 2 II h3 III h 4 IV
y(x) = y(x1) = y0 + y + y0 + y0 + y0 + ….
1! 0 2! 3! 4!

h I h 2 II h3 III h IV IV
i.e. y1 = y0 + y + y0 + y0 + y0 + ….. (2)
1! 0 2! 3! 4!
Similarly expanding y(x) in a Taylor’s series about x = x1. We will get.

h I h 2 II h3 III h 4 IV
y2 = y1 + y1 + y1 + y1 + y1 + ……. (3)
1! 2! 3! 4!

Similarly expanding y(x) in a Taylor’s series about x = x2 We will get.

h I h 2 II h3 III h 4 IV
y 3 = y2 + y2 + y2 + y2 + y2 + …... (4)
1! 2! 3! 4!

In general, Taylor’s expansion of y(x) at a point x= x n is

h I h 2 II h3 III h 4 IV
yn+1 = yn + yn + yn + yn + yn + ….. (5)
1! 2! 3! 4!

Example 1. Using Taylor’s expansion evaluate the integral of y  2 y  3e x , y (0)  0 at x  0.2 . Hence
compare the numerical solution obtained with exact solution .

Sol: Given equation can be written as 2 y  3e x  y, y (0)  0

Differentiating repeatedly w.r.t to ‘x’ and evaluating at x  0

y( x)  2 y  3e x , y(0)  2 y(0)  3e0  2(0)  3(1)  3

y( x)  2 y  3e x , y(0)  2 y(0)  3e0  2(3)  3  9
y( x)  2. y( x)  3e x , y(0)  2 y(0)  3e 0  2(9)  3  21
y iv ( x)  2. y( x)  3e x , y iv (0)  2(21)  3e0  45
y v ( x)  2. y iv  3e x , y v (0)  2(45)  3e0  90  3  93

In general, y ( n 1) ( x)  2. y ( n ) ( x)  3e x or y ( n 1) (0)  2. y ( n ) (0)  3e0

The Taylor’s series expansion of y ( x ) about x0  0 is

x2 x3 x4 x5
y ( x)  y (0)  xy(0)  y(0)  y(0)  y(0)  y(0)  ....
2! 3! 4! 5!

Substituting the values of y (0), y(0), y(0), y(0),..........

9 2 21 3 45 4 93 5
y ( x)  0  3 x  x  x  x  x  ........
2 6 24 120

9 2 7 3 15 4 31 5
y ( x)  3x  x  x  x  x  ........  equ1
2 2 8 40
Now put x  0.1 in equ1

9 7 15 31
y (0.1)  3(0.1)  (0.1) 2  (0.1)3  (0.1) 4  (0.1) 5  0.34869
2 2 8 40

Now put x  0.2 in equ1

9 7 15 31
y (0.2)  3(0.2)  (0.2) 2  (0.2) 3  (0.2) 4  (0.2) 5  0.811244
2 2 8 40

9 7 15 31
y (0.3)  3(0.3)  (0.3) 2  (0.3) 3  (0.3) 4  (0.3) 5  1.41657075
2 2 8 40

Analytical Solution:

dy
The exact solution of the equation  2 y  3e x with y (0)  0 can be found as follows
dx

dy
 2 y  3e x Which is a linear in y.
dx

Here P  2, Q  3e x

pdx 2 dx
I.F = 
e

e
 e 24


2 x x 2 x x
General solution is y.e  3e .e dx  c  3e  c

 y  3e x  ce 2 x where x  0, y  0 0  3  c  c  3

The particular solution is y  3e2 x  3e x or y ( x)  3e 2 x  3e x

Put x  0.1 in the above particular solution,

y  3.e0.2  3e0.1  0.34869

Similarly put x  0.2 , y  3e0.4  3e0.2  0.811265

put x  0.3 , y  3e0.6  3e0.3  1.416577

MATLAB CODE FOR THE IMPLEMENTATION OF TAYLOR SERIES METHOD:-

To find the numerical solution of the differential equation

 dy
 f ( x, y ) with the initial condition y ( x0 )  y0
dx

function [ soltable ] = odetaylor( f,x0,y0,xn,h,not )

% here f is a symbolic function function f(x,y) such that y=y(x)

syms x y(x);

xd=x0:h:xn;

n=length(xd);

yd=zeros(1,n);

yd(1)=y0;

for i=1:n-1

fold=f;

d=ones(1,not-1);

for j=1:not-1

d(j)=subs(fold,{x,y},{xd(i),yd(i)});

fold=subs(diff(fold,x),diff(y(x),x),fold);

end

yd(i+1)=yd(i)+sum(d.*(h.^(1:not-1)./factorial(1:not-1)));

end

soltable=[xd' yd'];

end

EULER’S METHOD

It is the simplest one-step method and it is less accurate. Hence it has a limited application.
 dy
Consider the differential equation = f(x,y) (1)
dx

With y(x0) = y0 (2)

Consider the first two terms of the Taylor’s expansion of y(x) at x = x0

y(x) = y(x0) + (x – x0) y1(x0) (3)

from equation (1) y1(x0) = f(x0,y(x0)) = f(x0,y0)

Substituting in equation (3)

y(x) = y(x0) + (x – x0) f(x0,y0)

At x = x1, y(x1) = y(x0) + (x1 – x0) f(x0,y0)

y1 = y0 + h f(x0,y0) where h = x1 – x0

Similarly at x = x2 , y2 = y1 + h f(x1,y1),

Proceeding as above, yn+1 = yn + h f(xn,yn)

This is known as Euler’s Method

dy
Example 1. Using Euler’s method solve for x = 2 from = 3x2 + 1,y(1) = 2,taking step size
dx
(I) h = 0.5 and (II) h=0.25

Sol: Here f(x,y) = 3x2 + 1, x0 = 1,y0 = 2

Euler’s algorithm is yn+1 = yn + h f(xn,yn), n = 0,1,2,3,….. (1)

h = 0.5 x1 = x0 + h = 1+0.5 = 1.5

Taking n = 0 in (1) , we have x2 = x1 + h = 1.5 + 0.5 = 2

y1 = y0 + h f(x0,y0)

i.e. y1 = y(0.5) = 2 + (0.5) f(1,2) = 2 + (0.5) (3 + 1) = 2 + (0.5)(4)

Here x1 = x0 + h = 1 + 0.5 = 1.5

y(1.5) = 4 = y1

Taking n = 1 in (1),we have

y2 = y1 + h f(x1,y1)

i.e. y(x2) = y2 = 4 + (0.5) f(1.5,4) = 4 + (0.5)[3(1.5)2 + 1] = 7.875

 Here x2 = x4 + h = 1.5 + 0.5 = 2

y(2) = 7.875
h = 0.25 x1 = 1.25, x2 = 1.50, x3 = 1.75, x4 = 2
Taking n = 0 in (1), we have

y1 = y0 + h f(x0,y0)

i.e. y(x1) = y1 = 2 + (0.25) f(1,2) = 2 + (0.25) (3 + 1) = 3

y(x2) = y2 = y1 + h f(x1,y1)

i.e. y(x2) = y2 = 3 + (0.25) f(1.25,3)

= 3 + (0.25)[3(1.25)2 + 1]

= 4.42188

Here x2 = x1 + h = 1.25 + 0.25 = 1.5

y(1.5) = 5.42188

Taking n = 2 in (1), we have

i.e. y(x3) = y3 = h f(x2,y2)

= 5.42188 + (0.25) f(1.5,2)

= 5.42188 + (0.25) [3(1.5)2 + 1]

= 6.35938

Here x3 = x2 + h = 1.5 + 0.25 = 1.75

y(1.75) =7. 35938

Taking n = 4 in (1),we have

y(x4) = y4 = y3 + h f(x3,y3)

i.e. y(x4) = y4 = 7.35938 + (0.25) f(1.75,2)

= 7.35938 + (0.25)[3(1.75)2 + 1]

= 8.90626
 Note that the difference in values of y(2) in both cases
(i.e. when h = 0.5 and when h = 0.25).The accuracy is improved significantly when h is reduced to 0.25
(Exact solution of the equation is y = x3 + x and with this y(2) = y2 = 10.

MATLAB CODE FOR THE IMPLEMENTATION OF EULER’S METHOD:-

To find the numerical solution of the differential equation

dy
 f ( x, y ) with the initial condition y ( x0 )  y0
dx

function [ soltable ] = odeeuler( f,x0,y0,xn,h )

x=x0:h:xn;

n=length(x);

y=zeros(1,n);

y(1)=y0;

for i=1:n-1

y(i+1)=y(i)+h*f(x(i),y(i));

end

soltable=[x' y'];

end

Modified Euler’s method

 i 1
It is given by y
i 
 yk  h / 2 f  xk , yk   f  xk 1 ,1k 1  , i  1, 2....., ki  0,1.....
k 1
 

Working rule :

i) Modified Euler’s method

  i 1
y i k 1  yk  h / 2 f  xk , yk   f  xk 1 ,1k 1  , i  1, 2....., ki  0,1.....
 

ii) When i  1 y 0k 1 can be calculated from Euler’s method

iii) k=0, 1……… gives number of iteration. i  1, 2...

gives number of times, a particular iteration k is repeated

Suppose consider dy/dx=f(x, y) -------- (1) with y(x0) =y0----------- (2)2

To find y(x1) =y1 at x=x1=x0+h

Now take k=0 in modified Euler’s method

We get y1
1
   
 y0  h / 2  f  x0 , y0   f x1 , y1i 1  ……………………… (3)

Taking i=1, 2, 3...k+1 in equation (3), we get

y1 0  y0  h / 2  f  x0 , y0  (By Euler’s method)


y11  y0  h / 2  f  x0 , y0   f x1 , y1 0 
  

y1 2  y0  h / 2  f  x0 , y0   f x1 , y11 
  
------------------------


y1 k 1  y0  h / 2  f  x0 , y0   f x1 , y1 k  
  
k   k 1
If two successive values of y1 , y1 are sufficiently close to one another, we will take the common
value as y2  y  x2   y  x1  h 

We use the above procedure again

Example1. Using modified Euler’s method, find the approximate value of x when x  0.3 given that
dy / dx  x  y and y  0   1

Sol: Given dy / dx  x  y and y  0   1

Here f  x, y   x  y, x0  0, and y0  1

Take h = 0.1 which is sufficiently small

Here x0  0, x1  x0  h  0.1, x2  x1  h  0.2, x3  x2  h  0.3

The formula for modified Euler’s method is given by


yk 1i   yk  h / 2  f  xk  yk   f xk 1 , yk 1i 1   1
  
Step1: To find y1= y(x1) = y (0.1)

Taking k = 0 in eqn(1)


yk 1i   y0  h / 2  f  x0  y0   f x1 , y1i 1    2 
  
when i = 1 in eqn (2)


y1i   y0  h / 2  f  x0 , y0   f x1 , y1 0 
  
(0)
First apply Euler’s method to calculate y 1
= y1

 y1 0  y0  h f  x0 , y0 

= 1+(0.1)f(0.1)

= 1+(0.1)

= 1.10

now  x0  0, y0  1, x1  0.1, y1  0   1.10

1
 
 y1   y0  0.1/ 2  f  x0 , y0   f x1 , y1  
0

= 1+0.1/2[f(0,1) + f(0.1,1.10)

= 1+0.1/2[(0+1)+(0.1+1.10)]

= 1.11

When i=2 in equation (2)


y1 2  y0  h / 2  f  x0 , y0   f x1 , y11 
  
= 1+0.1/2[f(0.1)+f(0.1,1.11)]

= 1 + 0.1/2[(0+1)+(0.1+1.11)]
 = 1.1105

 
y13  y0  h / 2  f  x0 , y0   f x1 , y1 2 
 

= 1+0.1/2[f(0,1)+f(0.1 , 1.1105)]

= 1+0.1/2[(0+1)+(0.1+1.1105)]

= 1.1105

 2
Since y1  y13

 y1 = 1.1105

Step:2 To find y2 = y(x2) = y(0.2)

Taking k = 1 in equation (1) , we get

 
y2i   y1  h / 2  f  x1 , y1   f x2 , y2i 1    3
 

I = 1,2,3,4,…..

For i = 1

 
y21  y1  h / 2  f  x1 , y1   f x2 , y2 0 
 

y2 0 is to be calculate from Euler’s method

y2 0  y1  h f  x1 , y1 

= 1.1105 + (0.1) f(0.1 , 1.1105)

= 1.1105+(0.1)[0.1+1.1105]

= 1.2316

= 1.1105  0.1/ 2  f  0.1,1.1105   f  0.2,1.2316 

(1)
 y 2

= 1.1105 +0.1/2[0.1+1.1105+0.2+1.2316]

= 1.2426
 
y2 2  y1  h / 2  f  x1 , y1   f x2 y2 1 
 
= 1.1105 + 0.1/2[f(0.1 , 1.1105) , f(0.2 . 1.2426)]

= 1.1105 + 0.1/2[1.2105 + 1.4426]

= 1.1105 + 0.1(1.3266)

= 1.2432


y23  y1  h / 2  f  x1 , y1   f x2 y2 2 
  
= 1.1105+0.1/2[f(0.1,1.1105)+f(0.2 , 1.2432)]

= 1.1105+0.1/2[1.2105+1.4432)]

= 1.1105 + 0.1(1.3268)

= 1.2432

 3
Since y2  y23

Hence y2 = 1.2432

Step:3

To find y3 = y(x3) = y y(0.3)

Taking k =2 in equation (1) we get

 
y31  y2  h / 2  f  x2 , y2   f x3 , y3i 1    4 
 

For i = 1 ,


y31  y2  h / 2  f  x2 , y2   f x3 , y3 0 
  
y3 0 is to be evaluated from Euler’s method .

y3 0  y2  h f  x2 , y2 

= 1.2432 +(0.1) f(0.2 , 1.2432)

= 1.2432+(0.1)(1.4432)
 = 1.3875

 y31 = 1.2432+0.1/2[f(0.2 , 1.2432)+f(0.3, 1.3875)]

= 1.2432 + 0.1/2[1.4432+1.6875]

= 1.2432+0.1(1.5654)

= 1.3997


y3 2  y2  h / 2  f  x2 , y2   f x3 , y31 
 
= 1.2432+0.1/2[1.4432+(0.3+1.3997)]

= 1.2432+ (0.1) (1.575)

= 1.4003

 
y33  y2  h / 2  f  x2 , y2   f x3 , y3 2 
 

= 1.2432+0.1/2[f(0.2 , 1.2432)+f(0.3 , 1.4003)]

= 1.2432 + 0.1(1.5718)

= 1.4004

 
y3 4  y2  h / 2  f  x2 , y2   f x3 , y33 
 

= 1.2432 + 0.1/2[1.4432+1.7004]

= 1.2432+(0.1)(1.5718)

= 1.4004

 3
Since y3  y3 4

 The value of y at x = 0.3 is 1.4004

MATLAB CODE FOR THE IMPLEMENTATION OF Modified Euler’s method METHOD:-

To find the numerical solution of the differential equation

 dy
 f ( x, y ) with the initial condition y ( x0 )  y0
dx

function [ soltable ] = odemodifiedeuler( f,x0,y0,xn,h,tol )

x=x0:h:xn;

n=length(x);

y=zeros(1,n);

y(1)=y0;

for i=1:n-1

yp=y(i)+h*f(x(i),y(i));

y(i+1)=y(i)+h*(f(x(i),y(i))+f(x(i),yp))/2;

while abs(yp-y(i+1)>tol)

yp=y(i+1);

y(i+1)=y(i)+h*(f(x(i),y(i))+f(x(i),yp))/2;

end

end

soltable=[x' y'];

end

Runge – Kutta Methods

I. Second order R-K Formula

yi+1 = yi+1/2 (K1+K2),

Where K1 = h (xi, yi)

K2 = h (xi+h, yi+k1) for i= 0,1,2-------

II. Third order R-K Formula

yi+1 = yi+1/6 (K1+4K2+ K3),

Where K1 = h (xi, yi)

 K2 = h (xi+h/2, y0+k1/2)

K3 = h (xi+h, yi+2k2-k1)

For i= 0,1,2-------

III. Fourth order R-K Formula

yi+1 = yi+1/6 (K1+2K2+ 2K3 +K4),

Where K1 = h (xi, yi)

K2 = h (xi+h/2, yi+k1/2)

K3 = h (xi+h/2, yi+k2/2)

K4 = h (xi+h, yi+k3)

For i= 0,1,2-------

Example 1. Apply the 4th order R-K method to find an approximate value of y when x=1.2 in stepsof 0.1,
given that y1 = x2+y2, y (1)=1.5

sol. Given y1= x2+y2,and y(1)=1.5

Here f(x,y)= x2+y2, y0 =1.5 and x0=1,h=0.1

So that x1=1.1 and x2=1.2

Step1:

To find y1 :

By 4th order R-K method we have

y1=y0+1/6 (k1+2k2+2k3+k4)

k1= h f(x0,y0) = (0.1) f(1,1.5) = (0.1) [12+(1.5)2] = 0.325

k2= hf (x0+h/2,y0+k1/2) = (0.1) f(1+0.05,1.5+0.325) = 0.3866

and

k3= h f((x0+h/2,y0+k2/2) = (0.1) f(1.05,1.5+0. 3866/2) = (0.1)[(1.05)2+(1.6933)2]

= 0.39698

k4= hf(x0+h,y0+k3) = (0.1)f(1.0,1.89698) = 0.48085

Hence
 1
y1  1.5  0.325  2  0.3866   2  0.39698   0.48085 
6
 1.8955

Step2:

To find y2, i.e., y  x2   y 1.2 

Here x1=0.1,y1=1.8955 and h=0.1

by 4th order R-K method we have

y2 = y1+(1/6) (k1+2k2+2k3+k4)

k1= hf(x1,y1) = (0.1)f(0.1,1.8955) = (0.1) [12+(1.8955)2] = 0.48029

k2 = hf (x1+h/2,y1+k1/2) = (0.1)f(1.1+0.1,1.8937+0.4796) = 0.58834

k3 = hf((x1+h/2,y1+k2/2) = (0.1)f(1.5,1.8937+0.58743) = (0.1)[(1.05)2+(1.6933)2]

= 0.611715

k4 = hf(x1+h,y1+k3) = (0.1)f(1.2,1.8937+0.610728) = 0.77261

Hence

y2 = 1.8937+(1/6) (0.4796+2(0.58834)+2(0.611715)+0.7726) = 2.5043

 y =2.5043 where x  0.2 .

MATLAB CODE FOR THE IMPLEMENTATION OF 4th order R-K METHOD:-

To find the numerical solution of the differential equation

dy
 f ( x, y ) with the initial condition y ( x0 )  y0
dx

function [ soltable ] = odeRK( f,x0,y0,xn,h )

n=length(x);

y=zeros(1,n);

y(1)=y0;

for i=1:n-1
 k1=h*f(x(i),y(i));

k2=h*f(x(i)+h/2,y(i)+k1/2);

k3=h*f(x(i)+h/2,y(i)+k2/2);

k4=h*f(x(i)+h,y(i)+k3);

y(i+1)=y(i)+(k1+2*(k2+k3)+k4)/6;

end

soltable=[x' y'];

end

Assignment-cum-Tutorial Questions

A. Questions testing the remembering / understanding level of students

I Objective Questions
dy
1. If = f(x,y),y(x0) = y0,the formula for fourth order Runge – Kutta method is _______
dx

2. In which of the following methods, successive approximations are used?

a) Picard’s method b) Taylor series method

c) Adams-Bashforth method d) None of these

3. Which among the following is the self-starting method?

a) Adams-Bashforth method b) Milne’s method

c) Runge-Kutta method d) Predictor method

4. Among the following, which is the best for solving initial value problem?

a) Modified Euler’s method b) Picard’s method

c) Runge-Kutta method of fourth order d) Taylor series method

5. Which of the following is a step-by-step method?

a) Picard’s method b) Taylor series method

c) Adams-Bashforth method d) None of these

B. Questions testing the ability of students in applying the
concepts

I) Multiple choice Questions

dy
1. If = -y, y(0) = 1, h = 0.01 then by Euler’s method, the value of y1 = __________
dx

a) 0.099 b) 0.0981 c) 0.99 d) none

dy y  x
2. If = , y(0) = 1 and h = 0.02, using Euler’s method the value of y1 = _________
dx y  x

a) 1.02 b) 1.04 c) 1.03 d) none

dy
3. If = x2+y2, y(0) = 0 using Taylor’s series method, the value of y(0.4)=___________
dx

a) 0.2133 b) 0.02133 c) 0.002133 d) None

4. The value of y at x= 0.1 using Runge – Kutta method of fourth order for the differential

dy
equation = x – 2y, y(0) = 1 taking h = 0.1 is __________
dx

a) 0.825 b) 0.0825 c) 0.813 d) None

5. The value of y at x = 0.1 using modified Euler’s method up to second approximation for

dy/dx = x – y, y(0) = 1 is _________

a) 0.909 b) 0.0909 c) 0.809 d) None

dy
6 If = 1 + y2, f(x0,y0) =1,h – 0.2,K1 = 0.2,K2 = 0.202,K3 = 0.20204,K4 = 0.20216, then the
dx

value of y1 by fourth order Runge – Kutta method is ________

a) 0.0202 b) 0.202 c) 0.102 d) None

dy
7. Using Runge – Kutta method, the approximate value of y(0.1) if = x + y2,y = 1
dx
 where x =0 and f(x0,y0) = 1 K1 = 0.1, K2 = 0.115, K3 = 0.116, K4 = 0.134 is

a) 1.116 b) 1.001 c) 1.211 d) None

dy
8. Using Runge-kutta method, to solve the differential equation  x  y , h=0.1 and
dx

y(0) = 1.

(i) The values of k1, k2, k3 and k4 respectively are

a) 0.11, 0.121, 0.1, 0.005

b) 0.1, 0.11, 0.1105, 0.12105
c) 0.111, 0.11105, 0.121005, 0.121
d) None of these
(ii) For the above problem y(0.1) = _______

a) 1.11034 b) 1.33011 c) 1.43001 d) None of these

II Problems

1. Solve y1 = x-y2, y(0) = 1 using Taylor’s series method and evaluate y(0.1), y(0.2).
2. Given y1 = x+ sin y, y(0)=1 compute y(0.2) and y(0.4) with h=0.2 using modified Euler’s
method

3. Using R-K method, find y(0.2) for the equation dy/dx=y-x , y(0)=1,take h=0.22.
dy
given that y=1 when x=0 and  x y
dx

dy
4. Using Taylor’s series method, solve the equation  x 2  y 2 for x  0.4 given that y0
dx
when x  0 4.
dy
5. Find the solution of = x-y , y(0)=1 at x =0.1 , 0.2 ,0.3 , 0.4 and 0.5 using modified
dx
Euler’s method.

6. Write the R-K method of 4th order formula for the solution of

7. Using R-K method, estimate y(0.2) and y(0.4) for the equation dy/dx=y2-x2/ y2+x2,y(0)=1,h=0.2.
8. Explain the Modified Euler's method.
 dy
9. Find y(0.1) & y(0.2) using Euler's modified formula given that dx = x2 – y, y(0) = 1
10. Find y(0.1) & y(0.2) using Runge - Kutta 4th order formula given that y1 = x2 – y & y(0) =1.
11. Write a code in MATLAB to find the numerical solution of first order ordinary differential equation
using R-K Method of fourth order.
12. Write a code in MATLAB to find the numerical solution of first order ordinary differential equation
using Euler’s Method.

C. Questions testing the analyzing / evaluating ability of students

dy
1. Solve by Taylor’s series method, the equation  log xy for y(1.1) and y(1.2), given
dx
y(1) = 2.
dy 2 xy  e x
2. Using R – K method, solve for y at x = 1.2, 1.4 from  given y(1) = 0.
dx x 2  xe x

GATE QUESTIONS

1. Consider an ordinary differential equation dx/dt=4t+4 if x =0 at t = 0, the increment

in x calculated using Runge-Kutta fourth order multi-step method with a step size of
Δt = 0.2 is __________ (GATE2014)
(A) 0.22 (B) 0.44 (C) 0.66 (D) 0.88
 UNIT-VI
CURVE FITTING

Objectives:
 To understand the application of 'Least Square Method'

Pre-requisite : Commands of MATLab

Syllabus: Fitting a straight line - Parabolic curve - exponential curve - power

curve by the method of least squares.

Learning Outcomes:
At the end of the unit, Students will be able to
 Fit a straight line to the given data.
 Fit a Parabola to the given data.
 Fit an exponential curve and a power curve the given data.

Least Square Method:

The principle of least squares is one of the popular methods for finding a curve fitting a given
data. Say be n observations from an experiment. We are interested
in finding a curve

Closely fitting the given data of size 'n'. Now at while the observed value of y is , the
expected value of y from the curve (1) is . Let us define the residual by

Likewise, the residuals at all other points , ..., are given by

....................(3)

...........................
Some of the residuals may be positive and some may be negative. We would like to find the
curve fitting the given data such that the residual at any is as small as possible. Now since
some of the residuals are positive and others are negative and as we would like to give equal
importance to all the residuals it is desirable to consider sum of the squares of these residuals,
say and thereby find the curve that minimizes . Thus, we consider

and find the best representative curve (1) that minimizes (4).

2.4.2 Least Square Fit of a Straight Line

Suppose that we are given a data set of 'n' observations from an

experiment. Say that we are interested in fitting a straight line

to the given data. Find the 'n' residuals 's by:

Now consider the sum of the squares of 's i.e

Note that is a function of parameters a and b. We need to find a,b such that is minimum.
The necessary condition for to be minimum is given by:

The condition yields:

i.e

Similarly the condition yields

Equations (5) and (6) are called as normal equations,which are to be solved to get desired values
for a and b.

The expression for i.e (3) can be re-written in a convenient way as follows:

Example: Using the method of least squares, find an equation of the form

that fits the following data:

Solution: Consider the normal equations of least square fit of a straight line i.e

Here n=5.
From the given data, we have,

Therefore the normal equations are given by:

30a +10b =243 ................(3)

10a+5b=76.......................(4)

On solving (3) and (4) we get

a = 9.1 , b= - 3 ................................................................(5)

Hence the required fit for the given data is

y=9.1x - 3 ...................... ..(6)

Remarks:

(1) Experimental data may not be always linear. One may be interested in fitting either a curve of
the form However, both of these forms can be linearized by
taking logarithms on both the sides. Let us look at the details:
On taking logarithms on both the sides we get:

Say

Using (3) in (2) we get

which is linear in X, Y.

On taking logarithms we get

we get

which is linear in Y, x.

Example: By the method of least square fit a curve of the form to the following data:

Solution.

Consider ----(1)

On taking logarithm on both the sides we get

Using (3) in (2) we get

Data in modified variables X,Y

Normal equations corresponding to the straight line fit (4) are:

From the modified data we get

normal equations take the form:

On solving for b & A, we obtain,

The desired curve is

Matlab code to fit exponential curve of the form y=a*e^(bx)

function [a,b]=exp1fit(x,y) % to fit y=a*e^(bx)

% on applying logarithm ln, we have lny=lna+bx

n=length(x);

A=[sum(x) n;sum(x.*x) sum(x)];

B=[sum(ln(y));sum(ln(y).*x)];

coef=A\B;

a=exp(coef(2)); % since coef(2)=lna

b=coef(1);

end

Matlab code to fit exponential curve of the form y=a*b^x

function [a,b]=exp2fit(x,y) % to fit y=a*(b^x)

% on applying log10 we have logy=loga+xlogb

n=length(x);
A=[sum(x) n;sum(x.*x) sum(x)];

B=[sum(log10(y));sum(log10(y).*x)];

coef=A\B;

a=10^coef(2);

b=10^coef(1);

end

Matlab code to fit power curve of the form y=a*x^b

function [a,b]=powerfit(x,y) % to fit y=a*(x^b)

% on applying logarithm ln, we have lny=lna+blnx

n=length(x);

A=[sum(ln(x)) n;sum(ln(x).*ln(x)) sum(ln(x))];

B=[sum(ln(y));sum(ln(y).*ln(x))];

coef=A\B;

a=exp(coef(2)); % since coef(2)=lna

b=coef(1);

end

Least Square fit of a parabola

Given a data set of n observations of an experiment .Now we try to

fit a best possible parabola
following the principle of least square. Finding the appropriate parabola amounts to determining
the constants a,b,c that minimize the sum of the squares of the residuals given
by

The necessary condition for E to be minimum is

Now the condition yields

i.e

Similarly yields

i.e

Finally yields
Equations (4), (5) and (6) are called as normal equations whose solution yields the values of the
constants a, b and c and thus the desired parabola.

Example: Given the following data from an experimental observation

y: 9.4 11.8 14.7 18.0 23.0

x: 1.0 1.6 2.5 4.0 6.0

fit a parabola in the form following the principle of least square.

Solution) Here n=5

The normal equations for finding a parabolic fit are:

(1)
 The normal equations are:

(2)

On Solving (2) for a,b,c we get

Matlab code to fit quadratic curve i.e., parabola y=ax^2+bx+c :-

function [coef]=parabolafit(x,y)

% this is the function file to fit parabola whose

% equation is y=ax^2+bx+c

n=length(x);

% Coefficient matrix of the normal equations to fit y=ax^2+bx+c

A=[sum(x.^2) sum(x) n;

sum(x.^3) sum(x.^2) sum(x);

sum(x.^4) sum(x.^3) sum(x.^2)];

B=[sum(y);sum(y.*x);sum(y.*(x.^2))];

disp([A,B]);

coef=A\B;

end

Assignment-Cum-Tutorial Questions

Section-A
Objective Questions:
1. Write the normal equation in method of least squares to fit a straight line.
2. Write the normal equation in method of least squares to fit a parabola.
3. Is it necessary that the curve due to method of least squares agree at the points in the given
data?
4. In method of least squares to fit a straight line, is the following statement is ture?
"The mean of the data points is always a point on the straight line."

Section-B
Subjective Questions

1. Derive the normal equations in fitting a straight line in Method of least squares.

2. Derive the normal equations in fitting a parabola in Method of least squares.

3. If p is the pull required to lift the weight by means of a ulley block, find a linear law of the
form p=a+bw, connecting p and w, using the following data:

w (lb): 50 70 100 120

p (lb): 12 15 21 25
Compute p, when w=150lb.

4. The results of measurement of electric resistance R of a copper bar at various temperatures

are listed below:
t: 19 25 30 36 40 45 50
R: 76 77 79 80 82 83 85
Find a relation R=a+bt when a and b are constants to be determined by you.

5. Write a code in MATLAB to implement Least Squares Method to fit a straight line.

6. Derive the normal equations in fitting a parabola in Method of least squares.

7. Write a code in MATLAB to implement Least Squares Method to fit a parabola curve.

8. Fit a parabola y= a+bx+cx2 to the following data:

x: 2 4 6 8 10
y: 3.07 12.85 31.47 57.38 91.29
9. The following table gives the results of the measurements of train resistances; V is the velocity
in miles per hour. R is the resistance in pounds per ton:

V: 20 40 60 80 100 120
R: 5.5 9.1 14.9 22.8 33.3 46.0
2
If R is related to V by the relation R=a+bV+cV , find a, b and c.

10. Fit the exponential curve to the following data:

x: 2 4 6 8
y: 25 38 56 84
11. Predict y at x=3.75 , by fitting a power curve to the given data:

x: 1 2 3 4 5 6
y: 298 4.26 5.21 6.10 6.80 7.50
12. Write a code in MATLAB to implement Least Squares Method to fit a power curve.