Transcript Slide 1

CEE203
NUMERICAL METHODS
IN ELECTRICAL
ENGINEERING
Assist. Prof. Dr. Çağatay ULUIŞIK
[email protected]
OUTLINE
 Introducing MatLab
 Fundamental Terms, Basic Definitions
 Solving Nonlinear Equations – Root Finding
 Solving Systems of Linear Equations
 Taylor Series
 Numerical Differentiation
 Numerical Integration
 Curve Fitting and Interpolation
 Ordinary Differential Equations
2
INTRODUCING MATLAB
Vectors, Arrays and Matrices
Assign a value to a variable
(observe the effect of the semicolon):
Creating a row vector of size (14)
Creating a column vector of size (41)
Creating an array of size (19) using the colon sign (:)
3
INTRODUCING MATLAB
Creating an array of size (19) using the for loop :
Creating a matrix of size (33) :
Creating a matrix of size (33) using 2 for loops :
4
BASIC MATRIX OPERATIONS IN MATLAB
5
GRAPHIC UTILITIES OF MATLAB
6
GRAPHIC UTILITIES OF MATLAB
7
COMPLEX NUMBERS IN MATLAB
Defining a Complex Number :
Conjugate :
Addition :
Magnitude :
Angle in Radians :
Subtraction :
Multiplication :
8
Angle in Degrees :
Division :
LOADING DATA FROM A FILE AND SAVING DATA TO A FILE
9
STATISTICS IN MATLAB
Generation of Uniformly Distributed Random Numbers :
Generation of Uniformly Distributed Random Integers :
Mean Value :
Variance :
Standard Deviation :
10
INTRODUCING MATLAB
Roots of a Polynom :
P( x )  x 3  2 x 2  3 x  4
Symbolic Integration :
Numeric Integration :
clc
: To clear the command window
clear all : To clear all the variables in memory
11
FORMAT OF DATA
round
floor
ceil
Round towards nearest integer
Round towards minus infinity
Round towards plus infinity
12
EXAMPLES OF MATLAB SCRIPTS
Example 1 : Calculating of the sum of integers from 1 to a user-specifed value n
a) using a formula b)using iterations
Example 2 : Creating and Calling a function
13
EXAMPLES OF MATLAB SCRIPTS
Example 3 : Error of the first n terms in the Taylor Series of the function ex
a) Calculating the error for a given n value
b) Calculating the minimum required number of terms n for a given error.
a)
b)
14
FUNDAMENTAL TERMS
 Accuracy: Closeness of agreement between a measured /
computed value and a true value.
 Measurement accuracy: the ability of an instrument to
measure the true value to within some stated error
specifications.
 Numerical accuracy: the degree to which the numerical
solution to the approximate physical problem
approximates the exact solution to the approximate
physical problem.
15
MEASUREMENT ACCURACY
 Comparison of a measured value with the real value
Accuracy =
True val. – Meas. val
True value
(%)
Value
interval
True value ?
Mean
16
Meas. Value
ERROR
 It is the difference between a measured or calculated
value of a quantity and its exact value.
 Systematic error plagues experiments or calculations
caused by negative factors. For example, a DC voltage
component, which unintentionally is present, e.g.,
because of a failure on the blockage capacitor, is a
systematic error. These can be removed once
understood/discovered via controls and calibration.
 Random error is an error which is always present, but
varies unpredictably in size and direction. They are related
to the scatter in the data obtained under fixed conditions
which determine the repeatability (precision) of the
measurement and follow well-behaved statistical rules.
17
ERROR
Absolute error:
It is an error that is expressed in physical units. It is the
absolute value of the difference between the measured
value and the true value (or the average value if the true
value is not known) of a quantity.
Relative error:
An error expressed as a fraction of the absolute error to
the true (or average) value of a quantity. It is always given
as a percentage.
18
ERROR
 A number is represented with a finite (fixed) number of
digits called word length.
 Precision of a measurement or the accuracy of a
computation bounds the number of significant digits. All
non-zero digits beyond the number of significant digits at the
right of a number are removed in one of two ways;
truncation or round-off.
 For example, 53.0534 has 6 significant digits. If it is going to
be represented only by 4 significant digits both truncation
and round-off processes yield 53.05. On the other hand, if
the number of significant digits will be 3, then truncation
and round-off processes yield, 53.0 and 53.1, respectively.
19
ERROR CALCULATIONS
 A measurement/calculation result should be given as
(a  a), which means the value may be anything between
(a - a) and (a + a).
 For example, if the measured speed is 98  3 km/hr, then
the real value may be anything between 95 and 101 km/hr.
 Total error is the sum of individual errors for the arithmetic
combination of two measurements.
 Total relative error is the sum of individual relative errors
for the mulitlicative combination of two measurements.
20
ERROR CALCULATIONS
 A parameter (say “A”) is going to be estimated / calculated
from two measurements (say “B” and “C”), with the
measurement errors of “b” and “c”, respectively.
 If A = B + C then total error will be
a bc
 If A = B  C then total error will be
b c
a  A   
B C
21
ERROR CALCULATIONS
 In general, the total (propagated) error is obtained from these
two properties as
m
C  A B
n

c
a
b
m
n
C
A
B
 For a multi-variable function, the total error is calculated from
y  f ( x1 , x2 , x3 ,...)
 f 
 f 
 f 
 x3  ...




y  
x1  
x2  


 x1 
 x2 
 x3 
22
UNCERTAINTY
 Uncertainty is a range that is likely to contain the true
value of a quantity being measured or calculated.
Uncertainty can be expressed in absolute or relative
terms.
 Modeling uncertainty is defined as the potential
deficiency due to a lack of information.
 Modeling error is the recognizable deficiency not due to
a lack of information but due to the approximations and
simplifications made there.
 Measurement error is the difference between the
measured and true values, while measurement
uncertainty is an estimate of the error in a measurement.
23
UNCERTAINTY
 Modeling and simulation uncertainties occur during the
phase of
 Conceptual modeling of the physical system





Mathematical modeling of the conceptual model
Discretization
Computer modeling of the mathematical model
Computer modeling of the discrete model
Numerical solution
 Numerical uncertainties occur during computations due to
round-off, turncation, non-convergence, artificial
dissipations, etc.
24
ERROR / UNCERTAINTY
Digital Voltmeter Accuracy:
(0.25 % Reading + 2 digit)
 Meaning; multiply the value you read on the multimeter by 0.25
% (i.e., by 0.0025) . This is called scaling error.
 Add 2 times the value of the least significant digit to the result.
This is called quantization error.
Example: You read 15.00 V from a 20 V scale.
• Measurement error : 0.0025  15.00 = 0.0375 V
• Quantization error : 2  0.01 = 0.02 V
:  0.0375 + 0.02 =  0.0575 V
• Total error
25
MODELING ERROR
 Truncation error is also defined for mathematical modeling.
For example, Taylor’s expansion or a Fourier-series
representation are used to replace a function in terms of an
infinite summation.
 Taking only a given number of low-order terms, and
neglecting the rest of the higher-order terms, introduces a
truncation error.
 For example, the Taylor expansion of the exponent function:
2
3
4
n
x
x
x
x
x
e x  1      ....   ...
1! 2! 3! 4!
n!
yields (for n=3)
Absolute error 


n4
26
xn
n!
FDTD MODELING ERROR
Derivative
f(x)
Mathematical definition

f ( x  x)  f ( x)
f ( x)  Lim
x
x
x  0
x
Finite difference approximation
2Pt

f ( x  x)  f ( x)
f ( x) 
x
x
+ O(∆x)
3Pt

