UoG Numerical method (CEng2073) Worku.

M (Structural Engineer )

Chapter-Three

3. Curve fitting

Introduction

The purpose of curve fitting is to specify an appropriate function and adjust its parameters in a
way that it matches, as close as possible, an experimental or historical data set.

3.1. Least square regression

The Least Squares Regression Line is the line that makes the vertical distance from the data
points to the regression line as small as possible. It is called a “least squares” because the best
line of fit is one that minimizes the variance (the sum of squares of the errors).

Let a set of data points ( x i , y i ) i = 1, 2, 3 … m, be given and let the curve given by y = f( x ) be
fitted to this data. The error of approximation at x =xi is

e i= y i−f ( x i )

S=(( y 1 −f ( x1 ))2 +( y 2 −f ( x 2 ))2 +( y 3 −f ( x 3 ))2 +. ..+( y m−f ( x m ))2 )

S=e21 +e 22 +e 23 +. . .+ e2m Then the method of least square consists in minimize for the best fit.

3.1.1. Linear regression

Let y = ax+b the straight line to be fitted to the given data then

S=(( y 1 −( a+bx 1 ))2 +( y 2 −(a+ bx 2 ))2 +( y 3 −(a+bx 3 ))2 +. ..+( y m−(a+bx m ))2 )

m
S=∑ ( y i −(a+bx i ))2 .. .. .. . .. .. . .. .. . .. .. ... .. .. ...... .. ... .. .. . .. .. . .. .. . .. .. .. .....(1 )
i=1

∂s ∂s
=0 ∧ =0
To minimize S set ∂a ∂b

m
∂s
=∑ −2∗( y i−( a+bx i ))=0
∂a i=1 divide by -2

1 2012 EC.
 UoG Numerical method (CEng2073) Worku.M (Structural Engineer )

m m
∂s
=∑ y −am−b ∑ x i=0. .. . .. .. . .. .. . .. .. . .. .. . .. .. . .. .. .. . .. .. . .. .. . .. .. . .. .. . .. .. . .. ..(2 )
∂a i=1 i i=1

m
∂s
=∑ −2 x i∗( y i −( a+bx i ))=0 divide by (−2)
∂b i=1

m m m
∂s
=∑ xi y i −a ∑ x i −b ∑ x 2i )=0 .. .. . .. .. . .. .. . .. .. . .. .. . .. .. . .. .. .. .(3 )
∂b i=1 i=1 i=1

From Eq. (2) and Eq. (3) we can get a system of equation

am+b ∑ x i =∑ y i
a ∑ x i +b ∑ x 2=∑ x i y i
i is called normal equation & can be written in matrix form

[∑ m
xi
∑ xi
∑ x i2 ]{ } { }
a = ∑ y i .. .. . .. .. . .. .. . .. .. . .. .. .. . .. .. . .. .. . .. .. . .. ..( 4 )
b ∑x y
i i

∂2 s ∂2 s
2
+ 2
Remark: ∂a ∂ b are both positive at point a and b hence this value provides a minimum
value of ‘S’

Example 1: Find a linear Equation for the following data and approximate the value of at x = 3

x 0 2 5 7
y -1 5 12 20
Solution:

x y x2 xy
0 -1 0 0
2 5 4 10
5 12 25 60
7 20 49 140

∑ 14 36 78 210
Using Equation (4)

2 2012 EC.
 UoG Numerical method (CEng2073) Worku.M (Structural Engineer )

[ 4 14 a
=
]{ } { }
36
14 78 8 210
⇒ 4 a+14=36 ∧ 14 a+78 b=210
Solving those equation simultinously
and we have a = -1.13 & b = 2.89, thus y = a + bx meaning

y = -1.13+2 . 89x is the desired equation

Y (3 ) = -1. 13+2. 89*3=7 .54

Therefore, the value of the function at x = 3 is 7.54

Study the error b/n the approximate function and exact output of y

