Engineering Numerical Analysis

Third Year / Surveying Eng. Dept "Engineering and Numerical Analysis" Assist. Prof. Dr.
Qais Fadhil
‫لطلبة الصف األول‬
‫إعداد‬
‫مىفق يحيى عثمان‬
‫أستاذ مساعد‬
‫قسم صحة المجتمع‬

‫المعهد التقني – كركىك‬
1
Third Year / Surveying Eng. Dept "Engineering and Numerical Analysis" Assist. Prof. Dr. Qais Fadhil
2
Subject Target students Hours in a week

Units
Engineering and Third year Theory Practical Total
Numerical Analysis students 3 - 3 6
The Aims
1. General aims:
This course will provide the student with a foundation in
mathematical principles and to solve the differential equations with
classic mathematics and numerical methods.
2. Special aims: The students can be able to;

A – Develop a mathematical relation for a given physical problem,
B – Solving the mathematical equations to get the unknown variables,
C – Choosing the appropriate solution method for the given differential
equations,
D – Solving first, second and higher order of differential equations,
using different method in classical mathematics and numerical
analysis, and,
E – Understanding the numerical methods that can help him in solving
non-standard mathematical equations using approximate methods
on computers, calculators or even on hands.
:‫تىزيع الدرجات خالل السنة‬
.)‫ (نصف السنة‬.‫ درجة امتحان فصلي نظري‬20 ‫ ويتضمن‬:‫ الفصل األول‬-
.‫ درجة امتحان فصلي نظري‬20 ‫ ويتضمن‬:‫ الفصل الثاني‬-
)‫ حضور الطالب‬+ ‫ درجات (امتحانات سريعة‬10 :‫ أعمال السنة‬-
)‫ درجة‬50 ‫(فيكون السعي السنوي من‬
‫ درجة‬50 :‫ االمتحان النهائي‬-
)‫ درجة‬011 ‫(فتكون الدرجة النهائية من‬
3
Syllabus: :‫المنهاج‬
week Topics
0 First Order Differential Equations
2 Methods of Solving First Order D.E.
3 Applications on First Order D.E.
4 Second Order Linear D.E. with Constant Coefficients
5 Applications on Second Order D.E.
6 Higher Order Linear D.E.
7 Applications on Higher Order D.E.
8 Integrating Factors
9 Fourier series for Periodic Functions
10 Even and Odd Functions
11 Applications on Fourier Series
12 Power Series
13 Solution of ODE by Power Series
14 Numerical Solution of Nonlinear D.E.
15 Simple Iteration Method
16 Newton-Raphson Method
17 Interpolation
18 Lagrange Interpolation
19 Solution of Linear Instantaneous D.E.
20 Direct and Indirect Method
21 Numerical Integration
22 Numerical Differentiation
23 Solving Partial D.E.
24 Applications
25 Newton Method for Curve Editing
26 Solving Ordinary D.E.
27 Runge-Kutta Method
28 Power Series
29 Exponential Equations
30 Frobenos Method
References: :‫المصادر‬
1- "Advanced Engineering Mathematics", by Erwin Kreyszig.

2- "Advanced Engineering Mathematics", by C. Ray Wylie.
3- "Mathematics - Schaum's Outline of Differential Equations", by Richard Bronson.
4
5
Most Important Trigonometric Relations:
Note: Review “Integration by Part”.
6
Engineering Analysis
Using Modeling
can be solved
using either or:
Classical Mathematics
(First term of year) Numerical Methods
(Second term of year)
Definitions:
# Differential Equation (D.E.): An equation involving one or more derivatives of a

function beside the function itself and the constants and the independent
variables.
#Solution of D.E.: An equation relates between the dependent and independent

variables, and is free of derivatives and satisfies the D.E.
d
For example; d ( ) is a D.E., and ∫ ( ) is a solution of this D.E. (where x
is independent variable and y is dependent variable).
Example: Is y(x)=c1sin2x + c2cos2x, where c1 and c2 are arbitrary constants, a

solution of the D.E.: y'' + 4y =0 ?
Solution: Differentiating y, we find y' = 2c1cos2x - 2c2sin2x and y'' = -4c1sin2x -
4c2cos2x. Therefore;
y''+4y = (-4c1sin2x -4c2cos2x) + 4(c1sin2x+c2cos2x)
= (-4c1+4c1)sin2x + (-4c2+4c2)cos2x = 0
(Thus y(x)=c1sin2x+c2cos2x satisfies the given differential equation for all values of
x and is a solution for it).
# Ordinary D.E.: (O.D.E.): A differential equation is an ordinary differential equation if

the unknown function depends on only one independent variable..
# Partial D.E.: (P.D.E.): If the unknown function depends on two or more independent
variables, the differential equation is a partial differential equation.
7
Examples on O.D.E. are:
while example on P.D.E. is:
Note: at present, we shall be concerned only with

ordinary differential equations (O.D.E.)
# Order of D.E.: Is the order of the highest derivative that appears in the D.E.
# First Order D.E.: A differential equation which contains only in addition to y and
given functions of x and constants.
# Linear D.E.: When the dependent variable (like y) and its derivatives ( ́ ) are
of first degree.
Therefore, the general linear ordinary D.E. of order n must take the form:
( ) ( ) ( ) ( ) ( )
and any D.E. cannot be put in this form is called "Nonlinear D.E.".
Examples:
…………………… Linear first order O.D.E.
………………… …Linear second order O.D.E.
( ) . . Non-linear third order O.D.E.(y'''y')
…… Linear second order O.D.E.
………..Non-linear second order O.D.E.(y y')
…………… …… Non-linear second order O.D.E.(sin y)
……………………Non-linear second order O.D.E.(y' y'')
( ) …………………Non-linear first order O.D.E.(y')2
………………Linear second order O.D.E.
8
# General solution of D.E.: (G.S.): When solving the D.E. to find the function that
satisfies the D.E. The solution with the unknown constants of (integration or
solution) is called "General Solution".
# Particular solution of D.E.: (P.S.): When the values of the unknown constants of
(integration or solution) are determined using initial or boundary conditions,
the solution is called "Particular Solution".
Therefore, are General Solutions,

while are Particular Solutions, because we
determined the values of c = 3, c1 = -4.2, and c2 = 5.8.
Note: The problem which contains initial conditions is called

"initial value problem", while the problem which contains
boundary conditions is called "boundary value problem".
Note that:
a) Method of Solving First-Order O.D.E.:
The standard form for a first-order D.E. in the unknown function y(x) is:
y' = f(x,y) or dy/dx = f(x,y) (like dy/dx = x3 + 2y = f(x,y))
where the derivative y′ appears only on the left side. Many, but not all, first-order
differential equations can be written in standard form by algebraically solving for y′ and
then setting f (x,y) equal to the right side of the resulting equation.
9
1) Separable Differential Equations:
A separable D.E. is an equation which can be written in the form:
A(x) dx + B(y) dy = 0 (i.e. separating x and y components)

If we integrate both sides; Then the solution will be;
∫ A(x) dx + ∫ B(y) dy = c where c is an arbitrary constant.
d
Example#1: Solve the D.E. :
d
Solution: This equation can be written in the form:
(x2+2) dx – y dy = 0 which is separable with A(x) = x2+2 and B(y) = -y.
Integrating both sides, gives: ∫(x2+2) dx - ∫y dy = ∫0 then solving gives;
solving for y, we obtain (with a new constant k = -2c)
OR √
Note: It is better to find the solution in terms of y = f(x)
d √
Example#2: Solve the D.E.:
d
Solution: separating variables and integrating both sides:
∫ ∫
√
| | ( | | ) G.S.

Solution: separating variables and integrating both sides:
∫9 y dy =∫ -4 x dx or 9 y2/2 = -2 x2 + c1
√
or G.S. (where c2 = c1/2)
11

d d
Solution: d , separating variables:
d
integrating both sides; ∫ ∫ ,
therefore; OR ( ) ………..(G.S.)
Example#5: Solve the initial value problem: ( )
Solution: , separating variables:
integrating both sides;
Now, taking e on both sides; .G.S. (where c2=ec1)

applying initial condition, that x = 0, y = 1, and substituting in the D.E.;
1=c2 e0 OR c2=1; therefore; ……P.S.
2) Exact First Order O.D.E.:
If the D.E. has the form: M(x,y) dx + N(x,y) dy =0 then this D.E. is called exact D.E.
if it satisfies the exactness condition that .
Note: If either N or M in the D.E. contains (sin, cos, tan, tanh,

sinh,….., log, ln, e) it may be means that this D.E. is not exact.
Now, if the exactness condition is satisfied, solve the D.E. using u=∫ M dx + k(y)
where, u=u(x,y),and k(y) is constant of integration, then after this integration derive
w.r.t y or ( ) to determine k(y) using .
Example#1: Solve the D.E. (x3 +3x y2)dx+(3x2 y +y3)dy=0.
Solution: first test for exactness condition: M(x,y)= x3+3xy2 and N(x,y)= 3x2y+y3
and;
11
which are equal, so the D.E. is exact.
so, ∫ ( ) ∫( ) ( ) ( ) …(1)
d ( )
To find k(y): d
d ( )
cancelling similar terms from the two sides; ( )
d
substituting in equation (1) we get; ( ) ( ) ... G.S.
Example#2: Solve the initial value problem:

(sinx coshy)dx - (cosx sinhy)dy =0, y(0)=3.
Solution: M(x,y) = sinx coshy and N(x,y) = - cosx sinhy
and;
which are equal, so the D.E. is exact.

so, ∫ ( ) ∫( ) ( ) ( ) …(1)
d ( )
To find k(y): d
d ( )
cancelling similar terms from the two sides; ( )
d
substituting in (1) we get; ( ) ….. G.S.

Now, using the given initial condition that x=0, y=3
cos0 cosh3 = c = 10.07
therefore, cosx coshy = 10.07 ….. P.S.
Example#3: Find if the D.E. y dx – x dy = 0 is exact or not.
Solution: Here; M(x,y) = y and N(x,y) = -x, and;
which are not equal, so the D.E. is not exact! Therefore, to change it to exact we need
to multiply it by a function called "Integrating Factor" (I.F.)
Now, multiply the given non-exact D.E. by I.F.=( ), to get a new equation:
12
O.K. exact D.E.
The question is "How can we get the I.F.??!"

The following equations are used to get I.F.;
( ) ( )
I.F.(x) = ∫ ( )
( ) ( )
OR I.F.(y) = ∫ ( )
Note: In mathematics: b (b)a and; a b

ba
Example#4: Solve the following D.E. using exact method:

2sin(y2)dx+ xy cos(y2)dy=0, y(2)=(π/2)0.5.
Solution: M(x,y)= 2 sin(y2), N(x,y)= xy cos(y2)
( ) ( ) …. NOT EXACT
( ) ( )
Try I.F.; I.F.(x) = ∫ ( )
( ) ∫ ( ( ) ( )) ∫ , Here ignore c
c ( )
( ) ( )
therefore,
multiply the D.E. by x3, we get; 2x3 sin(y2)dx + x4y cos(y2)dy = 0
( ( )) ( ) ( ( )) ( ) o.k.
Now, the equation is exact.
∫ ( ) ( ) ( ), ( ) ( ) ( )
therefore; k'(y) = 0 and k(y) = constant = c
and; ( ) ( ) …… G.S.
Substituting for the initial condition that x = 2, y = (π/2)0.5 into the G.S., we get:
therefore, ( ) ( ) …… P.S.
13
3) Linear First Order O.D.E.:

The first order O.D.E. is said to be linear if it can be written in the form:
y' + p(x) y = r(x) where p & r are any given functions of x.
If r(x) is zero: r(x)=0 "Homogeneous D.E."

If r(x) is not zero: r(x) ≠0 "Non-homogenous D.E."
1) Homogenous D.E.:
y' + p(x) y = 0 this equation can be easily solved by separating variables;
d
( ) ∫ ( )
taking e of both sides; ( ) ∫ ( )d G.S. (where c = e c* ) (Ready Formula)
2) Non-homogenous D.E.:
Here; y' + p(x) y = r(x)
d d
OR d ( ) ( )
d
comparing with Mdx + Ndy = 0, M = py - r, N = 1
, this means that this D.E. is not exact, to change it to exact multiply
it by I.F.;
( ) ∫ ( ) ∫ ( ) ∫
( ) ∫ d
∫ d ( ) ∫ d
multiply the two sides by I.F. (x);
∫ d ∫ d
or ( )
∫ d ∫ d
integrating w.r.t. x ∫
∫ d
divide both sides by and putting ∫
( ) [∫ ] ….. G.S. for non-homogenous equation. (Ready Formula)

14
2x
Example#1: Solve the following linear D.E.: y' - y = e
Solution: Here, p = -1, r = e 2x, h = ∫p dx = -x

( ) ∫ …. G.S.
Also, we can solve it using I.F.(= e-x). Take it as a H.W.
Example#2: Solve the following initial value D.E.:

y' + y tanx = sin2x , y(0) = 1
Solution: Here, p = tanx, r = sin 2x = 2 sinx cosx, h = ∫p dx = ∫tanx dx = ln(secx)
eh = secx, e-h = cosx, eh r = (secx)(2sinx cosx) = 2 sinx
therefore, y(x) = cosx [2 ∫sinx dx +c] = c cosx – 2 cos2x …. G.S.
Now, using the initial condition, x = 0, y = 1;
1 = c*1-2*12, thus c = 3, therefore, y(x) = 3 cosx – 2 cos2x …. P.S.
Note: sin(x+y) = sin(x)cos(y) + cos(x)sin(y), so

sin(2x) = sin(x+x) = sin(x)cos(x)+cos(x)sin(x) = 2sin(x)cos(x).
4) Bernoulli First Order O.D.E.:
A Bernoulli D.E. is a non-linear D.E. which can be written in the form;

y' + p(x) y = g(x) y n (where n is any real number)
If n = 0 or 1, this D.E. is linear (type 3 previously), otherwise change it to linear D.E. by
letting y1-n = z
d
Example#1: Solve the D.E.
d
Solution: This D.E. is not separable, nor exact, try to change it to Bernoulli eqn.
d d
( ) … (1)
d d
15
this is similar to y' + p(x) y = g(x) yn

d d
let y1-3=z OR z=y -2 d d
d
multiplying equation (1) by (-2y - 3), we get; d
d
OR ( ) … (2), linear in z, p(x)= -8/x
d
∫ d
multiply equation (2) by (x - 8), gets; (x - 8 dz – 8x - 9 z dx) = - 4 x - 4 dx

∫ d (x - 8.z) = ∫ -4 x - 4 dx
H.W.:
A) Find the G.S. for the following D.E.s:
1) 2 x  3 y  2 dx  3x  4 y  1 dy  0
2) y   y sin x  e
cos x
3) 2 x  3 y  2 dx  3x  4 y  1 dy  0
B) Find the P.S. for the following D.E.s:
1) sinh x cos y dx  coshx sin y dy  0 , y 0  0 .

2) 2 xy dy  x 2  y 2 dx , y 1  2
C) Solve the following initial value D.E. using exact method:


2 sin( y 2 ) dx  x y cos( y 2 ) dy  0 y ( 2) 
2
16
APPLICATIONS OF FIRST ORDER O.D.E.S:
Example#1: The tank in figure contains 200 gallons of water in which 40 lb of salt are
dissolved. Five gallons of brine, each containing 2 lb of dissolved salt, run into the
tank per minute, and the mixture, kept uniform by stirring, runs out at the same
rate. Find the amount of salt y(t) at any time t.
Solution: The time rate of change of salt is:

Brine of 5 gal/min
y' = dy/dt = salt inflow rate - salt outflow rate.
If y(t) is the amount of salt in the tank at any
time, and the amount of water is always 200 gal
in the tank, then each 1 gal. contains y(t)/200 Ib 200 gal. of water
of salt. And the 5 gals outflow contains:
5*y(t)/200 = 0.025 y(t) Ib of salt at each minute.
y' = 10 - 0.025y … (a)
by separating variables: , multiplying the left side by (-0.025/-0.025)
and integrating gives; ln(10-0.025y) = -0.025t + c*

taking exponentials for the two sides; 10-0.025y = ce-0.025t
OR …. G.S.
Note: the initial amount of salt in the tank is y(0)=40 Ib (given above);
Thus; y = 40 = 400 – 40 ce0, c = 360/40 substituting in the G.S.;
Now, y(t) = 400 - 360 e-0.025t (Ib) (P.S.) (amount of salt in the tank at any time t)
Example#2: Assume that the tank in Ex.#1 above contains 1000 gallons of water in
which 200 lb of salt are dissolved. Fifty gallons of brine, each containing (1+cos t)
lb of dissolved salt, run into the tank per minute, and the mixture, kept uniform by
stirring, runs out at the same rate. Find the amount of salt y(t) at any time t.
Solution: The time rate of change of salt is: y' = dy/dt = salt inflow rate –
salt outflow rate.
OR y' = In - Out
In = 50*(1+cos t) lb/min.
17
Out = 50(y(t)/1000) = 0.05 y lb/min.

so, y' = 50(1+cos t)-0.05y
y'+0.05y = 50(1+cos t) (non-homogenous linear D.E.)
Here p = 0.05, h = 0.05t, and r = 50(1+cos t)
(∫ ( ) )
solving and substituting for the initial condition;

y(t) = 1000 + 2.494 cost + 49.88 sint - 802.5 e-0.05t …. P.S. (Take it as a H.W.)
Example#3: Suppose that you turn off the heat in your home at night 2 hrs before you
go to bed; call this time (t = 0). If temperature T at t = 0 is 66 0F and at the time you
go to bed (t = 2) has dropped to 63 0F. What temperature can you expect in the
morning, say 8 hrs later (t = 10)? This process of cooling off will depend on the
outside temperature TA, which is assumed to be constant at 32 0F.
Solution: Newton's Law of Cooling shows that the time rate of change dT/dt of
temperature T of a body is proportional to the difference between T and the
outside temperature TA.
dT/dt = k (T-TA) = k (T-32), where k is unknown constant of proportionality.
d
separating variables; , integrating; ln(T-32) = k t + c1
taking e for the two sides; T(t) = 32 + cekt …. G.S.

using the initial condition that T(0) = 66
T(0) = 32+c = 66; c = 34, and T(t) = 32 + 34 ekt …. P.S.
we can also determine the constant k, using T(2) = 63 (another given condition)
T(2) = 32+34 ek.2 = 63, thus e2k = 0.911 765, and k = 0.5 ln(0.911 765) =-0.046 187
thus, T(10) = 53.4oF (temperature at morning)
H.W.: The tank shown contains 2000 gallons of water in which 350 Ib of salt are dissolved. Thirty
gallons of brine, each containing 1  e  0.1t

cos 2t Ib of dissolved salt, run into the tank per minute.
The mixture, which kept uniform by stirring, runs out at the same rate. Find the amount of salt y(t) in
the tank at any time t.
30 gal brine
per minute
2000 gal water 30 gal of mixture

350 Ib salt per minute
18
b) Method of Solving Second-Order O.D.E.:
1) Linear Second Order O.D.E.:

A second order D.E. is called linear if it can be written as;
( ) ( ) ( )
If r(x) = 0; ( ) ( ) (homogenous D.E.)
If r(x) ≠ 0; ( ) ( ) ( ) (non-homogenous D.E.)
1) Homogenous Equations:
For second order homogenous linear D.E., a general solution will be of the form:
y = c1 y1 + c2 y2 , where c1 & c2 are two arbitrary constants, and y1 & y2 are two
suitable solutions. Now, for the initial value problem, there must be two initial
conditions: y(xo) = ko, y'(xo) = k1
with constant coefficients:
Where p(x) and q(x) are constants and not functions of x, the D.E. is called "an
equation with constant coefficients";
Here a&b are "constant coefficients".
in operator form: D2 + aD + b = 0
characteristic equation: λ2 + aλ + b = 0
and the roots are: ( √ ), & ( √ )
Here, we have three cases:
In this case y1 = e λ1 x & y2 = e λ2 x

and y = c1y1+c2y2 = c1 e λ1 x + c2 e λ2 x
19

2y'' + 2y' - 12y = 0, y(0) = 10, y'(0) = 0
Solution: dividing by 2 (to make the y''-factor = 1) gives; y'' + y' - 6y = 0

In operator form: D2 + D – 6 = 0
Characteristic equation: λ2 + λ – 6 = 0
( √ b) ( √ )
( √ b) ( √ )
therefore, y = c1y1 + c2y2 = c1 e 2 x + c2 e -3 x …. G.S.
using the initial conditions: y(0) = 10, & y'(0) = 0 gives;
10 = c1e 2(0)+c2e - 3(0) = c1+c2; gives: c1 = 10-c2 …. eqn (a)
also; y' = 2c1 e 2 x - 3c2 e -3 x
0 = 2c1 e 2 (0) - 3c2 e -3 (0) = 2c1 – 3c2 …. eqn (b)
substituting for (a) into (b) gives c2 = 4, and c1 = 6
Now; y = 6 e 2 x + 4 e -3 x … P.S.
Here, λ1= λ2= λ = -a / 2, and y = c1y1+c2y2 = c1 eλ x + c2 x eλ x (i.e. multiply the second term by x)
Example#1: Solve the following D.E.: y'' + 8y' + 16y = 0
Solution: Characteristic equation: λ2 + 8λ + 16 = 0
( √ b) ( √ )
( √ b) ( √ )
(here a2 - 4b = 0) and λ1 = λ2 = λ = -a / 2 = -4
therefore, y = c1y1 + c2y2 = c1 e -4 x + c2 x e -4 x …. G.S.
21

y'' - 4y' + 4y = 0, y(0) = 3 & y'(0) = 1
Solution: Characteristic equation: λ2 - 4λ + 4 = 0
The roots are; λ1 = λ2 = λ = 2
y = c1y1+c2y2 = c1 e 2 x + c2 x e 2 x …. G.S.
using the initial conditions: y(0) = 3, y'(0) = 1
3 = c1e 2(0)+c2(0)e 2(0)
; gives: c1 = 3
also; y' = 2c1 e 2 x + c2 e 2x + 2c2 x e 2 x
1 = 2c1 e2 (0) + c2 e 2(0)+ 2c2 (0) e 2 (0) , 2c1+c2 = 1, c2 = -5
Now; y = 3 e 2 x - 5 x e 2 x … P.S.
Using √ and making the term a2- 4b positive by pulling out (-1) from the root
and putting √ under the roots, gives:
( √ ) ( √ ) ( √b )
and similarly for λ2 gives: ,where √b
( )
( )
a a
( )
( )
add and divide by 2 gives
subtract and divide by 2 gives
( ) …. G.S.

y'' + 0.2y' + 4.01y = 0, y(0) = 0 & y'(0) = 2
Solution: Characteristic equation: λ2 + 0.2λ + 4.01 = 0
21
( √ ) ( √( ) )
( √ ) ( √( ) )
and; w = 2 ( ) …. G.S.
using the first initial condition y(0) = 0 gives; y(0) = A = 0
using chain rule to get the derivative:

( )
from the second initial condition; y'(0) = 2, gives; y'(0) = 2B = 2; B = 1
…. P.S.
Summary of cases I,II,III

Case Roots Basis of solution General Solution
I Distinct real roots λ1&λ2 y = c1 eλ1 x + c2 eλ2 x
II Real double root λ= -a/2 a a y = (c1 +c2x) e -a x/2
Complex conjugate roots
a a
III (
a
)
Example#2: Solve the following boundary value D.E.:

y'' + y = 0, y(0) = 3 & y(π) = -3
Solution: The basis for the solution is y1 = cos x & y2 = sin x ,(because here a = 0)
and, y = c1 cos x + c2 sin x …. G.S.
first B.C. gives y(0) = c1 = 3, and the second B.C. gives y(π) = c1 cosπ + c2 (0) = -3
where, cos π = -1 and c1 = 3, therefore; y = 3 cos x + c2 sin x (No condition for c2)
H.W.#1: Solve the B.V. problems:

(1) y''+4y = 0, y(0) = 3, y(π/2) =-3. Ans. y = 3 cos 2x + c2 sin 2x
(2) y''+2y'+2y = 0, y(0) = 1, y(π/2) = 0 Ans. y = e –x cos x
22
APPLICATIONS OF SECOND ORDER HOMOGENOUS

LINEAR O.D.E.S WITH CONSTANT COEFFICIENTS:
- Free Oscillations (Mass-Spring Systems):
1) Un-damped System:
Hooke’s law: The restoring force F of a spring is
equal and opposite to the forces applied to the
spring and is proportional to the extension
(contraction) of the spring as a result of the
applied force; that is, F = −kl, where k denotes
the constant of proportionality, generally called
the spring constant, l is the extension (or Δ).
For example A steel ball weighing 128 lb is suspended from a spring, and the spring is
stretched 2 ft from its original length. Thus, F = -128 lb. Hooke’s law then gives -128 =
-k(2), or k = 64 lb/ft. OR simply k = F/Δ = 128/2 , so; k = 64 lb/ft.
To determine the equation of motion y(t) for this mechanical system;
From Newton's 2nd law; m*a = F OR m y'' = F
where F is the resultant of all forces acting on the body
at time t. Assuming downward direction as +ve, when Un-stretched
spring So
stretching the spring, an upward force is Fo = +k. So -a-
(Hooke's law) F=k.So y
System in static
Now, from equilibrium in vertical direction at static equilibrium
-b-
System in
case (b): ∑Fv = 0, -F + W = 0 motion
therefore, -k So + mg = 0, -c-
and, k So = mg W=mg
Also from equilibrium in stretched position (c): ∑Fv = m av k.y
W-k So-k y = m (d2y/dt2), but W = k So gives;

k.So
my''+ ky = 0 (governing D.E.)
(2nd order linear homogenous D.E.)
23 W=mg
The G.S. of this D.E. is a "harmonic oscillation": ()

a
where; √ √b
Example#1: If an iron ball of weight w = 89 N stretches a spring 10 cm in static case.

What would be the equation of motion if the spring is stretched 15 cm additional
downward distance and leaved to oscillate with this weight. Neglect the spring
weight.
Solution: Here, from Hooke's law, F = kδ

89 N = k*10 cm, OR k = 890 N/m, and mass m = W/g = 89/9.81 = 9.08 kg
and the D.E. is: my'' + ky = 0 9.08y'' + 890y = 0 (second order linear hom. D.E.)
characteristic equation: 9.08λ2 + 890 = 0 OR λ2 + 98.126 = 0
OR λ2 = -98.126 (complex roots).

and the general solution is: y = e -ax/2 (A cos wx + B sin wx)
a
y = A cos wt + B sin wt, where; √b √ (a = 0)
therefore, y = A cos 9.9t + B sin 9.9t …. G.S.

using initial conditions, y(0) = 0.15 m, gives A = 0.15 ,and y'(0) = 0, gives B = 0
therefore, y = 0.15 cos 9.9t (m) …. P.S.
2) Damped System:
If we connect the mass to a dashpot, we have to take the
corresponding viscous damping into account. The damping force
(Fd) has the direction opposite to the instantaneous motion, and it
is proportional to the velocity y' of the mass.
Now, Fd = -c y' (negative because it is opposite to motion)
Dashpot c
where, c is damping constant and has units of mass/second.
Damped System
24
k.y
now, ∑Fv = m av
+W - k So – ky - cy' = my'', but W = k So k.So
my''+ cy' + ky = 0 (governing D.E.)

(2nd order linear homogenous D.E.) W=mg
After dividing by (m), the corresponding c y'
characteristic equation is:

c
c
the roots are; √
c
and, √
Here we have three cases:
I/ Over damping; c2-4mk > 0 OR c2 > 4mk, gives "distinct real roots" (λ1, λ2)
t t
Here c is so large that c2 > 4mk, y(t) = c1 e λ1 + c2 e λ2 …G.S.
The body here do not oscillate.
II/ Critcal damping; c2-4mk = 0 OR c2 = 4mk, gives a "real double root" (λ)
Here B = 0 and λ1=λ2= -
and , y(t) = (c1 + c2 t) e –a t / 2 …. G.S.
The body here can have at most one passage through the equilibrium position.
III/ Under damping; c2-4mk < 0 OR c2 < 4mk, gives " conjugate complex roots".
This is the most interesting case; here c is so small that c2 < 4mk, then;
c
√ √
c a
gives;
and y(t) = e –αt (A cos w*t + B sin w*t) …. G.S.

(this equation represents the damped oscillation)
25
Example#1: Solve the previous example (iron ball) but if the system has damping
constant given by:
1. c = 200 kg/sec.
2. c = 179.8 kg/sec.
3. c = 100 kg/sec.
Solution: The governing D.E. is: my'' + cy' + ky = 0
1) For c = 200 kg/sec., (here c2 > 4mk) (over damping)

9.08 y'' + 200 y' + 890 y = 0, y(0) = 0.15, y'(0) = 0
the roots are λ1 = -11.01+4.882=-6.19, λ2 = -11.01-4.882 = -15.83
t t
therefore, y(t) = c1 e λ1 +c2 e λ2 = c1 e -6.19 t +c2 e -15.83 t …. G.S.
substituting for y(0) = 0.15; c1+c2 = 0.15, and y'(0) = 0, λ1c1+λ2c2 = 0
gives, c1 = 0.2463, c2 = -0.0963
t t
Now, y(t) = c1 e λ1 +c2 e λ2 = 0.2463 e -6.19 t – 0.0963 e -15.83 t …. P.S.
2) For c = 179.8 kg/sec., (here c2 = 4mk) (critical damping)

and, λ = -9.899, solving for real double roots and using initial condition:
c1 = 0.15, c2 = 1.485, and y(t) = (0.15+1.485 t) e -9.899 t …. P.S.
3) For c = 100 kg/sec., (here c2 < 4mk) (under damping)

and, λ1 = -5.506+8.227i λ2 = -5.506-8.227i, solving for complex conjugate roots and
using initial condition: A = 0.15, B = 0.1004,
and y(t) = e -5.506 t (0.15 cos8.227t+0.1004 sin8.227t) …. P.S.
2) Non-homogenous Equations:
y'' + p(x) y' + q(x) y = r(x) (non-homogenous 2nd order linear D.E.)
The general solution of this D.E. is of the form: y(x) = yh (x) + yp (x)
where [yh] is homogenous solution = c1 y1(x) + c2 y2(x)
and [yp] is any particular solution for the right hand side of the D.E. (r(x)),
and the final particular solution is obtained by getting values for c 1&c2 in yh (x).
26
# Method of Undetermined Coefficients:

If the D.E. has the form y''+ay'+by = r(x) with constant coefficients a&b, the particular
solution of the right hand side [r(x)] will be as follows:
Example#1: Solve the initial value problem for the non-homogenous equation:
y'' + 2y' + 101y = 10.4ex, y(0) = 1.1, y'(0) = -0.9
Solution: characteristic equation of the homogenous part is;

λ2+2λ+101= 0, and the roots: λ1 = -1+(1-101)0.5 = -1+10 i & λ2 = -1-(1-101)0.5 = -1-10 i
and the general solution of the homogenous part is:
yh = e –x (A cos10x+B sin10x), and to obtain the general solution of the whole D.E., we
have to assume any suitable particular solution yp of the non-homogenous part. Since ex
exist on the right hand side, assume;
yp=c ex substitute in the D.E. gives:(1+2+101) c ex = 10.4 ex, and 104 c ex = 10.4ex
so, 104c = 10.4, gives c = 0.1
now, y(x) = yh+yp = e –x (A cos10x+B sin10x)+0.1 ex …. G.S.
and from initial condition that y(0) = A+0.1 = 1.1, gives A = 1
now, y'(x) = e –x (-cos10x –B sin10x-10 sin10x +10B cos10x)+0.1e x
and y'(0) = -1+10B+0.1 = -0.9, gives B = 0
therefore, y'(x) = e –x cos10x+0.1 e x …. P.S.
Choice for yp in method of undetermined coefficients is:
Term in r(x) Choice for yp

k e λx c e λx
k xn (n=0,1,2,3,…..) cn xn +cn-1 xn-1 +…..+c1 x +ko
k cos wx
k sin wx L cos wx + M sin wx
k e λ x cos wx
k e λ x sin wx e λ x (L cos wx + M sin wx)
27
Example#2: Solve the non-homogenous equation: y'' + 4y = 8 x2
Solution: The general solution of the homogenous part is:
yh = A cos2x + B sin2x
For the non-homogenous part, suggest the choice: yp = k2x2+k1x+ko, then y''p = 2k2
Substitute in the D.E.: 2k2+4(k2 x2+k1 x+ko) = 8 x2
equating coefficients of x2, x1, and xo on both sides gives,
4k2 = 8, 4k1 = 0, 2k2 + 4ko = 0, thus, k2 = 2, k1 = 0, ko = -1
therefore, yp = 2x2-1, and the general solution is:
y(x) = yh+yp = A cos 2x + B sin 2x + 2x2-1
Example#3: Solve the non-homogenous equation: y'' - 3y' + 2y = ex
Solution: The characteristic equation for the homogenous part is;
λ2 -3λ+ 2 = 0, solving gives; λ1 = 1 & λ2 = 2, and yh = c1 ex + c2 e2x
For yp choose yp = k ex. But we see that ex is a solution of yh, therefore our new
suggestion must be; yp = k x ex, thus y'p = k(ex+x ex) & y''p = k(2 ex+x ex)
substituting in the D.E. gives; k(2+x)ex -3k(1+x)ex +2kxex =ex
and –k = 1 OR k = -1, and the general solution is;
y = c1 ex + c2 e2x – x ex …. G.S.
Example#4: Solve the initial value problem: y'' + 2y' + y = e–x, y(0) = -1, y'(0) = 1

λ2 +2λ+ 1 = 0, OR (λ+1)2 = 0, solving gives; λ = λ1 = λ2 = -1 (double root)
therefore, yh = (c1+c2 x) e –x
28
Assume yp= k e–x. But e–x is a solution of yh, therefore, assume yp = k x e–x which is
also a solution of yh. Therefore assume yp = k x2 e–x, and y'p = k(2x-x2)e-x, and y''p = k(2-
4x+x2)e-x
substituting gives; k(2-4x+x2)e-x +2k(2x-x2)e-x +kx2 e-x =2k e-x =e-x, gives k = 1/2
therefore, y = (c1+c2x)e-x + 0.5x2 e-x …. G.S.
the first initial conditions gives y(0) =c1 = -1, differentiating gives;
y'(0) = c2 - c1 = 1, OR c2 = 0
therefore, y = (0.5 x2 -1) e-x …. P.S.
Example#5: Solve the initial value problem:

y'' + 2y' + 5y = 1.25e 0.5x + 40 cos4x - 55 sin4x, y(0) = 0.2, y'(0) = 0

λ2 +2λ+ 5 = 0, and the roots are; λ1 = -1+2 i, λ1= -1-2 i
and yh = e-x (A cos2x+ B sin2x)
choosing yp=L e 0.5x + K cos4x + M sin4x,
thus; y'p = 0.5L e 0.5x - 4K sin4x + 4M cos4x
and; y''p = 0.25L e 0.5x - 16K cos4x - 16M sin4x, substituting in the left hand side;
(0.25+1+5) Le 0.5x+(-16K+8M+5K)cos4x+(-16M-8K+5M)sin4x = right hand side
This gives; 6.25L = 1.25, and L = 0.2
and; -11K+8M = 40 & -8K-11M = -55, this gives; K = 0 & M = 5
therefore, y = e –x (A cos2x + B sin2x) + 0.2e 0.5x + 5 sin4x …. G.S.
using the initial conditions; y(0) = A+0.2 = 0.2, gives, A = 0
Differentiating gives, y'(0) = e –x (-B sin2x + 2B cos2x)+ 0.1e 0.5x +20 cos4x = 0
Therefore, y = -10.05 e -x sin2x + 0.2 e 0.5x + 5 sin4x …. P.S.
29
APPLICATIONS OF SECOND ORDER NON-HOMOGENOUS

LINEAR O.D.E.S WITH CONSTANT COEFFICIENTS:
- Forced Oscillations (Mass-Spring Systems):
# For free motion: k.y

W-k So-k y-c y' = m y'' , but W = k So
k.So
so, m y'' + c y' + k y = 0 (Homogenous D.E.)
where; m y'' : force of inertia,
c y' : damping force, and, W=mg
k y : spring force.
c y'
# For forced motion:

W-k So-k y-c y' +F = m y'' , but W=k So
so, m y'' + c y' + k y = r(t) (Non-homogenous D.E.)
where; F = r(t) : driving force, or "input". The corresponding solution is called
"response" of the system of the driving force.
For periodic inputs at the cases of sinusoidal force, r(t)=Fo cos wt
and the governing D.E. will be: m y'' + c y' + k y = Fo cos wt
Here also, y(t)=yh(t)+yp(t), and using the method of undetermined coefficients:
yp(t) = a cos wt + b sin wt
and y'(p) = -w a sin wt + w b cos wt, y''(p) = -w2 a cos wt - w2 b sin wt
31
Example#1: A 16 lb object stretches a spring (8/9) ft by itself. The system is attached to

a damper which exerts a force of 12 lbs when the velocity is 2 ft/sec. The spring is
initially displaced 6 inches upwards from its equilibrium position and given an
initial velocity of 1 ft/sec downward. Find the displacement y(t) at any time t.
Solution: To find mass (m), and spring constant (k);

b
slug
f ec
lb/ft
To find the damping coefficient, (here we ignored any sign)

( ) ( ), then the governing D.E. is:
OR; ( ) ( )
Solving gives a real double root , and benefitting from the two initial conditions;
then the P.S. is: ( )
H.W: Repeat the same example above but with 5 lbs exerted force when the velocity is
√
2 ft/sec. Answer: ( ) ( )
31
Joseph Fourier: (1768-1830) French physicist and mathematician form Paris, accompanied
Napoleon to Egypt. He developed the theory of heat conduction (heat equation) using these
series which then became a most important tool in mathematical physics.
# Fourier Series: Is a technique to represent a complex periodic functions in

terms of series of simple functions like sine and cosine.
# Periodic Function: A function f(x) is called periodic if it is defined for all x such
that:
f(x+p) = f(x), where p any positive number and called a "period" of f(x) or
length of interval.
"sine and cosine are periodic functions"

"x, x2, x3, ex, coshx, lnx are non-periodic functions"
Therefore, for sinx & cosx the period is 2π, and for sin2x&cos2x the period is π,
and for sin3x&cos3x the period is 2π/3.
so, our problem is to represent
f(x)
various functions of period 2π in
terms of simple functions of
terms:(1, cosx, sinx, cos2x,
sin2x, cos3x, …., cosnx, sinnx).
x
summing up the series:
p
ao+a1cosx+b1sinx+a2cos2x+
Periodic function
b2sin2x+….=f(x)
OR ∑ ( b ) ( ) (complex periodic function of
period 2π)
- Determination of constant ao:
( ) ∑ ( b )
Integrating on both sides from –π to π:
∫ ( ) ∫ ∑ ( b )
∫ ( ) ∫ ∑ ( ∫ b ∫ )
This gives; ∫ ( )
32
- Determination of constant an of the cosine terms:
Again multiplying the equation by (cos mx) and integrating from –π to π:

∫ ( ) ∫ ∑ ( b )
∫ ∑ [ ∫
zero b ∫ ]
∫ ( ) n = 1,2,3,….
- Determination of constant bn of the sine terms:
Again multiplying the equation by (sin nx) and integrating from –π to π:

b ∫ ( ) n=1,2,3,….
Example#1: Find the Fourier series of the periodic function in figure below.
The formula is:
( ) , ( ) ( )
f(x
k )
x
-π 0 π 2π
-k
Solution: ∫ ( ) *∫ ∫ +
| |
This can be seen from that area under curve of f(x) between –π and π is zero.
∫ ( ) *∫ ( ) ∫ +
* | | +
because sin nx = 0 at –π, 0, π for all n = 1,2,3,…
b ∫ ( ) *∫ ( ) ∫ +
c c
* | | + , since cos(-α)=cos α and cos 0=1
b ( ( ) ) ( )
Now, cosπ = -1, cos2π = 1, cos3π = -1,…..etc., so in general:
33
, and thus ,
Therefore; the Fourier coefficients for bn:
b b b b b b
And the Fourier series of f(x) is:
( ) ( )
Now, ( ) ( )
In fact, un-necessary work in determining Fourier coefficients can be avoided

if a function is even or odd.
y
A function y = g(x) is: y = g(x)
 even if g(-x) = g(x) o

x
(here only cosine terms in Fourier Series)
even function (similar to cosine
∫ ( ) ∫ ( ) = area under curve function (symmetry about y-axis)
y
 odd if g(-x)=-g(x) y = g(x)
(here only sine terms in Fourier Series) o x

∫ ( ) = area under curve
odd function (similar to sine
function (symmetry about origin)
The Fourier series for an even function of period 2L is:

( ) ∑ (Fourier half-range cosine series)
where; ∫ ( ) ∫ ( ) , n = 1,2,3,…
And the Fourier series for an odd function of period 2L is:

( ) ∑ b (Fourier half-range sine series)
where; b ∫ ( ) , n = 1,2,3,…
34
Example#1: Find the two Fourier half-range expansions for the periodic
function in figure below. The formula is: f(x)
k
( ) {
( ) x
0 L/2 L
Solution: a) Even periodic extension: (half-range cosine series)

⁄ f(x)
∫ ∫ ( )
k
x
⁄ 0
∫ -L L/2 L
Even extension
∫ ( )
For the first integral, integrating by parts:

⁄
∫ | ∫
( )
Similarly, for the second integral;
∫ ( ) ( ) | ∫
( ) ( )
After cancelling sine terms and L2 terms, an will be:
( )
Thus a1(n=1) = 0, ( ) , a3(n=3) = 0, a4(n=4) = 0, a5(n=5) = 0, ( )
Hence, the cosine half-range expansion of f(x) is:

( ) ( )
b) Odd periodic extension: (half-range sine series)

solving b (H.W.)
Therefore sine half-range expansion of f(x) is:
( ) ( ) f(x)
k
-L 0 L/2 L x
-k
Odd extension
35
Power Series:
Is a standard method for solving linear ODEs with variable coefficients. It gives
solutions in the form of power series. These series can be used for computing values,
graphing curves, proving formulas, and exploring properties of solutions.
So the power series (in powers of x-x0) is:
Example#1: Using Power Series solve: ́

Solution: In the first step we use;
Then substitute
36
….. G.S.
Also, we can generalize the Power Series formula to solve the second order ODE s of
the form:
Then we use;
Example#2: Using Power Series solve;

Solution: we must do Power Series for all the terms;
Add terms of like powers of x (x0, x , x2,….) and equate the sum obtained to zero;
37
Numerical Analysis
Numeric methods: are methods that are suitable for solving, numerically and
approximately, the difficult problems on computers, calculators, or even on hands.
Numerical analysis requires the following steps:

1) Modeling We set up a mathematical model of our problem, such as an integral, a
system of equations, or a differential equation .
2) Choosing a numeric method and parameters (e.g., step size), perhaps with a
preliminary error estimation.
3) Programming We use the algorithm to write a corresponding program in a CAS
(Computer Algebra Systems), such as Maple, Mathematica, Matlab, or Mathcad, or,
in Java, C or FORTRAN, selecting suitable routines from a software system as
needed.
4) Doing the computation.
5) Interpreting the results in physical or other terms, also deciding to rerun if further
results are needed.
Degree of accuracy (D): is number of significant figures after the digit or number of
significant digits. For example the number 3.4148 has 4D (or 4S) accuracy, and the
number 6394.21 has 2D (or 2S) accuracy.
Formulas for Errors: If an approximate value is ̃ quantity whose exact value is a, we

call the difference:
the error of ̃
Hence;
Note: We always work with absolute value of error because its sign is not important.
38
1) Solution of Equation f(x) = 0 by Iteration

If f is a given function of x, then solving for the equation f(x) = 0 means finding x
values that satisfy this equation, this called:
a) The solution of the equation f or,
b) The roots of the equation f or,
c) The points of intersections of function f with x-axis.
Note: If f is of first degree then x has one value, and if f is polynomial

of second degree then x has two values, and so on.
Iteration Method: Is a method to find x-value in which we start from an initial guess
(which may be poor) and compute step by step (in general better and better)
approximations of an unknown solution of f.
This method is done by changing the equation f(x) = 0 to x = g(x), then we choose an xo
to compute x1 = g(xo), x2 = g(x1), and so on. In general:
Example#1: Set up an iteration process for the equation. ( ) ,

starting from x0 = 1, and using 3D accuracy.
Solution: The exact solution is:
The equation f(x) may be written as:
CONVERGES
S
DIVERGES
39
CONVERGES
CONVERGES
An iteration process defined by g1(x) is called convergent for an xo if the corresponding

sequence xo, x1, x2, ….. is convergent.
Note: At convergence, the solution goes into a specific value and more iterations give
more accurate value. Therefore we need (Degree of Accuracy) to stop the iterations.
41
2) Solution of Equation f(x) = 0 by Newton’s Method

Newton’s method, also known as Newton–Raphson’s method, is another iteration
method for solving equations f(x) = 0 where f is assumed to have a continuous
derivative ́ . The method is commonly used because of its simplicity and great speed.
In the second step we compute in the third step we compute

x3 from x2 again by the same formula, and so on.
Example#1: Set up a Newton iteration for computing the square root x of a given
positive number c and apply it to c = 2. Find the approximate solution to 6D
accuracy, starting from x0 = 1.
Solution:
x5 = 1.414214. Note that x4 and x5 are equal to the sixth significant figure, and this
means that the solution to 6D accuracy is x = 1.414214.
Example#2: Apply Newton’s method to the equation f(x)=x 3+x-1=0, use 6D accuracy.
Solution: From Newton’s formula, we get: Starting from x0 = 1
x5 = 0.682328. Note that x4 and x5 are the same to the sixth significant digit, therefore,
the approximate solution to 6D is x = 0.682328.
41
3) Solution of Equation f(x) = 0 by Secant Method

Newton’s method is very powerful but has the disadvantage that the derivative may
sometimes be difficult to get or computationally expensive, therefore the rid of
derivative is a matter of concern by:
and the new formula called secant method will be:
Example#1: Find, to 6D accuracy, the positive solution of f(x) = x-2sinx = 0 by the

secant method, starting from x0 = 2, x1 = 1.9.
Solution:
n xn-1 xn Nn Dn xn-1-xn
1 2.000000 1.900000 -0.000740 -0.174005 -0.004253
2 1.900000 1.895747 -0.000002 -0.006986 -0.000252
3 1.895747 1.895494
4 1.895494 1.895494
Therefore, the solution to 6D accuracy is x = 1.895494.
Now, let:
42
Given (x0,f0), (x1,f1), (x2,f2),…,(xn,fn) with arbitrarily spaced xj, Lagrange had the idea
of multiplying each fj by a polynomial that is 1 at xj and 0 at the other n nodes and then
taking the sum of n+1 polynomials.
1) Linear Lagrange Interpolation:

Is interpolation by a straight line p1(x) through
(x0,f0), (x1,f1)
p1 = L0 f0 + L1 f1
where,
Therefore,
Example#1: Compute a 4D-value of ln9.2 from ln9.0=2.1972, ln9.5=2.2513 by linear

Lagrange interpolation and determine the error, using ln9.2=2.2192 (4D).
Solution: we have x0 = 9.0, x1 = 9.5, f0 = ln9.0, f1 = ln9.5, therefore,
and,
43
2) Quadratic Lagrange Interpolation:

Is interpolation of given (x0,f0), (x1,f1), (x2,f2) by a second degree polynomial p2(x)
which by Lagrange’s idea is:
Example#2: Compute ln9.2 by Quadratic Lagrange Interpolation from the data given in
example#1 and the additional third value ln11.0 = 2.3979, use 4D accuracy.
Solution: Quadratic Lagrange Inter. formula is:
The error formula is:

Which means the error value = exact value – approximate value.
Percentage error is calculated by dividing error value by 100.
3) Newton’s Forward Divided Difference Interpolation: (for equal spacing h)
44
Now, let x = x0 + rh, then r = (x-x0) / h. Finally we get;

( ) ( )( ) ( )( ) ( )
( )
…... Forward Difference Interpolation Formula

where, called "Delta"
The error is:
Example#1: compute cosh0.56 from the four values in the following table and estimate
the error, use 6D accuracy.
Solution: After completing the above table;
The exact value to 6D is cosh0.56 = 1.160941.

Therefore, the error is ε = | |= | |= 0.000003
4) Newton’s Backward Divided Difference Interpolation: (for equal spacing h)
45
Also, let x = x0 + rh, then r = (x-x0) / h. Finally we get;

( ) ( )( ) ( )( ) ( )
( )
…... Backward Difference Interpolation Formula

where, called "Nabla"
Example#2: Compute a 7D-value of the Bessel function J0(x) for x =1.72 from the four
values in the following table, using (a) Newton’s forward formula,
(b) Newton’s backward formula.
Solution: Note that the two methods have the same table completion; the difference is
only in notations.
46
5) Spline Interpolation: (polynomial interpolation)

The given data (xo,fo), (x1,f1),…,(xn,fn) can be
interpolated by a polynomial Pn(x) of degree n or less so
that the curve of Pn(x) passes through these n+1 points.
But if n is large, there may be a trouble in calculating
Pn(x), which may oscillate between the nodes.
So, instead of using a single high-degree polynomial Pn over the entire interval
in which the nodes lie ( ), we use n low-degree one
over each subinterval between adjacent nodes, like qo(x) from xo to x1, then q1(x) from x1
to x2, and so on.
Then from this we compose an interpolation function g(x), called spline, by fitting these
small polynomials together into a single continuous polynomial curve passing through
the data points.
Now, if g(x) is the spline consists of n cubic polynomials qo,q1,…,qn-1. We wright these
polynomials in the form;
where j = 0,1,…,n-1. Using Taylor’s formula, we obtain;
Note that for equidistance nodes,
47
4
Example#1: Interpolate f(x) = x on the interval by the cubic spline g(x)
corresponding to the nodes xo = -1, x1 = 0, x2 = 1 and satisfying the conditions
( ) ( ) ( ) ( )
Solution: We can wright, ( ) ( ) ( ) , so

our spline consists of n =2 polynomials,
48
1A
1A
1B
(A) (B)
Figure (1) (A) Rectangular rule, (B) Trapezoidal rule.
49
Example#1: Using 6D accuracy:
Solution: It is important to do the following table:
…… to (6D) accuracy
51
Example#2: Using 6D acc., evaluate ∫ by Simpson’s rule with n = 10
Solution: Since h = 0.1, we can do the following table:
H.W. #1: Using Simpson's rule and a 5D accuracy, evaluate the integral:
 /2
x
J 
 /6 sin x
dx with n = 4.
H.W. #2: Using Simpson’s rule to 5D-accuracy, evaluate the integral:
 
1
J   1  2 e  x dx
2
with n = 10.
0
H.W. #3: Use the trapezoidal rule to 5D-accuracy to evaluate the integral:
 
1
J   1  e  x dx
2
with n = 10.
0
51
1) Euler (Euler-Cauchy) Method:

We know that the general form for the first order initial value ODE is:
where x0 and y0 are given.

Here, we will discuss methods of computing approximate numeric values of the solution
y(x) at equidistance points on x-axis:
Now, from Taylor series, which is:
And for small h value, h2, h3,… can be neglected, thus:
Example#1: Apply the Euler method to the following initial value problem, choosing
h=0.2 and computing y1,…..,y5, and using 3D accuracy then find the error.
́ ( )
Solution: Here f(x, y) = x + y; hence f(xn, yn) = xn + yn, and we see that:
yn+1 = yn + 0.2(xn + yn)
Now; y1 = yx=0.2 = y0 + 0.2(x0+y0) = 0 + 0.2(0+0) = 0
Also; y2 = yx=0.4 = y1 + 0.2(x1+y1) = 0 + 0.2(0.2+0) = 0.04
y3 = yx=0.6 = y2 + 0.2(x2+y2) = 0.04 + 0.2(0.4+0.04) = 0.128
y4 = 0.274, y5 = 0.489
52
Remember that: the error value = exact value – approximate value.

Percentage error is calculated by dividing error value by 100.
Now, to get the exact solution yexact:

p = -1, r = x , h = ∫p dx = -x
previously we learned that ∫ ∫
use integration by part; let u = x, therefore, du = dx, and let dv = e-x, therefore, v = -e-x
∫ ∫
∫ ∫ ( ) ( )
x -x x
Therefore, y = e [-e (x+1)+c] = c e – x - 1…….. G.S.
x
Now, substitute for x = 0, y = 0 gives c = 1 and y = e – x – 1 ….. is P.S.
Now, Do the following table:
n xn yn (Euler) yn (Exact) Error (%)

0 0.0 0.000 0.000 0.0
1 0.2 0.000 0.021 2.1
2 0.4 0.040 0.092 5.2
3 0.6 0.128 0.222 9.4
4 0.8 0.274 0.426 15.2
5 1.0 0.489 0.718 22.9
53
2) Runge-Kutta (RK) Method:

It is a fourth order method which has greater accuracy than of the Euler method. We
will see that in each step we first compute four auxiliary quantities k1, k2, k3, and k4.
Then the new value yn+1 is computed. Now,
Example#2: Apply Runge-Kutta method to the initial value problem given in the
previous example with the same h = 0.2, and computing five steps to 6D
accuracy, then find the percentage error.
Solution: Here, f(x,y) = x+y and x0 = y0 = 0, h = 0.2

To calculate y1, where x1 = x0 + h = 0 + 0.2 = 0.2:
k1 = 0.2(x0 + y0) = 0.2(0+0) =0
k2 = 0.2[(x0 + 0.5h) + (y0 + 0.5k1)] = 0.2[(0 + 0.1) + (0 + 0.5*0)] = 0.02
k3 = 0.2[(x0 + 0.5h) + (y0 + 0.5k2)] = 0.2[(0 + 0.1) + (0 + 0.5*0.002)] = 0.022
k4 = 0.2[(x0 + h) + (y0 + k3)] = 0.2[(0 + 0.2) + (0 + 0.022)] = 0.0444
Therefore, ( ) …to 6D accuracy.
It is better to do a table as below:
n xn yn (Runge-Kutta) yn (Exact) Error (%)

0 0.0 0.000000 0.000000 0.0000
1 0.2 0.021400 0.021403 0.0003
2 0.4 0.091818 0.091825 0.0007
3 0.6 0.222106 0.222119 0.0013
4 0.8 0.425521 0.425541 0.0020
5 1.0 0.718251 0.718282 0.0031
54
Curve Fitting:
1) Curve Fitting by Straight Line in Least Squares Method:
Where;
Example#1: Using the method of Least Squares fit a straight line for the following four
points:
Solution:
55
2) Curve Fitting by Polynomials of Second Degree in Least Squares Method:
Now, for second degree polynomial (parabola) p(x) = b0 + b1 x + b2 x2
Example#2: Fit a parabola through the data: (0,5), (2,4), (4,1), (6,6), and (8,7).
Solution: For the normal equations we need n = 5, Σxj = 20, Σxj2 = 120, Σxj3 = 800, Σxj4
= 5664, Σyj = 23, Σxjyj = 104, Σxj2 yj = 696. Hence these equations are:
Solving, we obtain the quadratic least squares parabola:
56
Addition and Scalar Multiplication for Matrices:

A Matrix: Is a rectangular array of numbers or functions which enclosed in brackets.
Example#1:
are matrices. The numbers (or functions) inside the matrix are called entries or, less
commonly, elements of matrix. The first matrix in up has two rows, which are the
horizontal lines of entries. Furthermore, it has three columns, which are the vertical
lines of entries. The second and third matrices are square matrices, which mean that
each has as many rows as columns 3 and 2, respectively. The entries of the second
matrix have two indices, signifying their location within the matrix. The first index is
the number of the row and the second is the number of the column, so that together the
entry’s position is uniquely identified. For example, (read a two three) is in Row 2 and
Column 3, etc.
Matrices having just a single row or column are called vectors. Thus, the fourth
matrix has just one row and is called a row vector. The last matrix has just one column
and is called a column vector.
Example#2: We are given a system of linear equations, briefly a linear system, such as:
57
where x1, x2, and x3 are unknowns. We form the coefficient matrix, call it A, by listing
the coefficients of the unknowns in the position in which they appear in the linear
equations.
[ ][ ] [ ]
OR: A x = b
Note: The symbol used for denoting a matrix such as A is either A or
Also from,
This new matrix is called augmented matrix of A. Here we just wrote the column
matrix b at the end of matrix A.
General Notation af a matrix:
Matrix A has m rows and n columns which are called size of the matrix.
Now, for the matrices in Example#1, the

sizes are 2*3, 3*3, 2*2, 1*3, and 2*1.
If m=n, we call A as n*n square matrix.

58
A vector is a matrix with only one row or column. Its entries are called the
components of the vector.
Thus
,
Example#3:
59
Example#4:
Example#5:
Also,
61
Matrix Multiplication:
Example#6:
Example#7:
61
So,
Transposition of a Matrix (T):

We obtain the transpose of a matrix by writing its rows as columns (or
equivalently its columns as rows). This also applies to the transpose of vectors. Thus, a
row vector becomes a column vector and vice versa.
Example#8:
Also,
62
Special Matrices:
63
1) Gauss Elimination Method:
Now, consider the following system of linear equations of triangular form (upper):
Back Substitution: means that we solve the last equation, x2 = -26/13 = -2, and then
work back to obtain x1 = (2-5x2)/2 = 6
Example#1:
(1)
Solution:
(2)
64
Second and Third Order Determinants:
65
Inverse of a Matrix:
Example#1:
Example#2: Find the inverse of the matrix,
66
Solution:
Therefore;
2) Gauss-Seidel Iteration Method:

Example#1:
We can write it as follows,
Least Square Method for Curve Fitting
67
68
Determination of Eigenvalues and Eigenvectors:

Consider multiplying the nonzero vector (x) by a given square matrix (A),
* +* + * +
and this can be written as; * +* + * + , where called
“Eigenvalue” and x is called “Eigenvector” and the problem is called “Eigenvalue

Problem”. Therefore;
Example 1: Determine the Eigenvalues and the Eigenvectors of the matrix;
* +
Solution: (a) Eigenvalues, which must be determined first;
This can be written in matrix form;
because;
The solution of this quadratic (2 nd order) equation gives two Eigen values;
.
(b1) Eigenvectors of A corresponding to . This vector is obtained by substituting
with (arbitrary), and solve for x1 and x2:
This is done by assuming x1 = t (arbitrary assumption) gives x2 = 2t, so the x’s corresponding to is:
69
(b2) Eigenvectors of A corresponding to . This vector is obtained by substituting

with , and solve for x1 and x2:
And also assuming x1 = 2t gives x2 = -1 gives the x’s corresponding to :
Example 2: Determine the Eigenvalues and the Eigenvectors of the matrix;
[ ]
Solution: the determinant gives the cubic equation;
To find Eigenvectors, we apply Gauss elimination method to the system ( ) , first with λ = 5,
and then with λ = -3,
”Good Luck”
71

Engineering Numerical Analysis

Uploaded by

Copyright:

Available Formats

Engineering Numerical Analysis

Uploaded by

Document Information

Original Description:

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Engineering Numerical Analysis

Uploaded by

Copyright:

Available Formats

Third Year / Surveying Eng. Dept "Engineering and Numerical Analysis" Assist. Prof. Dr.

‫لطلبة الصف األول‬

‫قسم صحة المجتمع‬

Subject Target students Hours in a week

2. Special aims: The students can be able to;

:‫تىزيع الدرجات خالل السنة‬

1- "Advanced Engineering Mathematics", by Erwin Kreyszig.

Most Important Trigonometric Relations:

Note: Review “Integration by Part”.

# Differential Equation (D.E.): An equation involving one or more derivatives of a

#Solution of D.E.: An equation relates between the dependent and independent

Example: Is y(x)=c1sin2x + c2cos2x, where c1 and c2 are arbitrary constants, a

# Ordinary D.E.: (O.D.E.): A differential equation is an ordinary differential equation if

Examples on O.D.E. are:

while example on P.D.E. is:

Note: at present, we shall be concerned only with

Therefore, are General Solutions,

Note: The problem which contains initial conditions is called

a) Method of Solving First-Order O.D.E.:

1) Separable Differential Equations:

A separable D.E. is an equation which can be written in the form:

A(x) dx + B(y) dy = 0 (i.e. separating x and y components)

solving for y, we obtain (with a new constant k = -2c)

Note: It is better to find the solution in terms of y = f(x)

Example#3: Solve the D.E. :

Example#4: Solve the D.E. :

Example#5: Solve the initial value problem: ( )

Solution: , separating variables:

integrating both sides;

Now, taking e on both sides; .G.S. (where c2=ec1)

2) Exact First Order O.D.E.:

if it satisfies the exactness condition that .

Note: If either N or M in the D.E. contains (sin, cos, tan, tanh,

w.r.t y or ( ) to determine k(y) using .

Example#1: Solve the D.E. (x3 +3x y2)dx+(3x2 y +y3)dy=0.

which are equal, so the D.E. is exact.

substituting in equation (1) we get; ( ) ( ) ... G.S.

Example#2: Solve the initial value problem:

Solution: M(x,y) = sinx coshy and N(x,y) = - cosx sinhy

which are equal, so the D.E. is exact.

substituting in (1) we get; ( ) ….. G.S.

Example#3: Find if the D.E. y dx – x dy = 0 is exact or not.

Solution: Here; M(x,y) = y and N(x,y) = -x, and;

O.K. exact D.E.

The question is "How can we get the I.F.??!"

Note: In mathematics: b (b)a and; a b

Example#4: Solve the following D.E. using exact method:

Solution: M(x,y)= 2 sin(y2), N(x,y)= xy cos(y2)

Now, the equation is exact.

therefore; k'(y) = 0 and k(y) = constant = c

3) Linear First Order O.D.E.:

If r(x) is zero: r(x)=0 "Homogeneous D.E."

taking e of both sides; ( ) ∫ ( )d G.S. (where c = e c* ) (Ready Formula)

comparing with Mdx + Ndy = 0, M = py - r, N = 1

( ) [∫ ] ….. G.S. for non-homogenous equation. (Ready Formula)

Solution: Here, p = -1, r = e 2x, h = ∫p dx = -x

Also, we can solve it using I.F.(= e-x). Take it as a H.W.

Example#2: Solve the following initial value D.E.:

1 = c1-212, thus c = 3, therefore, y(x) = 3 cosx – 2 cos2x …. P.S.

and y(t) = e –αt (A cos wt + B sin wt) …. G.S.