f ( x  x)  f ( x  x)
f ( x) 
x
2x
+ O(∆x2)
27
SOLVING NONLINEAR EQUATIONS
 Fixed Point Iteration Method
 Bisection Method
 Secant Method
 Newton-Raphson Method
 MatLab Built-in Functions
28
FIXED POINT ITERATION METHOD
The equation f(x)=0 is rewritten in the form :
x  g( x )
The intersection point of the graphs of the functions y=x and y=g(x) is the
solution and is called the fixed point .
The numerical solution is obtained by an iterative process.
First a xi value near the fixed point is chosen and substituted into g(x).
The solution is substituted back into g(x). So the iteration formula is given by:
xi 1  g( xi )
29
FIXED POINT ITERATION METHOD
y=x
g(x2)
g(x4)
g(x6)
g(x7)
g(x5)
g(x3)
g(x)
g(x1)
x
x2=
x4=
x6= x7= x5=
x3=
g(x1)
g(x3) g(x5) g(x6) g(x4)
g(x2)
30
x1
FIXED POINT ITERATION METHOD
Example:
Convergence Test:
6
g( x ) 
x 1
x2  x  6  0
x( x  1 )  6
g ( x ) 
xi 1  g ( xi ) 
x  6 /( x  1 )
6
xi  1
5
y=x
3
2
g(x)=6/(x+1)
1
0
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
31
( x  1 )2
g ( 2 ) 
6
4
6
5.5
6
6
( 2  1 )2

2
1
3
i
xi
1
2
3
4
5
6
7
8
9
10
11
6.0000
0.8571
3.2308
1.4182
2.4812
1.7235
2.2030
1.8732
2.0882
1.9429
2.0388
FIXED POINT ITERATION METHOD
Choosing the appropriate iteration function g(x)
The fixed point iteration method converges if in the neighborhood of
the fixed point the derivative of g(x) has an absolute value smaller
than 1. (also called Lipschitz continious)
Example:
x e0.5x  1.2x  5  0
The plot of the function shows that the
15
equation has a solution between 1 and 2.
10
f(x)=xe0.5x+1.2x-5
The actual solution is at x=1.5050 .
5
0
-5
0
g ( x )  1
x
0.5
1
1.5
2
2.5
3
32
FIXED POINT ITERATION METHOD
Case A
Case B
5  x e 0.5 x
g( x )  x 
1.2
 e0.5 x ( 1  0.5 x )
1.2
g ( x ) 
 e 0.5 ( 1  0.5 )
g ( 1 ) 
 2.0609
1.2
g ( 2 ) 
 e( 1  1 )
 4.5305
1.2
Divergent
i
1
2
3
4
5
6
7
8
9
10
11
Case C
5
g( x )  x 
g( x )  x 
e 0.5 x  1.2
 5e0.5 x
g ( x ) 
2( e0.5 x  1.2 )2
 5e0.5
g ( 1 ) 
2( e
0.5
2
 1.2 )
 5e
g ( 2 ) 
2( e  1.2 )2
 0.5079
 0.4426
Convergent
xi
1.0000
2.7927
-5.2367
4.4849
-31.0262
4.1667
-23.7195
4.1668
-23.7223
4.1668
-23.7223
g ( x ) 
g ( 1 ) 
g ( 2 ) 
5  1.2 x
e0.5 x
 3.7  0.6 x
e0.5 x
 3.7  0.6
e0.5
 1.8802
 3.7  0.6  2
 0.9197
e
Divergent
i
xi
i
xi
1
2
3
4
5
6
7
8
9
10
11
1.0000
1.7552
1.3869
1.5622
1.4776
1.5182
1.4987
1.5080
1.5035
1.5057
1.5046
1
2
3
4
5
6
7
8
9
10
11
1.0000
2.3048
0.7057
2.9183
0.3482
3.8500
0.0554
4.7986
-0.0688
5.2606
-0.0946
33
BISECTION METHOD
STEP 1: Choose two points a and b such that a solution exists between them :
f(a) f(b) < 0
STEP 2: Find c=(a+b)/2
STEP 3: Determine whether the true solution is between a and c or between c and b.
Then select the subinterval that contains the true solution
If f(a) f(c) < 0 , solution is between a and c.
anew=a,
bnew=c
If f(b) f(c) < 0 , solution is between c and b.
anew=c,
bnew=b
STEP 4: Repeat steps 2 and 3 until (b-a) is less than a specified tolerance
34
BISECTION METHOD
80
80
f(b)
f(a)  f(b) < 0
60
60
c=(a+b) /2
40
40
20
0
True Root
a
c
b
0
a
c
b
f(a)
f(a)
-20
-20
5
f(b)
20
6
7
8
9
10
11
12
5
13
80
80
60
60
40
40
f(b)
20
6
7
8
9
10
11
12
13
8
9
10
11
12
13
20
f(b)
a
0
c
b
a
0
f(a)
f(a)
-20
b
-20
5
6
7
8
9
10
11
12
13
5
35
6
7
BISECTION METHOD
f ( x )  x 2  3x  28
a
b
c
f(c)
abs(b-a)
5.0000
12.0000
8.5000
18.7500
7.0000
5.0000
8.5000
6.7500
-2.6875
3.5000
6.7500
8.5000
7.6250
7.2656
1.7500
6.7500
7.6250
7.1875
2.0977
0.8750
6.7500
7.1875
6.9688
-0.3428
0.4375
6.9688
7.1875
7.0781
0.8655
0.2188
6.9688
7.0781
7.0234
0.2584
0.1094
6.9688
7.0234
6.9961
-0.0430
0.0547
6.9961
7.0234
7.0098
0.1075
0.0273
6.9961
7.0098
7.0029
0.0322
0.0137
6.9961
7.0029
6.9995
-0.0054
0.0068
6.9995
7.0029
7.0012
0.0134
0.0034
6.9995
7.0012
7.0004
36
0.0040
0.0017
NEWTON RAPSON METHOD
f(x1)
f(x)
Slope: f’(x1)
f(x2)
Exact Root
Slope:
f’(x4)
x5
x4
Slope:
f’(x3)
xi 1  xi 
37
Slope:
f’(x2)
x3
f ( xi )
f ( xi )
x
x2
x1
NEWTON RAPSON METHOD
f ( x )  0.5  2 x  1
xi
f(xi)
6.00000000
31.00000000
4.60238918
11.14583002
1.3976
3.27847524
3.85164919
1.3239
2.13314198
1.19335902
1.1453
1.34820293
0.27297398
0.7849
1.03883431
0.02728345
0.3094
1.00051801
0.00035912
0.0383
1.00000009
0.00000006
0.0005
f ( xi )
xi 1  xi 
f ( xi )
38
xi -xi-1
SECANT METHOD
f(x1)
f(x)
f(x2)
f(x3)
Exact Root
f(x5)
f(x4)
x
x6
xi 1  xi 
x5
x4
x3
xi  xi 1
f ( xi )
f ( xi )  f ( xi 1 )
39
x2
x1
SECANT METHOD
f ( x )  0.5  2 x  1
xi
f(xi)
6.00000000
31.00000000
5.00000000
15.00000000
1.0000
4.06250000
7.35419026
0.9375
3.16075727
3.47149501
0.9017
2.35451439
1.55711029
0.8062
1.69873759
0.62308391
0.6558
1.26127246
0.19853535
0.4375
1.05669682
0.04008167
0.2046
1.00494836
0.00343583
0.0517
1.00009654
0.00006692
0.0049
1.00000017
0.00000011
0.0001
xi 1  xi 
xi -xi-1
xi  xi 1
f ( xi )
f ( xi )  f ( xi 1 )
40
MATLAB BUILT-IN FUNCTIONS
The roots(p) command finds the roots of a polynomial, whose coefficients
are defined in the row vector p. The following MatLab script find the roots of
the polynom
P(x)=x3-3x2-13x+15
>> coef=[1 -3 -13 15];
>> xi=roots(coef)
xi =
5.0000
-3.0000
1.0000
The fzero(‘function’,x0) command finds the roots of the ‘function’ near x0.
The following MatLab script find the roots of the function f(x)=x3-3x2-13x+15
>> x=fzero('x^3-3*x^2-13*x+15',7)
x=5
>> x=fzero('x^3-3*x^2-13*x+15',3)
x=1
41
SYSTEMS OF LINEAR EQUATIONS
 Gauss Elimination Method
 Gauss-Jordan Elimination Method
 LU Decomposition Method