Exact
x output(y) y= f ( x )=-1. 13+2 .89x
0 -1 -1.13
2 5 4.65
5 12 13.32
7 20 19.1
o
Example 2: The table below gives the temperature (T) in degree centigrade ( c) and the length of
o
heated concrete cube, which have 150 mm at 20 c . Fined a linear equation L = a+bT that best fit
the given data.

o 20 30 50 70 90 110
T( c )
L(mm 150 150.05 150.1 150.5 150. 151
) 4 8
Solution

o
T( c ) L(mm) T2 TL
20 150 400 3000
30 150.05 900 4501.5
50 150.14 2500 7507
70 150.5 4900 10535
90 150.8 8100 13572
110 151 12100 16610

3 2012 EC.
 UoG Numerical method (CEng2073) Worku.M (Structural Engineer )

55725.
sum 370 902.49 28900 5
Using Eq.(4) and m = 6

[ 3706 370
]{ } {
a 902 . 49
=
28900 b 55725 . 5 }
⇒ a=149. 6856 ∧ b=0 .0118 therefore L=149 .6856 +0 .0118 T

is the linear equation which fit the given data.

Exercise 1: The table below gives the force applied on steel road and a total length of the road
after the application of force. (Submitted)

a) Fined a linear equation L = a+bF that best fit the given data.
b) Estimate the length after the 20000N application of force

F(kN) 50 55 70 75 90 110 150

L(mm) 1000 1000.4 1000. 1001 1001.2 1001.4 1002
9
 Practice also using MAT-LAB: Curve fitting tool using fit (curve fit)

3.1.2. Fitting a polynomial

Let the polynomial of nth degree

y=ao +a1 x +a 2 x 2 +a3 x 3 +.. .+a n x n be fitted to the data points ( x i , y i ) , i = 1, 2, 3, 4…,m

m
∂S
S=∑ ( y i −( a o + a1 x i +a2 x 2 +a3 x 3 +. ..+a n x n ))2 =0
i i i i
Setting
∂a i we obtain

ma0 +a1 ∑ x i + a2 ∑ x 2i + a3 ∑ x 3i + a4 ∑ x 4i + .. .. . .+ an ∑ x ni =∑ yi

a 0 ∑ x i +a1 ∑ x 2i +a2 ∑ x 3i +a3 ∑ x i4 +. .. . ..+an ∑ x ni =∑ x i y i

……………………………………………………………………

…………………………………………………………………….

a 0 ∑ x ni + a1 ∑ x ni +1 + a2 ∑ x ni +2 +a 3 ∑ x ni +3 +. .. . ..+ an ∑ x ni +n =∑ x ni y i

4 2012 EC.
 UoG Numerical method (CEng2073) Worku.M (Structural Engineer )

In matrix form, we can write as follow:

[ ]{ } { }
m ∑ xi ∑ x 2i .. . ∑ xni a0 ∑ yi
∑ xi ∑ x 2i ∑ x 3i .. . ∑ x n+1
i a1 ∑ xi yi
. . . . . . =.
. . . . . . .
∑ x ni ∑ x n+1
i ∑ x n+2
i ∑ x n+3
i ∑ x n+n
i
an ∑ xn yi i
……………. (5)

Example 3, Find a polynomial of 2nd degree for the data given below.

x 0 1 2
y 1 6 7
Solution: Let y = a0+a1x+a2x2

We require the quantity ∑ ∑ x 2 , ∑ x3 , ∑ x 4 , ∑ y i , ∑ y i x i , ∑ y i x 2

xi , i i i i

x y x2 x3 x4 xy x2y
0 1 0 0 0 0 0
1 6 1 1 1 6 6
2 17 4 8 16 34 68

∑ 3 24 5 9 17 40 74
Then eq. (5) becomes

[ ]{ } { } [
∑ x i ∑ x2
]{ } { }
m i ao ∑ yi 3 3 5 ao 24
∑ xi ∑ x2 ∑ x3i i
a1 = ∑ yi xi ⇒ 3 5 9 a1 = 40
∑ x2 i
∑ x3 ∑ x 4
i i
a2 ∑ yi x2 i
5 9 17 a2 74