42
GAUSS ELIMINATION METHOD
2 x1  x2  x3  4
4 x1  2 x2  5 x3  22
5 x1  2 x2  x3  15
 2  1 1  4
4 2 - 5 22 


5 2  1 15 
1  0.5 0.5  2
4
2
- 5 22 

5
2
 1 15 
1  0.5 0.5  2
0

1

1.75
7.5


0 4.5  3.5 25 
1  0.5 0.5  2
0
4
- 7 30 

5
2
 1 15 
2 
1  0.5 0.5
0

1

1.75
7.5


0
0
4.375  8.75
1  0.5 0.5  2
0

1

1.75
7.5


0
0
1
 2
Backward Substitution :
x3  2
x1  1
x2  1.75 x3  7.5
x2  1.75  ( 2 )  7.5
x1  0.5  4  0.5( 2 )  2

x2  4
x2  4

x3  2
x1  1
43
1  0.5 0.5  2
0
4
 7 30 

0 4.5  3.5 25 
LU DECOMPOSITION METHOD
LU x  b
Ux y
Ly  b
A x  b 

A  LU 
• A is the coefficients matrix.
• U is an upper triangular matrix and is obtained at the end of the Gauss
elimination procedure.
• L is a lower triangular matrix which has all 1’s on the diagonal and the
elements below the diagonal are the multipliers mij that multiply the pivot
equation when it is used to eliminate the elements below the pivot
coefficient at the Gauss elimination procedure.
• L is a lower triangular matrix and y is obtained from Ly=b using forward
substitution.
• U is an upper triangular matrix and x is obtained from
backward substitution.
44
Ux=y
using
LU DECOMPOSITION METHOD
2 x1  x2  x3  4
2  1 1 
4 2 - 5 
4 x1  2 x2  5 x3  22 
 m21 =2 
5 2  1



5 x1  2 x2  x3  15
A
1 
2  1
U  0 4
 7 
0 0 4.375
0
 1
L  m 21
1
 m31 m32
2  1 1 
0 4 - 7  m =2.5 

 31
5 2  1
0  1
0
0
0   2
1
0
1 2.5 1.125 1
y
A  LU
1   x1    4 
2  1
0 4
  x    30 

7

 2  

0 0 4.375  x3   8.75
  


0
0  y1   4
1
2
  y    22 
1
0

 2   
2.5 1.125 1  y3   15 
  

L
1 
2  1
0 4
 m =1.125

7

 32
0 4.5  3.5
b
U
x
y
Forward Substitution :
Backward Substitution :
y1  4
4.375x3  8.75 
2 y1  y2  22
4 x2  7 x3  30
2  ( 4 )  y2  22 
4 x2  7  (-2)  30 
y2  30
x1  1
x2  4
2 x1  x2  x3  4
2.5 y1  1.125y2  y3  15
2.5  ( 4 )  1.125 30  y3  15 
x3  2
y3  8.75
45
2 x1  4  2  4 
x1  1
x2  4
x3  2
GAUSS-JORDAN ELIMINATION METHOD
2 x1  x2  x3  4
4 x1  2 x2  5 x3  22
5 x1  2 x2  x3  15
 2  1 1  4
4 2 - 5 22 


5 2  1 15 
1  0.5 0.5  2
4

2
5
22


5
2
 1 15 
1  0.5 0.5  2
0

4
7
30


5
2
 1 15 
1  0.5 0.5  2
0

1

1.75
7.5


0 4.5  3.5 25 
2 
1  0.5 0.5
0

1

1.75
7.5


0
0
4.375  8.75
1  0.5 0.5  2
0

1

1.75
7.5


0
0
1
 2
1  0.5 0  1
0

1
0
4


0
0
1  2
1 0 0 1 
0 1 0 4 


0 0 1  2
x1  1
x2  4
x3  2
46
1  0.5 0.5  2
0

4

7
30


0 4.5  3.5 25 
1  0.5 0.5  2
0

1
0
4


0
0
1  2
GAUSS-JORDAN ELIMINATION METHOD
x1  x2  2 x3  3 x4  16
 1 - 1 2 - 3 16 
 2 6 4 -1 3 


 3 4 - 7 5 - 39


- 2 - 3 4 - 1 22 
2 x1  6 x2  4 x3  x4  3
3 x1  4 x2  7 x3  5 x4  39
 2 x1  3x2  4 x3  x4  22
1 - 1 2 - 3 16 
0 8
0
5 - 29

0 7 - 13 14 - 87 


0 - 5 8 - 7 54 
-3
16 
1 - 1 2
0 1
0 0.625 - 3.625

0 7 - 13 14
- 87 


-7
54 
0 - 5 8
1 - 1
0 1

0 0

0 0
2
-3
16 
0
0.625
- 3.625
1 - 0.74038 4.7404

8 - 3.875 35.875
1 - 1
0 1

0 0

0 0
2 - 3 16
0 0 - 3
1 0 4

0 1 - 1
1 - 1
0 1

0 0

0 0
1 - 1
0 1

0 0

0 0
-3
16 
1 - 1 2
0 1
0 0.625 - 3.625 

0 0 - 13 9.625 - 61.625


-7
54 
0 - 5 8
2
-3
16 
0
0.625
- 3.625 
1 - 0.74038 4.7404 

0 2.0481 - 2.0481
2
0
1
0
0 13
0 - 3
0 4

1 - 1
 1 - 1 2 - 3 16 
 0 8 0 5 - 29


 3 4 - 7 5 - 39


- 2 - 3 4 - 1 22 
1 - 1
0 1

0 0

0 0
1 - 1
0 1

0 0

0 0
47
0
0
1
0
 1 - 1 2 - 3 16 
0 8
0
5 - 29

 0 7 - 13 14 - 87


- 2 - 3 4 - 1 22 
-3
16 
1 - 1 2
0 1
0
0.625 - 3.625 

0 0 - 13 9.625 - 61.625


8 - 3.875 35.875 
0 0
2
-3
16 
0
0.625
- 3.625
1 - 0.74038 4.7404

0
1
-1 
0 5
0 - 3
0 4

1 - 1
1
0

0

0
0
1
0
0
1 - 1
0 1

0 0

0 0
0
0
1
0
2
-3
16 
0 0.625 - 3.625
1
0
4 

0
1
-1 
0 2
0 - 3
0 4

1 - 1
x1  2
x2  -3
x3  4
x4  -1
GAUSS-JORDAN ELIMINATION METHOD
MATLAB SCRIPT
A=[2 -1 1 -4; 4 2 -5 22; 5 2 -1 15];
[n,col]=size(A);
%------------------------------------------------------------------for i=1:n
A(i,:)=A(i,:)/A(i,i);
disp('A ='); disp(A); pause;
for k=i+1:n
dd=A(k,i);
A(k,:)=A(k,:)-dd*A(i,:); A
pause;
end
end
%------------------------------------------------------------------for i=n:-1:2
for k=i-1:-1:1
dd=A(k,i);
A(k,:)=A(k,:)-dd*A(i,:); A
pause;
end
end
x=A(:,col);
disp(' The unknowns are :'); disp(x);
rref(A) % MatLab built-in function
48
TAYLOR SERIES
Taylor Series :
f ( x0 )
f ( x0 )
f ( 4 )( x0 )
2
3
f ( x )  f ( x0 )  f ( x0 )( x  x0 ) 
( x  x0 ) 
( x  x0 ) 
( x  x0 )4  
2!
3!
4!
Maclaurin Series :
f ( 0 ) 2 f ( 0 ) 3 f ( 4 )( 0 ) 4
f ( x )  f ( 0 )  f ( 0 )x 
x 
x 
x  
2!
3!
4!
Example :
f ( x )  cos( x )
f (0 ) 1
f ( x )   sin( x )
f ( 0 )  0
f ( x )   cos( x )
f ( 0 )  1
f ( x )  sin( x )
f ( 0 )  0
f ( 4 ) ( x )  cos( x )
f ( 4 )( 0 )  1
f ( 5 ) ( x )   sin( x )
f ( 5 )( 0 )  0
f ( 6 ) ( x )   cos( x )
f ( 6 ) ( 0 )  1
f ( 7 ) ( x )  sin( x )
f ( 7 )( 0 )  0
f ( 8 ) ( x )  cos( x )
f ( 8 )( x )  1
cos( x )  1 
1 2 1 4 1 6 1 8 1 10 1 12
x  x  x  x 
x 
x  
2
4!
6!
8!
10!
12!
49
TAYLOR SERIES
2
y=cos(x)
2 Terms
3 Terms
4 Terms
5 Terms
6 Terms
7 Terms
8 Terms
1.5
1
0.5
0
-0.5
-1
-1.5
-2
-8
-6
-4
-2
0
50
2
4
6
8
NUMERICAL DIFFERENTIATION
 Two-Point Forward Difference Method
 Two-Point Backward Difference Method
 Two-Point Central Difference Method
 Three-Point Forward Difference Method
 Three-Point Backward Difference Method
51
NUMERICAL DIFFERENTIATION
Forward Difference
Approximated
derivative
Backward
Difference
f(xi+1)
f(xi)
f(x)
f(xi)
Approximated
derivative
True
derivative
f(xi-1)
f(x)
h
xi
h
xi+1
xi-1
Central Difference
f(x)
f(xi+1)
xi
Forward Difference:
f ( xi 1 )  f ( xi )
df

dx x x
xi 1  xi
True derivative
i
Backward Difference:
f ( xi )  f ( xi 1 )
df

dx x x
xi  xi 1
Approximated
derivative
i
f(xi-1)
Central Difference:
f ( xi 1 )  f ( xi 1 )
df

dx x x
xi 1  xi 1
2h
xi-1
x i +1
52
True
derivative
i
NUMERICAL DIFFERENTIATION
xi 1  xi  xi  xi 1  h
Truncation
Error:
2 Point Forward Difference:
f ( xi 1 )  f ( xi )
df

dx x x
h
i
O( h )
2 Point Backward Difference:
f ( xi )  f ( xi 1 )
df

dx x x
h
i
O( h )
2 Point Central Difference:
f ( xi 1 )  f ( xi 1 )
df

dx x x
2h
i
O( h 2 )
3 Point Forward Difference:
 3 f ( xi )  4 f ( xi 1 )  f ( xi 2 )
df

dx x x
2h
i
O( h 2 )
3 Point Backward Difference:
f ( xi 2 )  4 f ( xi 1 )  3 f ( xi )
df

dx x x
2h
i
O( h 2 )
53
NUMERICAL DIFFERENTIATION
Example : For f(x)=ln(x) estimate the value of f’(0.2) taking h=0.01.
xi-2 0.18
f(xi-2) ln(0.18)
xi-1 0.19
f(xi-1) ln(0.19)
f ( x ) 
xi
f(xi)
f ( 0.2 )  5
0.2
ln(0.2)
xi+1 0.21
f(xi+1) ln(0.21)
xi+2 0.22
f(xi+2) ln(0.22)
1
x
Method Used
f’(0.2)
Abs.
Error
4.879
0.121
f ( 0.2 ) 
f ( xi1 )  f ( xi ) ln( 0.21)  ln( 0.2 )

 4.879
h
0.01
2 Point Forward
Difference
f ( 0.2 ) 
f ( xi )  f ( xi 1 ) ln( 0.2 )  ln( 0.19 )

 5.1293
h
0.01
2 Point Backward
5.1293
Difference
f ( 0.2 ) 
f ( xi1 )  f ( xi1 ) ln( 0.21)  ln( 0.19 )

 5.0042
2h
0.02
2 Point Central
Difference
f ( 0.2 ) 
f ( 0.2 ) 
3 Point Forward
 3 f ( xi )  4 f ( xi1 )  f ( xi2 )  3ln( 0.2 )  4 ln( 0.21 )  ln( 0.22 )

 4.9925
Difference
2h
0.02
3 Point Backward
f ( xi2 )  4 f ( xi1 )  3 f ( xi ) ln( 0.18 )  4 ln( 0.19 )  3ln( 0.2 )

 4.9906
Difference
2h
0.02
54
0.1293
5.0042 0.0042
4.9925 0.0075
4.9906 0.0094
NUMERICAL DIFFERENTIATION
h
f’(0.8)
Abs. Error
0.7500
0.5000
0.4000
0.3000
0.2000
0.1000
0.0500
0.0100
0.0050
0.0010
0.0005
0.0001
0.8819
0.9710
1.0137
1.0615
1.1157
1.1778
1.2125
1.2423
1.2461
1.2492
1.2496
1.2499
0.3681
0.2790
0.2363
0.1885
0.1343
0.0722
0.0375
0.0077
0.0039
0.0008
0.0004
0.0001
h
f’(0.8)
Abs. Error
0.7500
0.5000
0.4000
0.3000
0.2000
0.1000
0.0500
0.0100
0.0050
0.0010
0.0005
0.0001
2.2893
1.4663
1.3733
1.3141
1.2771
1.2566
1.2516
1.2501
1.2500
1.2500
1.2500
1.2500
1.0393
0.2163
0.1233
0.0641
0.0271
0.0066
0.0016
6.51e-005
1.63e-005
6.51e-007
1.63e-007
6.51e-009
2 Point
Forward
Difference:
2 Point
Backward
Difference:
f ( x )  ln( x )
f ( x )  1 / x
f ( 0.8 )  1.25
2 Point
Central
Difference:
3 Point
Forward
Difference:
55
h
f’(0.8)
0.7500
0.5000
0.4000
0.3000
0.2000
0.1000
0.0500
0.0100
0.0050
0.0010
0.0005
0.0001
3.6968
1.9617
1.7329
1.5667
1.4384
1.3353
1.2908
1.2579
1.2539
1.2508
1.2504
1.2501
Abs. Error
2.4468
0.7117
0.4829
0.3167
0.1884
0.0853
0.0408
0.0079
0.0039
0.0008
0.0004
0.0001
h
f’(0.8)
Abs. Error
0.7500
0.5000
0.4000
0.3000
0.2000
0.1000
0.0500
0.0100
0.0050
0.0010
0.0005
0.0001
1.0597
1.1311
1.1609
1.1903
1.2178
1.2399
1.2472
1.2499
1.2499679
1.2499987
1.24999967
1.249999987
0.1903
0.1189
0.0891
0.0597
0.0322
0.0101
0.0028
0.0001
3.21e-005
1.298e-006
3.25e-007
1.30e-008
NUMERICAL DIFFERENTIATION
3
2.8
2.6
2.4
2.2
2
1.8
1.6
f ( x )  x2  4x  2
1.4
f ( x )  2 x  4
1.2
1
Slope  f ( 1.5 )  1
1.3
1.4
1.5
1.6
1.7
1.8
1.9
2
56
2.1
2.2
2.3
2.4
2.5
2.6
2.7
NUMERICAL DIFFERENTIATION
f ( x )  x2  4x  2
f ( x )  2 x  4
f ( 1.5 )  1
h
f’(1.5)
Abs.
Error
1.0000
0
1.0000
0.8000
0.2000
0.8000
0.5000
0.5000
0.5000
0.2000
0.8000
0.2000
0.1000
0.9000
0.1000
0.0500
0.9500
0.0500
0.0100
0.9900
0.0100
0.0050
0.9950
0.0050
0.0010
0.9990
0.0010
57
NUMERICAL INTEGRATION
 Rectangle Methods
 Midpoint Method
 Trapezoidal Method
 Simpsons Methods