2
Then a0 = 1, a1 = 2, a3 = 3 the required polynomial is given by y=1+2 x +3 x

Exercise 2:

x 10 20 40
y 9 15 55
a. Find 2nd degree polynomial of the given data.

5 2012 EC.
 UoG Numerical method (CEng2073) Worku.M (Structural Engineer )

b. find 3rd degree polynomial of the given data (Hint use mat lab or excel for
computation)

3.2. Interpolation

Given a set points (xi, yi) I = 1, 2, 3, …n satisfying the relation y = f(x) where the explicit nature
of f(x) is not known. It is required to find a function J(x) agree at the set of tabulated points such
process is called interpolation.

3.2.1. Forward interpolation

If yo, y1, y2, y3, …,yn a set of value of y, then

y 1 − y 0 ; y 2 − y 1 ;.......... y n − y n−1 are called the
difference of y.

Thus
Δy 0 = y 1− y 0 ; Δy 1= y 2 − y 1 ; Δy 2= y 3− y 2 ;......... Δy n−1 = y n − y n−1 are called the first
forward difference of y.

Δ 2 y 0 =Δy1 −Δy 0 ; Δ2 y 1 =Δy2 −Δy 1 ;.. . .. .. .. ; Δ 2 y n−1 =Δy n−1 −Δy n−2 ; this is called 2nd forward

difference of y.

Similarly we can defined 3rd and higher order difference.

Forward difference table is shown bellow

x y Δ Δ2 Δ3 Δ4 Δ5
x0 y0 Δy0 Δ2 y 0 Δ3 y 0 Δ4 y 0 Δ5 y 0
x1 y1 Δy1 Δ2 y 1 Δ3 y 1 Δ4 y 1
2 3
x2 y2 Δy2 Δ y2 Δ y2
x3 y3 Δy3 Δ2 y 3
x4 y4 Δy 4
x5 y5

Example 4: find the forward difference of (2, 0), (6, -1), (4, 2), (8, 5)

Solution forward difference becomes

6 2012 EC.
 UoG Numerical method (CEng2073) Worku.M (Structural Engineer )

x y Δ Δ2 Δ3
2 0 −1 −4 −4
6 −1 3 0
4 2 3
8 5

Exercise 3: find the forward difference of (2, 4), (6, -7), (4, 20), (8, 50), (10, 34);

3.2.2. Backward difference

The difference
∇ y 1 = y 1− y 0 ; ∇ y 2= y 2− y 1 ; ∇ y 3 = y 3 − y 2 ;.........∇ y n = y n − y n−1 are called the
first backward difference

∇ 2 y 2 =∇ y 2 −∇ y 1=( y 2− y1 )−( y 1 − y 0 )

∇ 2 y 3 =∇ y 3 −∇ y 2=( y 3 − y 2 )−( y 2 − y 1 )

Called 2nd backward difference formula

x y Δ ∇2 ∇3 ∇4 ∇5
x0 y0
x1 y1 ∇ y1
2
x2 y2 ∇ y2 ∇ y2
x3 y3 ∇ y3 ∇ 2 y3 ∇ 3 y3
x4 y4 ∇ y4 ∇2 y 4 ∇ 3 y4 ∇4 y4
x5 y5 ∇ y5 ∇ 2 y5 ∇ 3 y5 ∇ 4 y5 ∇5 y5

Example 5: find the backward difference of (3, 4), (-1,2), (5,0), (8,10)

Solution: by using back difference table

x y ∇ ∇ 2 ∇3
3 4
−1 2 −2
5 0 −2 0
8 10 10 12 12

7 2012 EC.
 UoG Numerical method (CEng2073) Worku.M (Structural Engineer )

Exercise 4: find the backward difference of (2, 4), (6, -7), (4, 20), (8, 50), (10, 34);

3.3. Newton’s derived difference interpolating formula

Given n+1, points (x0, y0), (x1, y1), … (xn,yn) where each xi and equally spaced. It is required to
find a polynomial of nth degree yn(x)such that y and y(x) agree at a tabulated points

Let
x i=x o +h , since y (x) is a polynomial of nth degree it can be written as
n

y n ( x )=a 0 +a1 ( x−x 0 )+a 2 ( x−x 0 )( x−x 1 )+a3 ( x−x 0 )( x−x 1 )( x−x 2 )+. .+a n ( x−x 0 )( x−x 1 )∗. ..( x−x n−1 )

Since y and y n ( x ) are agree at a tabulated point xi, i =1, 2, 3, 4, 5, …….

We have:

2 3 n
y1−y0 Δy0 Δ y0 Δ y0 Δ y0
y0=a lignl¿0;¿¿¿a1= = ¿a2= 2¿a3= 3 . . .∧an= n ¿
x 1− x 0 h 2 ! h 3 ! h n ! h
Δy 0 Δ2 y 0 Δn y 0
y n ( x )= y 0 + ( x−x 0 )+ ( x−x 0 )( x −x1 )+. ..+ ( x−x 0 )( x−x 1 )∗.. .∗( x−x n−1 )
h 2 ! h2 n ! hn

This equation is called newton forward difference formula

Remark if x = x0 + ph

p ( p−1) Δ 2 y 0 p ( p−1)( p−2) Δ 3 y 0 p ( p−1)( p−2)( p−3) Δ 4 y 0

y n ( x )= y 0 + pΔy 0 + + + +. ..
2! 3! 4!
( p ( p−1)( p−2)∗.. .∗( p−n+1 )∗( Δ n y 0 )
+
n!

Example 6: Using newton forward difference formula find the interpolating formula y(x) and y
(1/2)

Where

8 2012 EC.
 UoG Numerical method (CEng2073) Worku.M (Structural Engineer )

x 0 1 2 3
y 2 3 6 9
Solution

The forward difference table is x, y

x y Δ Δ2 Δ3
0 2 1 2 −2
1 3 3 0 Here h =1; therefore the formula
2 6 3
3 9
Δy 0 Δ2 y 0 Δn y 0
y n ( x )= y 0 + ( x−x 0 )+ ( x−x 0 )( x −x1 )+. ..+ ( x−x 0 )( x−x 1 )∗.. .∗( x−x n−1 )
h 2 ! h2 n ! hn

1 2 −2
y n ( x )=2+ ( x−0)+ ( x−0 )( x−1)+ (x −0)( x−1 )( x−2)
1 2∗1 3∗2∗1

−1 2
y n ( x )= x +2 x 2 − x +2
3 3 is the interpolating formula

−1 2 51
y n (1/2)= ∗(1 /2 )+2∗(1/22− (1/2)+2=
3 3 21

Exercise 5: Using newton’s forward formula find the interpolating formula for

x 1 3 5 7 9
y 15 50 76 150 400
and approximate y (10).

Exercise 6: find the missing term in the following

x 0 1 2 3 4
y 1 3 9 ? 81

 Newton’s backward interpolating formula

In this case, we replace x0 by xn, x1 by xn-1, and xk by xn-k in the newton’s forward difference
formula

9 2012 EC.
 UoG Numerical method (CEng2073) Worku.M (Structural Engineer )

∇ yn ∇2 yn ∇ n yn
y ( x )= y n + ( x−x n )+ ( x−x n )( x −x n−1 )+ .. .+ ( x−x n )( x−x n−1 )∗. ..∗( x−x 1 )
h 2 ! h2 n ! hn

x−x n
p=
Remark: if we let h the above formula becomes

p ( p+1 )∇ 2 y n p ( p+ 1)( p+2) ∇ 3 y n p( p+1 )( p+ 2)( p+3 )∇ 4 y n

y ( x )= y n + p ∇ y n + + + +.. .
2! 3! 4!
p( p+2 )( p +2)∗.. .∗( p+ n−1 )∗( ∇ n y n )
+
n!

Example 7: using Newton’s backward difference formula find polynomial for

x 0 1 2 3
y 1 2 1 10
Solution: using backward difference table

x y ∇ ∇2 ∇ 3
0 1
1 2 1
2 1 −1 −2
3 10 9 10 12 Here h=1 , xn = 3
x=x n + ph ⇒ x=3+ p ⇒ p=x−3 ; substitute

( x−3)( x−3+1)∗10 ( x−3 )(( x−3+1)(( x−3+2)∗12

y ( x )=10+( x−3 )∗9+ +
in the formula 2! 3!

y ( x )=2 x 3 −7 x 2 +6 x +1

Remark: Newton forward interpolating formula is used for interpolating near beginning of
tabulated values, but Newton backward interpolating formula is used for interpolating near end
of set of tabulated values.

Exercise 7: using Newton’s backward difference formula find polynomial for

x 0 3 6 9
y 3 10 2 30
3.4. Lagrange Interpolating Polynomial

10 2012 EC.
 UoG Numerical method (CEng2073) Worku.M (Structural Engineer )

Given n+1, points (x0, y0), (x1, y1), … (xn,yn) where each xi not equally spaced. It is required to
find a polynomial of nth degree yn(x) such that y and y(x) agree at a tabulated points. We find a
polynomial of degree n.

Ln (x )=a o + a1 x+ a2 x 2 + a3 x 3 +. ..+ an x n be fitted to the data points ( x i , y i ) , i =0, 1, 2, 3,…,n

Such that Ln(xi) = yi for each i, = 0, 1, 2, 3, 4, ….n

Let,
n
( x−x 1 )( x−x 2 )( x−x 3 )( x−x 4 )∗.. .. .. . ..∗( x−x n ) (x −x j )
Lo ( x )= ⇒∏
( x o −x 1 )( x o−x 2 )( x o −x 3 )( x o− x 4 )∗. .. . .. .. .∗( x o −x n ) j=1 ( x o− x j )
n
( x−x o )( x−x 2 )( x−x 3 )(x −x 4 )∗. .. . .. .. .∗( x−x n ) ( x−x j )
L1 ( x )=
( x1 −x o )( x 1−x 2 )( x 1−x 3 )( x 1 −x 4 )∗. .. . .. .. .∗( x 1−x n )
⇒ ∏ ( x 1 −x j )
j=0∧ j≠1

………………………………………………………………………………………………………
……………………………………………………………………………………………………..

n−1
( x−x o )( x−x 1 )( x−x 2 )( x−x 3 )( x−x 4 )∗.. . .. .. . .∗( x− xn−1 ) ( x−x j )
Ln (x )= ⇒∏
( x n −x o )( x n −x2 )( x n −x 3 )(x n −x 4 )∗.. . .. .. . .∗(x n −x n−1 ) j=0 ( xn −x j )
m
Set , L( x )=∑ Li ( x ) y i =L0 y 0 + L1 y 1 + L2 y 2 +. .. .+ Ln y n
i=0

L(x) is a polynomial of degree n is Called Lagrange Interpolating Polynomial

Example 8: use Lagrange formula to find an interpolating formula for

x 0 1 3 4
y -12 0 12 24

Solution: degree of polynomial =3

11 2012 EC.
 UoG Numerical method (CEng2073) Worku.M (Structural Engineer )

( x−1 )( x−3)( x−4 ) −1

Lo ( x )= ⇒ ( x 3−8 x 2 +19 x−12
(0−1 )(0−3)(0−4 ) 2x

( x−0 )(x −3)( x−4 ) 1

L1 ( x )= ⇒ (x 3 −7 x 2 +12 x )
(1−0 )(1−3 )(1−4 ) 6

( x−0 )( x−1)( x−4 ) −1

L2 ( x )= ⇒ (x 3 −5 x 2 +4 x )
(3−0 )(3−1 )(3−4 ) 6

( x−0)( x−1 )( x−3) 1

L3 (x )= ⇒ ( x 3 −4 x 2 +3 x )
(4−0 )(4−1 )(4−3 ) 12
m
⇒ L n ( x )=∑ Li ( x ) yi =L0 y 0 + L1 y 1 + L2 y 2 +. . ..+ Lm y m
i=0
Ln ( x )=x −6 x 2 +17 x−12
3

Exercise 8: use Lagrange formula to find an interpolating formula for

x 0 3 7 8 9
y -12 0 12 30 25
Exercise 9: The compressive strength of concrete is given in the table for the same ratio of
ingredient ratio. (Hint. Newton interpolation formula)

Number of blows 5 15 25 35 45 55
Compressive strength of concrete 10 12 18 20 26 28
(MPa)

a. Calculate the increasing in compressive strength from the number of blow from 20 to 50.
b. Calculate the maximum compressive strength we obtain using the number of blow less
than 39.

12 2012 EC.
 UoG Numerical method (CEng2073) Worku.M (Structural Engineer )

3.5. Fourier approximations

Fourier showed that an arbitrary periodic function can be represented by an infinite series of
sinusoids of harmonically related frequencies. For a function with period T, a continuous Fourier
series can be written:

f (t )=a 0 +a 1 cos( ω0 t )+b 1 sin(ω 0 t )+a2 cos(2 ω 0 t )+b2 sin (2ω 0 t )+⋯
more concisely
∞
f (t )=a 0 + ∑ [ ak cos(kω 0 t )+b k sin(kω 0 t ) ]
k=1

Where w0=2 π /T is called fundamental frequency and its constant multiples 2w0, 3w0, etc., are
called harmonics.

The coefficients of the equation can be calculated as follows:

T
2
a k = ∫ f (t )cos(kω 0 t )dt
T 0
T
2
b k = ∫ f (t )sin (kω0 t )dt
T 0
for k=1 , 2 , …
T
1
a 0= ∫ f (t )dt
T 0

Example 9: consider a rectangular pulse

13 2012 EC.
UoG Numerical method (CEng2073) Worku.M (Structural Engineer )

 Fundamental period (T) = 2 sec

 Frequency (f) = 1/T = 0.5 Hz
2π 2 π
= =π rad /sec
 Angular frequency (ωo) = T 2 sec

{
1 , 0<t <0 .5
f (t )= 0 , 0. 5<t <1 .5
1 , 1 .5<t <2

First, calculate the average value

ao
T 2
1 1
a o = ∫ f (t )dt= ∫ f (t )dt
T 0 T 0
0. 5 1.5 2
1
a o=
2
∫ 1 dt=12 ∫ 0 dt+12 ∫ 1 dt=0 .25+0+0 . 25=0 . 5
0 0.5 1. 5

 Determine the cosine coefficient

ak
T 0 .5 2
2 2
a k = ∫ f (t )cos(kω o t )dt=
T 0 2
∫ cos(kπt )dt+ 22 ∫ cos(kπt )dt
0 1 .5

ak=
1
kπ
1 2
( sin(kπt )|00 .5 + (sin(kπt )|12 . 5= sin k
kπ kπ
π
2( ) k =1 , 2 , 3 , 4 .. .
The first few values of
ak
14 2012 EC.
UoG Numerical method (CEng2073) Worku.M (Structural Engineer )

2 −2 2
a 1= , a 2=0 , a 3= a4 =0 , a3 = . .. .
π 3π 5π

 a k Zero for all even values of k

 Determine the sine coefficient

bk
T 0. 5 2
2 2 2
b k = ∫ f (t )sin (kω o t )dt= ∫ sin(kπt )dt + ∫ sin (kπt )dt
T 0 2 0 2 1. 5
1 1
b k =− (cos(kπt )|0.0 5 − (cos( kπt )|21. 5 =0 k=1, 2, 3, 4 ...
kπ kπ
∞
f (t )=a 0 + ∑ [ ak cos(kω 0 t )+b k sin(kω 0 t ) ]
Thus: k =1

[ ( ) ]
∞
2 π
f (t )=0 . 5+ ∑ sin k cos( kω0 t )
k =1 kπ 2

15 2012 EC.
UoG Numerical method (CEng2073) Worku.M (Structural Engineer )

Exercise 9: formulate Fourier series for the given equation

{
1 , 0<t <0 .5
f (t )= −1 ,0 . 5<t <1. 5
1 , 1. 5<t <2

16 2012 EC.