Simpsons 1/3 Method
Simpsons 3/8 Method
58
NUMERICAL INTEGRATION
1. Rectangle Method
f(b)
f(b)
f(a)
f(a)
a
a
b
b
a
a f ( x )dx  ( b  a ) f ( b )
Trapezoidal Method
f(x)
f(b)
f(b)
f(a)
f(a)
a
b
a
a
b
f ( x )dx  ( b  a ) f (
ab
)
2
b
b
f ( x )dx  ( b  a ) f ( a )
Midpoint Method
f(x)
2. Rectangle Method
f(x)
b
59
a
f ( x )dx  ( b  a )
f(x)
b
 f ( a )  f ( b )
2
NUMERICAL INTEGRATION
1. Rectangle Method
f(x)
f(x2)
f(x1)
I1
I2
I3
Ii
Ii+1
IN-1
IN
h
h
h
h
h
h
h
x1
a
ba
h
N
x2
x3

xi
b
N
a
i 1
xi+1 
N
 f ( x )dx   I i  h f ( xi )
60
i 1
xN-1
xN
xN+1
b
NUMERICAL INTEGRATION
2. Rectangle Method
f(x)
f(x2)
f(x1)
x1
a
I1
I2
I3
Ii
Ii+1
IN-1
IN
h
h
h
h
h
h
h
x2
x3

xi
xi+1 
N 1
b
N
a
i 1
 f ( x )dx   I i  h  f ( xi )
61
i 2
xN-1
xN
xN+1
b
NUMERICAL INTEGRATION
Midpoint Method
f(x)
 x  x3 
f 2

 2 
x x 
f 1 2
 2 
I1
I2
I3
Ii
Ii+1
IN-1
IN
h
h
h
h
h
h
h
x1
x2
x3
(x +x )
)
a (x +x
2
2
1
2
2
b

a

xi+1 
xi
xN-1
3
 xi  xi 1 
f ( x)dx   I i  h f 

 2 
i 1
i 1
N
N
62
xN
xN+1
b
NUMERICAL INTEGRATION
Trapezoidal Method
f(x)
f(x2)
f(x1)
x1
a
I1
I2
I3
Ii
Ii+1
IN-1
IN
h
h
h
h
h
h
h
x2
b
x3

xi+1 
xN-1
h N
f ( x )dx 
Ii 
f ( xi )  f ( xi 1 )
2 i 1
i

1
a
63

N
xi


xN
xN+1
b
NUMERICAL INTEGRATION
Simpson 1/3 Method
f(x)
A second order polynomial g(x) is used
to approximate the integrand f(x)
g(x)
f(x1)
h/2=k
h/2=k
x1
a
x2
b

a
k
ab

f ( x)dx   f (a)  4 f (
)  f (b)
3
2

64
x3
b
NUMERICAL INTEGRATION
Simpson 3/8 Method
f(x)
A third order polynomial g(x) is used
to approximate the integrand f(x)
g(x)
f(x1)
h/3=k
x1
a
b

a
h/3=k
x2
h/3=k
x3
x4
b

3k 
f ( x)dx 
f (a)  3 f ( x2 )  3 f ( x3 )  f (b) 

8 

65
NUMERICAL INTEGRATION
200
h  2.5
200
150
150
100
100
50
50
0
0
200
1
2
3
4
5
6
7
8
9
10 11 12
h 1
0
0
200
150
150
100
100
50
50
0
0
1
2
3
4
5
6
7
8
9
10 11 12
0
0
h2
1
2
3
4
5
6
7
8
9
10 11 12
4
5
6
7
8
9
10 11 12
h  0.5
1
2
3
f ( x )  x3 -6618x2  96x
NUMERICAL INTEGRATION
1. Rectangular Method
Num
Int.
5
1115
2.5
1277.5
2
1310
1
1375
0.5
1407.5
0.1
1433.5
0.05
1436.8
0.01
1439.3
0.005 1439.7
0.001 1439.9
0.0005 1440
0.0001 1440
h
Abs.
Error
325
162.5
130
65
32.5
6.5
3.25
0.65
0.325
0.065
0.0325
0.0065
2. Rectangular Method
Rel.
Error %
22.569
11.285
9.0278
4.5139
2.2569
0.45139
0.22569
0.045139
0.022569
0.0045139
0.0022569
0.00045139
Num
Int.
5
1765
2.5
1602.5
2
1570
1
1505
0.5
1472.5
0.1
1446.5
0.05
1443.3
0.01
1440.6
0.005 1440.3
0.001 1440.1
0.0005 1440
0.0001 1440
h
11 3
x - 18x 2
1

 96x dx  1440
67
Abs.
Error
325
162.5
130
65
32.5
6.5
3.25
0.65
0.325
0.065
0.0325
0.0065
Rel. Error %
22.569
11.285
9.0278
4.5139
2.2569
0.45139
0.22569
0.045139
0.022569
0.0045139
0.0022569
0.00045139
NUMERICAL INTEGRATION
h
8
4
2
1
0.5
0.1
0.05
0.01
0.005
0.001
0.0005
0.0001
1. Rect. Method
Num
Abs.
Int.
Error
400
352
360
372
380
387.36
388.34
389.13
389.23
389.31
389.32
389.33
10.667
37.333
29.333
17.333
9.3333
1.9733
0.99333
0.19973
0.099933
0.019997
0.009999
0.002
2. Rect. Method
Num
Abs.
Int.
Error
720
512
440
412
400
391.36
390.34
389.53
389.43
389.35
389.34
389.34
330.67
122.67
50.667
22.667
10.667
2.0267
1.0067
0.20027
0.10007
0.020003
0.010001
0.002
10 3
x - 16x2
2

Midpoint Method
h
Num
Int.
8
4
2
1
0.5
0.1
0.05
0.01
0.005
0.001
0.0005
0.0001
304
368
384
388
389
389.32
389.33
389.33
389.33
389.33
389.33
389.33
Abs.
Error
85.333
21.333
5.3333
1.3333
0.33333
0.013333
0.0033333
0.00013333
3.3333e-005
1.3333e-006
3.3333e-007
1.3332e-008
Trapezoidal Method
h
Num
Int.
Abs.
Error
8
4
2
1
0.5
0.1
0.05
0.01
0.005
0.001
0.0005
0.0001
560
432
400
392
390
389.36
389.34
389.33
389.33
389.33
389.33
389.33
170.67
42.667
10.667
2.6667
0.66667
0.026667
0.0066667
0.00026667
6.6667e-005
2.6667e-006
6.6667e-007
2.6663e-008
1168
 73x - 40 dx 
 389.33
3
68
NUMERICAL INTEGRATION
Numerical
Absolute
Integration
Error
1. Rectangle Method
387.36
1.9733
2. Rectangle Method
391.36
2.0267
Midpoint Method
389.32
0.013333
Trapezoidal Method
389.36
0.026667
Simpsons 1/3 Method
389.33
5.6843e-014
Simpsons 3/8 Method
389.33
5.6843e-014
(h=0.1)
10 3
x - 16x2
2

1168
 73x - 40 dx 
 389.33
3
69
CURVE FITTING WITH A LINEAR FUNCTION
Curve Fitting with a Linear Function
(xi ,yi)
yi
yN
y1
f(x)=a1x+a0
ri
rN
(xN ,yN)
(x1 ,y1)
r1
y2
r2
(x2 ,y2)
x1
x2
xi
xN
A linear function f(x)=a1x+a0 is used to best fit the given data points (xi ,yi).
A residual ri at point (xi , yi) is the difference between the value yi of the data point
and the value of the function f(xi) used to approximate the data points : ri= yi -f(xi)
70
CURVE FITTING WITH A LINEAR FUNCTION
The overall error E is defined as the sum of the squares of the residuals :
E
N

i 1
ri2 
N
2




y

a
x

a
 i 1i 0
i 1
Linear least-squares regression :
The coefficients a1 and a0 of the linear function f(x)=a1x+a0 are determined such that the
error E has the smallest possible value. E is minimum if:
E
0
a0
E
0
a1
The coefficients a1 and a0 are found as:
a1 
Sx 
nS xy  S x S y
a0 
nS xx  S x 
2
N

i 1
xi
Sy 
S xxS y  S xyS x
nS xx  S x 2
N

i 1
S xy 
yi
71
N

i 1
xi yi
S xx 
N

i 1
xi2
CURVE FITTING WITH A LINEAR FUNCTION
Example : Use linear least-squares regression to determine the coefficients a1 and a0
in the function f(x)=a1x+a0 that best fits the data given below.
Sx 
Sy 
1
2
3
y
6
9.2
13
4
5
6
7
8
9
10
14.7 19.7 21.8 22.8 29.1 30.2 32.2
N
 xi  1  2  3  4  5  6  7  8  9  10  55
i 1
N
 yi  6  9.2  13  14.7  19.7  21.8  22.8  29.1 30.2  32.2  198.7
i 1
S xy 
S xx 
a1 
x
N
 xi yi  1 6  2  9.2  3 13  4 14.7  5 19.7  6  21.8  7  22.8  8  29.1 9  30.2  10  32.2  1337.7
i 1
N
 xi2  12  22  32  42  52  62  72  82  92  102  385
i 1
nS xy  S x S y
nS xx  S x 2

10 1337.7  55 198.7
10  385  552
 2.9679 a0 
72
S xxS y  S xyS x
nS xx  S x 2

385 198.7  1337.7  55
10  385  552
 3.5467
CURVE FITTING WITH A LINEAR FUNCTION
yi
40
35
f(xi)
(ri)2
ri= yi -f(xi)
6.0000
6.5145
-0.5145
0.2648
9.2000
9.4824
-0.2824
0.0798
13.0000
12.4503
0.5497
0.3022
14.7000
15.4182
-0.7182
0.5158
19.7000
18.3861
1.3139
1.7264
21.8000
21.3539
0.4461
0.1990
22.8000
24.3218
-1.5218
2.3159
29.1000
27.2897
1.8103
3.2772
30.2000
30.2576
-0.0576
0.0033
32.2000
33.2255
-1.0255
30
f(x)=2.9679x+3.5467
25
20
15
10
5
0
0
1
2
3
4
5
6
7
8
9
10
11
E=
73
+
1.0516
9.736
CURVE FITTING WITH A LINEAR FUNCTION
Example : Use linear least-squares regression to determine the coefficients a1 and a0
in the function f(x)=a1x+a0 that best fits the data given below.
x
2
5
6
8
9
13
15
y
7
8
10
11
12
14
15
Sx 
Sy 
N
 xi  2  5  6  8  9  13  15  58
i 1
N
 yi  7  8  10  11  12  14  15  77
i 1
S xy 
S xx 
a1 
N
 xi yi  2  7  5  8  6 10  8 11  9 12  13 14  15 15  717
i 1
N
 xi2  22  52  62  82  92  132  152  604
i 1
nS xy  S x S y
nS xx  S x 2

7  717  58  77
7  604  582
 0.64
a0 
74
S xxS y  S xyS x
nS xx  S x 2

604  77  717  58
7  604  582
 5.6968
CURVE FITTING WITH A LINEAR FUNCTION
16
14
f(x)=0.64x+5.6968
12
10
8
6
4
0
1
2
3
4
yi
7.0000
8.0000
10.0000
11.0000
12.0000
14.0000
15.0000
5
6
7
f(xi)
6.9769
8.8970
9.5370
10.8171
11.4572
14.0174
15.2975
8
9
10
11
12
13
14
ri= yi -f(xi)
(ri)2
0.0231
-0.8970
0.4630
0.1829
0.5428
-0.0174
-0.2975
0.0005
0.8046
0.2143
0.0334
0.2947
0.0003
0.0885
+
E=
75
1.4363
15
16
CURVE FITTING WITH A LINEAR FUNCTION
Example : Use linear least-squares regression to determine the coefficients a1 and a0 in the
function f(x)=a1x+a0 that best fits the data given below. Use f(x) to determine the
interpolated value n for x=4 and the extrapolated value m for x= -8.
Sx 
Sy 
-8
-7
-5
-1
0
2
4
5
6
y
m
15
12
5
2
0
n
-5
-9
N
 xi  -7 - 5 - 1  0  2  5  6  0
i 1
N
 yi  15  12  5  2  0 - 5 - 9  20
i 1
S xy 
S xx 
a1 
x
N
 xi yi  (-7) 15  (-5) 12  (-1)  5  0  2  2  0  5  (-5)  6  (-9)  -249
i 1
N
 xi2  (-7) 2  (-5)2  (-1) 2  02  22  52  62  140
i 1
nS xy  S x S y
nS xx  S x 2

7  (-249)  0  20
7 140  02
 - 1.7786
a0 
76
S xxS y  S xyS x
nS xx  S x 2

140  20  (-249)  0
7 140  02
 2.8571
CURVE FITTING WITH A LINEAR FUNCTION
20
m
15
f(x)=1.7786x+ 2.8571
10
5
0
n
-5
-10
-10
-8
-6
-4
0
-2
2
4
6
8
Interpolation : Using the data points for estimating the expexted values between the known points.
n  f ( 4 )  1.7786  (4)  2.8571 4.2571
Extrapolation : Using the data points for estimating the expexted values beyond the known points.
m  f ( 8 )  1.7786  (8)  2.8571 17.0857
77
LAGRANGE POLYNOMIALS
For n+1 data points (x0 , y0), (x1 , y1), (x2 , y2), ….. (xk , yk), ….. (xn , yn) a unique polynomial
of order n can be found, which passes through these n data points. A Lagrange interpolating
polynomial Pn(x) that passes through n points is defined as :
Pn ( x ) 
Lk ( x ) 
n
 yk Lk ( x )  y0 L0 ( x )  y1L1( x )   yk Lk ( x )    yn Ln ( x )
k 0
( x  x0 )( x  x1 ).....(x  xk 1 )( x  xk 1 ).....(x  xn )
( xk  x0 )( xk  x1 ).....(xk  xk 1 )( xk  xk 1 ).....(xk  xn )
1 , i  k
Lk ( xi )  
0 i  k
Pn ( xi )  yi
78
LAGRANGE POLYNOMIALS
Example : Find the Lagrange interpolation polynomial that passes through the points :
x
y
2
2.5
4
0.5 0.4 0.25
L0 ( x ) 
( x  x1 )( x  x2 )
( x  2.5 )( x  4 )

 ( x  2.5 )( x  4 )
( x0  x1 )( x0  x2 ) ( 2  2.5 )( 2  4 )
x0  2,
L1( x ) 
( x  x0 )( x  x2 )
( x  2 )( x  4 )
( x  2 )( x  4 )


( x1  x0 )( x1  x2 ) ( 2.5  2 )( 2.5  4 )
3/ 4
x2  4 ,
L2 ( x ) 
( x  x0 )( x  x1 )
( x  2 )( x  2.5 ) ( x  2 )( x  2.5 )


( x2  x0 )( x2  x1 ) ( 4  2 )( 4  2.5 )
3
P2 ( x ) 
y0  0.5
x1  2.5, y1  0.4
y2  0.25
2
 yk Lk ( x )  y0 L0 ( x )  y1L1( x )  y2 L2 ( x )
k 0
( x  2 )( x  4 )
( x  2 )( x  2.5 )
 0.5( x  2.5 )( x  4 )  0.4
 0.25
3/ 4
3
 0.5( x 2  6.5 x  10 ) 
1.6 2
1
( x  6 x  8 )  ( x 2  4.5 x  5 )
3
12
 0.05 x 2  0.425x  1.15
79
MATLAB SCRIPT
x=[2 2.5 4];
y=[0.5 0.4 0.25];
polyfit(x,y,2)
LAGRANGE POLYNOMIALS
x
2
2.5
4
The Lagrange interpolation polynomial that passes through the points :
y
0.5 0.4 0.25
P2( x )  0.05x2  0.425x  1.15
1.2
1
0.8
0.6
0.4
0.2
0
0
1
2
3
80
4
5
6
LAGRANGE POLYNOMIALS
Example : Find the Lagrange interpolation polynomial that passes through the points :
x
0
1
2
y
-1
2
7
L0 ( x ) 
( x  x1 )( x  x2 )
( x  1 )( x  2 )

 0.5( x  1 )( x  2 )
( x0  x1 )( x0  x2 ) ( 0  1 )( 0  2 )
x0  0, y0  1
x1  1,
y1  2
L1( x ) 
( x  x0 )( x  x2 ) ( x  0 )( x  2 )

  x( x  2 )
( x1  x0 )( x1  x2 ) ( 1  0 )( 1  2 )
x2  2 ,
y2  7
L2 ( x ) 
( x  x0 )( x  x1 )
( x  0 )( x  1 )

 0.5 x( x  1 )
( x2  x0 )( x2  x1 ) ( 2  0 )( 2  1 )
P2 ( x ) 
2
 yk Lk ( x )  y0 L0 ( x )  y1L1( x )  y2 L2 ( x )
k 0
 1  0.5( x  1 )( x  2 )  2 x( x  2 )  7  0.5 x( x  1 )
 0.5( x 2  3x  2 )  2 x 2  4 x  3.5 x 2  3.5 x
 x2  2x 1
81
LAGRANGE POLYNOMIALS
The Lagrange interpolation polynomial that passes through the points :
x
0
1
2
y
-1
2
7
P2 ( x )  x 2  2x  1
25
20
15
10
5
0
-5
-2
-1
0
1
82
2
3
4
LAGRANGE POLYNOMIALS
Example : Find the Lagrange interpolation polynomial that passes through the points :
x
0
1
3
4
y
3
0
0
15
L1( x ) 
L2 ( x ) 
L0 ( x ) 
( x  x1 )( x  x2 )( x  x3 )
( x  1 )( x  3 )( x  4 )  1


( x  1 )( x  3 )( x  4 )
( x0  x1 )( x0  x2 )( x0  x3 ) ( 0  1 )( 0  3 )( 0  4 ) 12
L3 ( x ) 
( x  x0 )( x  x1 )( x  x2 )
( x  0 )( x  1 )( x  3 ) 1

 x( x  1 )( x  3 )
( x3  x0 )( x3  x1 )( x3  x2 ) ( 4  0 )( 4  1 )( 4  3 ) 12
P3 ( x ) 
3
y1 L1( x )  y2 L2 ( x )  y3 L3 ( x )
 yk Lk ( x )  y0 L0 ( x )  

k 0
0
0
1
5

( x  1 )( x  3 )( x  4 )  x( x  1 )( x  3 )
4
4

1
1
( x  3 )( x  1 )( x  4 )  5 x( x  1 ) 
( x  3 )( 4 x 2  4 )
4
4
 ( x  3 )( x 2  1 )  ( x  3 )( x  1 )( x  1 )  x 3  3x 2  x  3
83
x0  0 , y 0  3
x1  1,
y1  0
x 2  3,
y2  0
x3  4 ,
y3  15
LAGRANGE POLYNOMIALS
The Lagrange interpolation polynomial that passes through the points :
x
0
1
3
4
y
3
0
0
15
P3( x )  x3  3x 2  x  3
30
25
20
15
10
5
0
-5
-10
-15
-2
-1
0
1
2
84
3
4
5
LAGRANGE POLYNOMIALS
Example : Find the Lagrange interpolation polynomial that passes through the points :
x
0
3
4
5
y
8
-4
0
18
L0 ( x ) 
( x  x1 )( x  x2 )( x  x3 )
( x  3 )( x  4 )( x  5 )  1


( x  3 )( x  4 )( x  5 )
( x0  x1 )( x0  x2 )( x0  x3 ) ( 0  3 )( 0  4 )( 0  5 ) 60
L1( x ) 
( x  x0 )( x  x2 )( x  x3 )
( x  0 )( x  4 )( x  5 ) 1

 x( x  4 )( x  5 )
( x1  x0 )( x1  x2 )( x1  x3 ) ( 3  0 )( 3  4 )( 3  5 ) 6
L3 ( x ) 
( x  x0 )( x  x1 )( x  x2 )
( x  0 )( x  3 )( x  4 ) 1

 x( x  3 )( x  4 )
( x3  x0 )( x3  x1 )( x3  x2 ) ( 5  0 )( 5  3 )( 5  4 ) 10
P3 ( x ) 
L2 ( x ) 
3
y2 L2 ( x )  y3 L3 ( x )
 yk Lk ( x )  y0 L0 ( x )  y1L1( x )  
k 0
x0  0,
0
8
4
18

( x  3 )( x  4 )( x  5 )  x( x  4 )( x  5 )  x( x  3 )( x  4 )
60
6
10

1
2( x  3 )( x  4 )( x  5 )  10 x( x  4 )( x  5 )  27 x( x  3 )( x  4 )
15

1 3
2 x  24 x 2  94 x  120  10 x 3  90 x 2  200 x  27 x 3  189 x 2  324 x
15

1
 15 x 3  75 x 2  30 x  120  x 3  5 x 2  2 x  8
15
85




y0  8
x1  3,
y1  4
x2  4,
y2  0
x3  5,
y3  18
LAGRANGE POLYNOMIALS
The Lagrange interpolation polynomial that passes through the points :
x
0
3
4
5
y
8
-4
0
18
P3( x )  x3  5x2  2x  8
35
30
25
20
15
10
5
0
-5
-1
0
1
2
3
86
4
5
6
ORDINARY DIFFERENTIAL EQUATIONS
A differential equation is an equation that contains derivatives of an unknown
function.
A differential equation that has one independent variable is called an
ordinary differential equation (ODE).
If an ODE involves only first derivatives of the dependent variable (y) with
respect to the independent variable (x), it is a first – order ordinary
differential equation.
A first –order ODE is linear, if it is a linear function of y and dy/dx.
dy
 ax 2  by  0
dx
( linear )
dy
 ayx  b y  0
dx
( nonlinear )
87
EULER’S METHOD
Numerical Solution
yi+1
Exact Solution
y(x)
Slope =f(xi, yi)
yi
h
xi+1
xi
dy
 f ( x, y )
dx
Slope 
dy
 f ( xi , yi )
dx x  x
i
xi 1  xi  h
yi 1  yi  h f ( xi , yi )
88
EULER’S METHOD
Example : Use Euler’s method to solve the Ordinary Differential Equation below from x=0 to x=1
with the initial conditions x=0 and y=2. (Take h=0.2) Compute the errors in each step
for the exact solution y=3ex-x-1
dy
 x y
dx
f ( x, y )  x  y
x0  0
0  x 1
y( 0 )  2
y0  2
xi 1  xi  h
yi 1  yi  h f ( xi , yi )
x1  x0  h  0  0.2  0.2
y1  y0  h f ( x0 , y0 )  2  0.2 f ( 0,2 )  2  0.20  2  2.4
x2  x1  h  0.2  0.2  0.4
y2  y1  h f ( x1 , y1 )  2.4  0.2 f ( 0.2, 2.4 )  2.4  0.20.2  2.4  2.92
x3  x2  h  0.4  0.2  0.6
y3  y2  h f ( x2 , y2 )  2.92  0.2 f ( 0.4, 2.92 )  2.92  0.20.4  2.92  3.584
x4  x3  h  0.6  0.2  0.8
y4  y3  h f ( x3 , y3 )  3.584  0.2 f ( 0.6, 3.584 )  3.584  0.20.6  3.584  4.4208
x5  x4  h  0.8  0.2  1
y5  y4  h f ( x4 , y4 )  4.4208  0.2 f ( 0.8, 4.4208 )  4.4208  0.20.8  4.4208  5.465
89
EULER’S METHOD
h=0.25
xi
0
0.2500
0.5000
0.7500
1.0000
h=0.1
yi
yi
Abs.
(Euler) (Exact) Error
2.0000
2.5000
3.1875
4.1094
5.3242
2.0000
2.6021
3.4462
4.6010
6.1548
0
0.1021
0.2587
0.4916
0.8306
h=0.2
xi
0
0.2000
0.4000
0.6000
0.8000
1.0000
xi
0
0.1000
0.2000
0.3000
0.4000
0.5000
0.6000
0.7000
0.8000
0.9000
1.0000
h=0.05
yi
yi
Abs.
(Euler) (Exact) Error
2.0000
2.2000
2.4300
2.6930
2.9923
3.3315
3.7147
4.1462
4.6308
5.1738
5.7812
2.0000
2.2155
2.4642
2.7496
3.0755
3.4462
3.8664
4.3413
4.8766
5.4788
6.1548
yi
yi
Abs.
(Euler) (Exact) Error
2.0000
2.4000
2.9200
3.5840
4.4208
5.4650
2.0000
2.4642
3.0755
3.8664
4.8766
6.1548
0
0.0642
0.1555
0.2824
0.4558
0.6899
90
0
0.0155
0.0342
0.0566
0.0832
0.1146
0.1517
0.1951
0.2459
0.3050
0.3736
xi
0
0.0500
0.1000
0.1500
0.2000
0.2500
0.3000
0.3500
0.4000
0.4500
0.5000
0.5500
0.6000
0.6500
0.7000
0.7500
0.8000
0.8500
0.9000
0.9500
1.0000
yi
yi
(Euler) (Exact)
2.0000
2.1000
2.2075
2.3229
2.4465
2.5788
2.7203
2.8713
3.0324
3.2040
3.3867
3.5810
3.7876
4.0069
4.2398
4.4868
4.7486
5.0261
5.3199
5.6309
5.9599
2.0000
2.1038
2.2155
2.3355
2.4642
2.6021
2.7496
2.9072
3.0755
3.2549
3.4462
3.6498
3.8664
4.0966
4.3413
4.6010
4.8766
5.1689
5.4788
5.8071
6.1548
Abs.
Error
0
0.0038
0.0080
0.0126
0.0177
0.0232
0.0293
0.0359
0.0431
0.0510
0.0595
0.0687
0.0788
0.0897
0.1015
0.1142
0.1280
0.1429
0.1590
0.1763
0.1950
EULER’S METHOD
6.5
6
5.5
5
4.5
4
3.5
3
2.5
2
0
0.1
0.2
0.3
0.4
0.5
91
0.6
0.7
0.8
0.9
1
EULER’S METHOD
Example : Use Euler’s method to solve the Ordinary Differential Equation below from x=0 to x=1.8 with
the initial conditions x=0 and y=1. (Take h=0.6) Compute the errors in each step for the exact solution.
dy
 yx  x3
dx
0  x  1.8
f ( x , y )  yx  x 3
x0  0
2
yexact  x 2  e0.5x  2
y( 0 )  1
y0  1
xi 1  xi  h
yi 1  yi  h f ( xi , yi )
x1  x0  h  0  0.6  0.6


y1  y0  h f ( x0 , y0 )  1  0.6 f ( 0,1 )  1  0.6 1 0  03  1
x2  x1  h  0.6  0.6  1.2


y2  y1  h f ( x1 , y1 )  1  0.6 f ( 0.6,1 )  1  0.6 1 0.6  0.63  1.2304
x3  x2  h  1.2  0.6  1.8


y3  y2  h f ( x2 , y2 )  1.2304  0.6 f ( 1.2,1.2304 )  1.2304  0.6 1.2304  1.2  1.23  1.0795
92
EULER’S METHOD
h=0.6
xi
0
0.6000
1.2000
1.8000
h=0.1
yi
yi
Abs.
(Euler) (Exact) Error
1.0000
1.0000
1.2304
1.0795
1.0000
1.1628
1.3856
0.1869
0
0.1628
0.1552
0.8926
h=0.3
xi
0
0.3000
0.6000
0.9000
1.2000
1.5000
1.8000
yi
yi
Abs.
(Euler) (Exact) Error
1.0000
1.0000
1.0819
1.2118
1.3203
1.2773
0.8395
1.0000
1.0440
1.1628
1.3107
1.3856
1.1698
0.1869
0
0.0440
0.0809
0.0989
0.0652
0.1075
0.6526
xi
0
0.1000
0.2000
0.3000
0.4000
0.5000
0.6000
0.7000
0.8000
0.9000
1.0000
1.1000
1.2000
1.3000
1.4000
1.5000
1.6000
1.7000
1.8000
yi
yi
(Euler) (Exact)
1.0000
1.0000
1.0099
1.0293
1.0575
1.0934
1.1355
1.1821
1.2305
1.2778
1.3199
1.3518
1.3675
1.3587
1.3157
1.2255
1.0718
0.8337
0.4841
1.0000
1.0050
1.0198
1.0440
1.0767
1.1169
1.1628
1.2124
1.2629
1.3107
1.3513
1.3787
1.3856
1.3620
1.2955
1.1698
0.9634
0.6481
0.1869
93
Abs.
Error
0
0.0050
0.0099
0.0147
0.0192
0.0235
0.0272
0.0303
0.0323
0.0329
0.0314
0.0269
0.0181
0.0033
0.0201
0.0557
0.1084
0.1855
0.2972
EULER’S METHOD
1.4
1.2
1
0.8
0.6
0.4
0.2
0
0
0.2
0.4
0.6
1
0.8
94
1.2
1.4
1.6
1.8