View
225
Download
0
Category
Preview:
Citation preview
8/2/2019 Dan Um Jay A
1/560
Numerical Analysis
Course No. : AAOC C341
Dr. P. Dhanumjaya
8/2/2019 Dan Um Jay A
2/560
Instructor-Incharge
Dr. P. Dhanumjaya
Instructor
Dr. Anil Kumar
8/2/2019 Dan Um Jay A
3/560
Text BookAn Introduction to Numerical Analysis
Author: Devi Prasad
Publisher: Narosa Publishing House, 3
8/2/2019 Dan Um Jay A
4/560
Reference BooksAn Introduction to Numerical Analy
Kendall E. Atkinson, John Wiley &(2004).
Elementary Numerical Analysis, S
Conte, Carl de Boor, McGraw-Hill,(1981).
Numerical Analysis, R. L. Burden,
Faires, 7th edition, Thomson Lear
8/2/2019 Dan Um Jay A
5/560
Evaluation SchemTest-1: Weightage 25%,
Date & Time: 20-09-2008 @ 2-3PTest-2: Weightage 25%,Date & Time: 01-11-2008 @ 2-3P
Tut-Test/Quiz/Assignments: WeightagDate & Time:
Comprehensive Exam: Weightage 40Date & Time: 06-12-2008 @ 2-5P
8/2/2019 Dan Um Jay A
6/560
Topics to be coveRoots of equations.
Solve
f(x) = 0,for
x.x)
8/2/2019 Dan Um Jay A
7/560
Linear algebraic equations.
Given the as and the cs, solve
a11x1 + a12x2 = c1,
a21x1 + a22x2 = c2,
for the xs.x2
8/2/2019 Dan Um Jay A
8/560
Curve fitting (Interpolation).
x
x( )f
..
. . .
Interpolation
8/2/2019 Dan Um Jay A
9/560
Integration.
I =ba
f(x) dx.
Find the area under the curve.x( )f
8/2/2019 Dan Um Jay A
10/560
Ordinary differential equations.
Givendy
dt
= f(t, y).
Solve for y as a function of t.y
8/2/2019 Dan Um Jay A
11/560
Why you should study Numerical M
Numerical methods are extremely
problem solving tools. They are cahandling large systems of equationonlinearities and complicated geo
that are often arises in engineeringThey are ingeneral impossible to sanalytically.
8/2/2019 Dan Um Jay A
12/560
During your career, you may often
occasion to use commercially avapackages that involve numerical mThe intelligent use of these packa
predicted on knowledge of the basunderlying the methods.
Many problems cannot be approac
commercially available packages.conversant with numerical method
8/2/2019 Dan Um Jay A
13/560
Finite Digit
Arithmetic and EWe consider how numbers are repres
computers.
Most computers have an integer mod
floating-point mode for representing nThe integer mode is used only to rintegers.
The floating-point form is used to rreal numbers
8/2/2019 Dan Um Jay A
14/560
Note 1Scientific calculations are usually
in floating-point arithmetic.
8/2/2019 Dan Um Jay A
15/560
An n-digit floating-point number in ba
the form
x = (.d1d2d3 dn) e
where (.d1d2d3 dn) is a -fraction cmantissa, and e is an integer called th
A floating-point number is said to be nif d1 = 0 or else d1 = d2 = d3 = = 0
8/2/2019 Dan Um Jay A
16/560
Note 2Most digital computers use the ba
(binary) number system or base 8base 16 (hexadecimal).
Pocket calculators use base 10 (d
8/2/2019 Dan Um Jay A
17/560
Example 1The conversion of a base 2 number to
(11011.01)2 = 1 24 + 1 23 + 0 2+ 1 20 + 0 21 + 1 = 27.25
8/2/2019 Dan Um Jay A
18/560
Floating-Point representations on v
computers
Computer n M =
............................IBM 7094 2 27 2
HP 67 10 10 9
Intel 8087 2 24 12
............................
8/2/2019 Dan Um Jay A
19/560
Note 3The exponent e is limited to a rang
m < e < M
for certain integers m and M.
8/2/2019 Dan Um Jay A
20/560
Chopping and
RoundingMost real numbers x cannot be repre
exactly by the floating-point represent
We approximate by a nearby number
representable in the machine.Given an arbitrary real number x, we denote its floating-point representatio
8/2/2019 Dan Um Jay A
21/560
There are two ways of producing f
chopping and rounding.
Let a real number x be written in the f
x = (0.d1d2d3 dndn+1 )with d1 = 0.
8/2/2019 Dan Um Jay A
22/560
The chopped machine representation
given by
f l(x) = (0.d1d2d3 dn)
8/2/2019 Dan Um Jay A
23/560
The rounded representation of x is giv
f l(x) = (0.d1d2d3 dn)
whenever 0
dn+1 0.
8/2/2019 Dan Um Jay A
105/560
The first approximation to the locationis then
p1 = b1 f(b1) b1 a1f(b1)
f(a1)
=
To determine whether the zero is con(a1, p1) = (1, 1.1) or on (p1, b1) = (1.1,
8/2/2019 Dan Um Jay A
106/560
We calculate
f(p1) = 0.549 < 0.Since f(a1) and f(p1) are of the same
Intermediate Value Theorem tells thabetween p1 and b1.
For the next iteration, we take(a2, b2) = (p1, b1) = (1.1, 2).
8/2/2019 Dan Um Jay A
107/560
The second approximation to the locazero is
p2 = b2 f(b2) b2 a2f(b2)
f(a2)
= 1.15
We calculate
f(p2)
0.274 < 0,
which is of the same sign as f(a2)
8/2/2019 Dan Um Jay A
108/560
Hence, the Intermediate Value Theorthe zero is now between p
2and b
2.
We take
(a3, b3) = (p2, b2) = (1.15174363
In the third iteration, we get
p3 = 1.17684091, and f(p3) = 0
8/2/2019 Dan Um Jay A
109/560
Newtons Methody
Tangentline
y=f(x)
8/2/2019 Dan Um Jay A
110/560
Let pn denote the most recent approxzero p of the function f(x).
Replace f(x) by its tangent line approbased at the location x = pn and take
x-intercept of the tangent line as the napproximation pn+1 to the root.
8/2/2019 Dan Um Jay A
111/560
Since the tangent line approximation x = p
nis given by
y f(pn) = f(pn)(x pn)
Now substituting x = pn+1 and y = 0,
pn+1 = pn f(pn)f(pn)
,
where 0 1 2
8/2/2019 Dan Um Jay A
112/560
NoteWe observe that each iteration of Newmethod, we require two separate funcevaluations:
The evaluation of the function.
The evaluation of the function deri
8/2/2019 Dan Um Jay A
113/560
Example
Find the root of the function
f(x) = x3 + 2x2 3x 1 = on the interval (1, 2) by using the New
method.
Solution. To apply Newtons method, tderivative of f(x) is needed. We get
f ( ) 3 2 4 3
8/2/2019 Dan Um Jay A
114/560
p1 = p0 f(p0)
f(p0) = 1.25,
p2 = p1
f(p1)
f
(p1)
= 1.2009345
p3 = p2 f(p2)f(p2)
= 1.1986958
p4 = p3 f(p3) = 1 1986912
8/2/2019 Dan Um Jay A
115/560
NoteIn the preceding example, Newton
achieved an accuracy of 1.937
1only eight function evaluations, fouevaluations of f(x) and four evaluf(x).
For comparison, starting from the (1, 2), the bisection method needsevaluations and the method of fals
8/2/2019 Dan Um Jay A
116/560
Example
Find the root of the function
f(x) = x3 + 2x2 3x 1,on the interval (3,2) using (i) Bisemethod (ii) Secant method (iii) Regulamethod and (iv) Newtons method.
Compare all the methods withp = 2.9122291785.
8/2/2019 Dan Um Jay A
117/560
Example
Consider the function
f(x) = x3 + 2x2 3x 1.The sequences generated by Newton
for p0 = 1, p0 = 2 and p0 = 3.
8/2/2019 Dan Um Jay A
118/560
Influence of p0
n p0 = 1 p0 = 2
1 1.2500000000 1.4705882353 2.0
2 1.2009345794 1.2471326788 1.4
3 1.1986958411 1.2006987324 1.24 1.1986912435 1.1986949265 1.2
5 1.1986912435 1.1
6 1.1
8/2/2019 Dan Um Jay A
119/560
NoteIn the preceding example, the charesults only in a variation in the nuiterations needed to achive conve
8/2/2019 Dan Um Jay A
120/560
Example
Consider the function
f(x) = tan(x) x 6.Here, if we take p0 = 0.48, the Newton
converges to 0.4510472613 in five itera
If we take p0 = 0.4, then the sequenceby Newtons method fails to converge
5000 iterations.
8/2/2019 Dan Um Jay A
121/560
NoteThus, unlike simple enclosure metNewtons method is not guaranteeconverge for an arbitrary starting p
If we take p0 = 0, then the sequen
converges after 42 iterations to 69This is infact, one of the many zer
f(x) = tan(x)
x
6
8/2/2019 Dan Um Jay A
122/560
NoteThis shows even when Newtons mconverges, it may converge to a vfrom p0.
Order ofC
8/2/2019 Dan Um Jay A
123/560
Convergence
Theorem. Let p be the root of the functf(x) = 0, and assume that f(p)
= 0.
f(x), f(x), and f(x) is continuous inneighborhood of the root p of f(x) the
limn
|en+1||en|2 ,
where =12|f(p)||f(p)| .
W k h
8/2/2019 Dan Um Jay A
124/560
Now taking
en+1 = ppn+1 = ppn + ff
= en +
f(pn)
f(pn).
We calculate
f( ) f ( ) f ( ) f ( )
8/2/2019 Dan Um Jay A
125/560
Since f(p) = 0, we get
f(pn) = f(p en) = enf(p) + en2
We note thatf(pn) = f
(p) + n,
where n 0 as n .
8/2/2019 Dan Um Jay A
126/560
On substituting for f(pn) and f(pn), w
en+1 = en +enf(p) + e
2n
2 f(
f(p) + n
We rewrite the above equation as
en+1 = en +enf(p) + e
2n
2 [f(p)
f(p) + n
8/2/2019 Dan Um Jay A
127/560
We have
en+1 =enn + e
2
n
2 [f(p) + nf(p) + n
Now taking
limn
en+1e2n
=1
2
f(p)
f(p).
8/2/2019 Dan Um Jay A
128/560
Hence,
limn
en+1e2n
= limn
|en+1||en|2 =
12|f(|f(p
Therefore, Newtons method convergquadratically.
E l
8/2/2019 Dan Um Jay A
129/560
Example
The equation
f(x) = x3 + x2 3x 3 = has a root on the interval (1, 2), name
For n 1, compute the ratio |pnp||pn1p|2 athat this value approaches 12
|f(p)||f(p)| .
8/2/2019 Dan Um Jay A
130/560
M
ullers Method
y=f(x)
x3
y
8/2/2019 Dan Um Jay A
131/560
NoteThe methods discussed up to thisus to find a zero of the function onapproximation to that zero is know
These methods are not very satisf
when all the zeros of a function ar
Mullers method is useful for obtai
real and complex roots of a functio
8/2/2019 Dan Um Jay A
132/560
Mullers method is an extension of themethod.
The secant method begins with two inapproximations x0 and x1 and determ
next approximation x2 as the intersecx-axis with the line through (x0, f(x0))(x1, f(x1)).
8/2/2019 Dan Um Jay A
133/560
Mullers method uses three initial appx0, x1 and x2 and determines the nexapproximation x3 by considering the iof the x-axis with the parabola throug(x0, f(x0)), (x1, f(x1)) and (x2, f(x2))
8/2/2019 Dan Um Jay A
134/560
Let
y = a(x
x0)2 + b(x
x0) +
be the equation of the parabola passithe points (x0, f(x0)), (x1, f(x1)) and
The constants a, b and c are determinconditions
f(x2) = (x2 x0)2 + b(x2 x0
8/2/2019 Dan Um Jay A
135/560
Let x1 x0 = h1, x0 x2 = h2, we ob
a = h2f(x1) (h1 + h2)f(x0) + hh1h2(h1 + h2)
b =
f(x1)
f(x0)
ah21
h1 ,
and
c = f(x0).
8/2/2019 Dan Um Jay A
136/560
Now, the points of intersection of the
y = a(x x0)2 + b(x x0) +with x-axis, i.e., y = 0 are
a(x x0)2 + b(x x0) + c =
8/2/2019 Dan Um Jay A
137/560
Due to roundoff errors caused by the of nearly equal numbers, we apply th
x x0 = 2cb
b2
4ac
The sign in the denominator should bthat the denominator will be largest in
8/2/2019 Dan Um Jay A
138/560
With this choice the next approximatioclosest to the zero of f(x).
Therefore,
x3 = x0 2cb + sgn(b) b2 4After finding x3, to begin the next iterause the three recent approximations x
8/2/2019 Dan Um Jay A
139/560
NoteAt each step, the method involves the
b2
4ac, so the method gives appro
complex roots when b2 4ac < 0.
Example
8/2/2019 Dan Um Jay A
140/560
Example
Do three iterations of Mullers methodroot of
f(x) = ex + 1 = 0,
starting with the values 1, 0, 1.
Solution. Clearly f(x) has no real root
We start three initial guesses x0 =
1
and x2 = 1, then we derive the iteratio
8/2/2019 Dan Um Jay A
141/560
The iteration table:
n xn0 1.00
1 0.00
2 1.003 1.0820 + 1.5849i4 1.2735 + 2.8505i
Example
8/2/2019 Dan Um Jay A
142/560
Example
Find the root of the function
f(x) = x3 + 2x2 3x 1,by using the Mullers method. Here ta
starting points x0 = 0.5, x2 = 1 and x1
Solution. Plz try
8/2/2019 Dan Um Jay A
143/560
Fixed Point MethA fixed point of the function g(x) is annumber p for which g(p) = p.
i.e., whose location is fixed by g(x).
8/2/2019 Dan Um Jay A
144/560
Mean Value TheoIf the function f(x) is continuous on tinterval [a, b] and differentiable on theinterval (a, b), then there exists a real (a, b) such that
f() = f(b) f(a)b a .
Proof. Refer any calculus book.
Example
8/2/2019 Dan Um Jay A
145/560
Example
sin 0 = 0,
this implies x = 0 is said to be a fixthe function sin x.
g(x) = x2,
here x = 0, 1 are two fixed points o
Let
8/2/2019 Dan Um Jay A
146/560
0
2
4
6
8
10
Example
8/2/2019 Dan Um Jay A
147/560
Example
Let
g1(x) = ex,
and
g2(x) = x2 + 1.
These two functions do not have any
8/2/2019 Dan Um Jay A
148/560
Question?Under what conditions a function is gto have an unique fixed point?.
8/2/2019 Dan Um Jay A
149/560
TheoremLet g(x) be continuous on the closed [a, b] with g : [a, b]
[a, b]. Then g(x)
point p [a, b].
Furthermore, if g(x) is differentiable o
interval (a, b) and there exists a positik < 1 such that
|g
(x)| k < 1, for all x (
8/2/2019 Dan Um Jay A
150/560
ExistenceAssume that g(x) is continuous on thinterval [a, b] with g : [a, b]
[a, b].
Define the auxiliary function
h(x) = g(x) x.
8/2/2019 Dan Um Jay A
151/560
Note that h(x) satisfies the following p
(i) h(x) is the difference of two functiocontinuous on [a, b], this implies h(x)continuous on that interval.
(ii) the roots of h(x) are precisely the of g(x).
8/2/2019 Dan Um Jay A
152/560
Now, since
minx[a,b] g(x) a, and maxx[a,b] g(xit follows that
h(a) = g(a) a 0, and h(b) = gIf either h(a) = 0 or h(b) = 0, then we
a root of h(x) which is a fixed point of
8/2/2019 Dan Um Jay A
153/560
If neither h(a) = 0 nor h(b) = 0, thenh(b) < 0 < h(a).
Since h(x) is continuous on [a, b] theIntermediate Value Theorem guarant
existence of p [a, b] such that h(p) =implies g(p) = p.
Existence of a fixpoint
8/2/2019 Dan Um Jay A
154/560
point
a
b
y
y= x
y= g
Uniqueness
8/2/2019 Dan Um Jay A
155/560
Uniqueness
This part of the proof will proceed bycontradiction.
Suppose that p and q are both fixed pfunction g(x) on the interval [a, b] with
By the definition of the fixed point g(pg(q) = q.
8/2/2019 Dan Um Jay A
156/560
Then
|p
q|
=|g(p)
g(q)
|.
Using the Mean Value Theorem, we g
|p
q
|=
|g()(p
q)
|= |g()| |p q| k|p q| < |p q
which is a contradiction.
8/2/2019 Dan Um Jay A
157/560
NoteThe hypotheses of this theorem are sconditions but not necessary conditio
It is possible for a function to violate oof the hypotheses, yet still have a fixe
Example
8/2/2019 Dan Um Jay A
158/560
Example
Consider the function
g(x) = 4x(1 x),on the interval [0.1,).Here,
limx
g(x) ,here g does not map [0.1,
) onto its
lim |g(x)| +
8/2/2019 Dan Um Jay A
159/560
Fixed Point IteratIf it is known that a function g(x) has point, one way to approximate the va
fixed point is to use what is known aspoint iteration scheme.
8/2/2019 Dan Um Jay A
160/560
DefinitionA fixed point iteration scheme to apprfixed point p of a function g(x), genera
sequence {pn} by the rule pn = g(pnn 1, given a starting approximation
8/2/2019 Dan Um Jay A
161/560
Connection between fixed point problroot finding problems.
Every root finding problem can be trainto any number of different fixed poin
Some of these fixed point problems wrapidly, some will converge slowly an
not converge at all.
8/2/2019 Dan Um Jay A
162/560
The conversion process is actually qu
Take the root finding equation f(x) =algebraically transform it into an equaform
x = the expression on the right-hand sideresulting function is a corresponding
function g(x).
Example
8/2/2019 Dan Um Jay A
163/560
Example
Find a suitable interval and corresponiterative function g(x) such that the fix
iteration converges to the solution of t
f(x) = ex x 2 = 0.Perform four iterations.
Solution. We now find the suitable inte
8/2/2019 Dan Um Jay A
164/560
We find the function values at differenWe have
f(0) = 1 < 0,f(1) =
0.28 < 0,
f( 2 ) = 3.389 > 0.
Therefore, root lies on the interval I =
Rewrite the given equation as
8/2/2019 Dan Um Jay A
165/560
We check the conditions for the existeuniqueness of the fixed point method
g( 1 ) = 1.0986 > 1,
g( 2 ) = 1.3862 < 2.
This shows that
g : [1, 2]
[1, 2]
and also ( ) is continuous on [1 2]
8/2/2019 Dan Um Jay A
166/560
Now
g
(x) =
1
x + 2 ,and
maxxI |g
(x)| = g
(1) = 0.333 kTherefore, the function g(x) satisfies
properties of fixed point method.
8/2/2019 Dan Um Jay A
167/560
We denote fixed point iteration schem
xn = g(xn1), n = 1, 2, 3, We choose x0 = 1, we find
x1 = g(x0) = g(1) = 1.0986,x2 = g(x1) = g(1.0986) = 1.13
x3 = g(x2) = g(1.13095) = 1.1
x4 = g(x3) = g(1 14134) = 1 1
Example
8/2/2019 Dan Um Jay A
168/560
Example
Consider the function
f(x) = x3
+ x2
3x 3,has a unique zero on the interval [1, 2
We approximate this root using fixed iteration.
We start the equation
8/2/2019 Dan Um Jay A
169/560
We express
x = x3 + x2 33
= g1(x).
Similarly, we can also find
g2(x) = 1 + 3x + 3x2
,
g3(x) =
3 + 3x x21/3
8/2/2019 Dan Um Jay A
170/560
Now, choose suitable gi(x), i = 1, 2, 3
it has a unique fixed point and usin
point iteration find it.
8/2/2019 Dan Um Jay A
171/560
TheoremLet g(x) be continuous on the closed [a, b] with g : [a, b]
[a, b].
Furthermore, Suppose that g(x) is difon the open interval (a, b) and there e
positive constant k < 1 such that |g(xfor all x (a, b).
Then the sequence {pn} generated b( ) t fi d i
8/2/2019 Dan Um Jay A
172/560
Proof. We need to show
|pn p| 0 as n for any starting value p0 [a, b]. Takin
|pn p| = |g(pn1) g(p)= |g()| |pn1 p
k|pn1
p| k2|p p|
8/2/2019 Dan Um Jay A
173/560
Now, since k < 1, we get
limn |pnp| limn kn
|p0p| = |p0p
Order ofConvergence
8/2/2019 Dan Um Jay A
174/560
Convergence
Let g(x) be a continuous function on tinterval [a, b] with g : [a, b]
[a, b] and
that g(x) is continuous on the open inwith |g(x)| k < 1 for all x (a, b).
If g(p) = 0 then for any p0 [a, b], thepn = g(pn1) converges only linearly tpoint p.
8/2/2019 Dan Um Jay A
175/560
We need to prove if g(p) = 0, then thconvergence is only linear.
In other words, it must be shown that
limn |pn+1 p||pn p| = ,
where is a asymptotic error constan
8/2/2019 Dan Um Jay A
176/560
Note that pn p, we get cn p.
Since g(x) is continuous on (a, b), we
limn
|g(cn)| = |g( limn cn)| = |gHence,
limn
|pn+1 p|
|pn p|=
|g(p)
|.
Thi h h h d f
Theorem
8/2/2019 Dan Um Jay A
177/560
Theorem
If the iteration function g(x) is such thcontinuous in some neighborhood of
point p and g(p) = 0, g(p) = 0, then point method converges quadratically
Proof. We calculate
en+1 = ppn+1 = g(p) g(
= g(p) g(p en).
8/2/2019 Dan Um Jay A
178/560
Using Taylor series expansion, we ge
en+1 = g(p) g(p) eng(p) + e2n
2!
where n lies between pn and p. We g
en+1 = eng(p) e
2n
2g(n)
U i th f t ( ) 0 fi d
8/2/2019 Dan Um Jay A
179/560
en+1 =
e2n
2
[g(p) + n] ,
where n 0 as n .
We calculate |en+1|
|en
|2
=1
2|g(p) + n|.
Finally we get
Example
8/2/2019 Dan Um Jay A
180/560
Example
Use the fixed-point iteration to find a function
f(x) = ex sin x = 0.Solution. We write f(x) = 0 in the form
x = x + ex sin x = g(x)Since
f (0 5) 0 127 > 0
8/2/2019 Dan Um Jay A
181/560
To verify g(x) is a convergent iterationwe need to check the existence and u
conditions:
g(0.5 ) = 0.6271 > 0.5,
g(0.7 ) = 0.552367 < 0.7
we obtain
0 5 < g(x) < 0 7 for all x
8/2/2019 Dan Um Jay A
182/560
We find
g
(0.5) = 0.48,g(0.7) = 0.26.
Since g(x) is a monotone function on|g(x)| < 1, for all x I
Hence, the fixed point iteration will coi h I
Newtons MethodMultiple Roots
8/2/2019 Dan Um Jay A
183/560
Multiple Roots
A root p of f(x) is called a multiple rom, if
f(x) = (x p)m q(x),where limxp q(x) = 0, and m is a posinteger.
Theorem
8/2/2019 Dan Um Jay A
184/560
Theorem
Newtons method converges linearly multiple root.
But if pn+1 is taken as
pn+1 = pn m f(pn)f(pn) ,
where p is a multiple root of f(x) of o
it converges quadratically.
8/2/2019 Dan Um Jay A
185/560
Proof. We define
g(x) = x
f(x)f(x) x = p
p x = p,
where f(x) = (x p)m
q(x), p is multiporder m.
We prove continuity of g(x) at p
8/2/2019 Dan Um Jay A
186/560
Taking
limxp
g(x) = limxp
x
= limxp x limxp(x
p)mq(x)
m(xp)m1q(x) + (xThis implies that
lim g(x) = p = g(p)
8/2/2019 Dan Um Jay A
187/560
We now find g(p):
g(p) = limxp
g(x) g(p)x p .
Substituting the expression forg(x)
, w
g(p) = limxp
x f(x)f(x)
p
x p
8/2/2019 Dan Um Jay A
188/560
We now prove the continuity of g(x) a
We know that
g(x) = x
f(x)
f
(x)
,
when x = p. Then
g(x) = f (x) f
(x)[f ( )]2
8/2/2019 Dan Um Jay A
189/560
Now using the expression for f(x), w
limxp
g(x) = 1 1m
= g(p)
Therefore, g(x) is continuous at p.
We conclude that the convergence ra
8/2/2019 Dan Um Jay A
190/560
For obtaining the quadratic convergenthe modified Newtons iteration
pn+1 = pn m f(pn)f(pn)
.
To prove that this converges quadratisufficient to prove that
g
(p) = 0
Example
8/2/2019 Dan Um Jay A
191/560
p
Find the root of the equation
f(x) = x2
+ 2x + 1 = 0,
using Newtons method for simple roomultiple roots.
Solution. Given
f (x) = x2
+ 2x + 1
8/2/2019 Dan Um Jay A
192/560
The iteration formula for finding simplNewtons method is given by
xn+1 = xn f(xn)f(xn)
= xn x2n + 2xn + 12(xn + 1)
=
xn2
1
2 n = 0 1 2 3
8/2/2019 Dan Um Jay A
193/560
x1 =
0.5,
x2 = 0.75,x3 = 0.875,
x4 = 0.9375,...x20 =
0.999999,
x21 = 1
8/2/2019 Dan Um Jay A
194/560
Using Newtons iterative formula for mis given by
xn+1 = xn m
f(xn)
f(xn)
,
= xn 2x2n + 2xn +
2(xn + 1)
= xn (xn + 1)1
8/2/2019 Dan Um Jay A
195/560
Here, we obtain the root in the first ite
Hence, Newtons method for multiple much faster than Newtons method foroot.
Example
8/2/2019 Dan Um Jay A
196/560
p
Do three iterations of Newtons methothe double root of
x3 2x2 0.75x + 2.25 = 0which is close to 1 such that iterationsquadratically.
Solution. We use modified Newtons m
f(p )
8/2/2019 Dan Um Jay A
197/560
We have
f(pn) = p3
n 2p2
n 0.75pn +f(pn) = 3p
2n 4pn 0.75
Thenpn+1 =
p3n + 0.75pn 4.53p2n 4pn 0.75
.
8/2/2019 Dan Um Jay A
198/560
We calculate
n = 0, p1 = 1.5714n = 1, p2 = 1.5009
n = 2, p3 = 1.5
System of NonlinEquations
8/2/2019 Dan Um Jay A
199/560
q
We consider the system of n non-lineequations
f1(x1, x2, x3, , xn) = 0f2(x1, x2, x3, , xn) = 0
...
fn(x1, x2, x3, , xn) = 0
where x1 x2 xn are n unknowns
Newtons Method
8/2/2019 Dan Um Jay A
200/560
Newton s Method
We discuss the Newtons method for two non-linear equations:
f(x, y) = 0,
g(x, y) = 0.
Let (x0, y0) be an initial approximationof the system.
8/2/2019 Dan Um Jay A
201/560
If (x0 + h, y0 + k) is the root of the syswe must have
f(x0 + h, y0 + k) = 0,
g(x0 + h, y0 + k) = 0.
Assuming that f and g are sufficientlydifferentiable, we expand by Taylors sobtain
f f
8/2/2019 Dan Um Jay A
202/560
g0
+ hg
x0+ k
g
y0+ higher order t
where f0 = f(x0, y0) andfx0
=fxx=x
Neglecting higher order terms, we obsystem of equations
h f k f f
8/2/2019 Dan Um Jay A
203/560
We now define the Jacobian
J(f, g) =
f
x
f
ygx
gy
If the Jacobian does not vanish, then unique solution given by
h =gfy
f gyJ(f g)
8/2/2019 Dan Um Jay A
204/560
The new approximations are given by
x1 = x0 + h,y1 = y0 + k.
The process is to be repeated till we oroots to the desired accuracy.
8/2/2019 Dan Um Jay A
205/560
We now write the Newtons iteration f
xn+1 = xn f gy gfyJ(f, g)(xn,y
andyn+1 = yn
gfx f gx
J(f, g)
(xn,y
where n = 0 1 2 3
Example
8/2/2019 Dan Um Jay A
206/560
Solve the system of nonlinear equatio
x2
y3
x + y 1.1 2 5 = x3 + y2 + x y + 0.8 7 5 =
by Newtons method starting withx0 = 0.2, y0 = 0.2.
Fixed Point Meth
8/2/2019 Dan Um Jay A
207/560
We consider the system of equations
F(x, y) = 0,G(x, y) = 0.(2)
We assume that the system (2) have (, ).
We rewrite the above system (2)
( )
http://-/?-http://-/?-http://-/?-http://-/?-8/2/2019 Dan Um Jay A
208/560
Let (xn, yn) be nth approximations tocomputed from the iteration scheme
xn+1 = f(xn, yn),
yn+1 = g(xn, yn), n = 0, 1, 2(4)
Note that
= f(, ),(5)
= g( )(6)
8/2/2019 Dan Um Jay A
209/560
Now substracting (4) from (6), we get
xn+1 = f(, ) f(xn, yn)= f(xn + xn, yn +
f(xn, yn).
http://-/?-http://-/?-http://-/?-http://-/?-8/2/2019 Dan Um Jay A
210/560
Using the Taylors series expansion, w
xn+1 = f(xn, yn) + ( xn)fx(xn+ ( yn)fy(xn, yn) + high
f(xn, yn)
= ( xn)fx(xn, yn) + ( + higher order terms.
Similarly we find
8/2/2019 Dan Um Jay A
211/560
Neglecting the higher order terms, wethe matrix form
xn+1 yn+1
=
fx(xn, yn) fy(xn, yn)
gx(xn, yn) gy(xn, yn)
In compact form
en+1 = Ben,
h B i th ffi i t t i
Theorem
8/2/2019 Dan Um Jay A
212/560
If fx and fy are continuous in R (rectaregion, containing (, ) and (xn, yn) f
n = 0, 1, 2, ) then it is proved that thscheme is given by
xn+1 = f(xn, yn),yn+1 = g(xn, yn), n = 0, 1, 2
converges to (, ) if
8/2/2019 Dan Um Jay A
213/560
B k < 1,for some number k, i.e.,
|fx| + |fy| k < 1,
and|gx| + |gy| k < 1.
Example
8/2/2019 Dan Um Jay A
214/560
Consider the nonlinear system
x = 1 + sin
8 (x + y)
,
y = 1
cos
8
(x
y)
8/2/2019 Dan Um Jay A
215/560
(i) Find a region D 2 such that the
xn+1 = 1 + sin
8
(xn + yn)
yn+1 = 1
cos
8
(xn
yn)
is guaranteed to converge to a uniquefor any (x0, y0)
D.
8/2/2019 Dan Um Jay A
216/560
(ii) Determine a bound on the rate of convergence of the method in max no
Solution. We have
xy
= H(x, y) = f(x, y)
g(x, y)
= 1 + sin
8 (x + y)
1 cos 8 (x y)
8/2/2019 Dan Um Jay A
217/560
Notice that for any (x, y), we haveH(x, y)
[0, 2]
[0, 2].
Thus, if R is any closed and boundedsuch that [0, 2] [0, 2] R, then H(D
8/2/2019 Dan Um Jay A
218/560
We define
B = fx(x, y) fy(x, y)
gx(x, y) gy(x, y)
=
8
cos 8
(x + y) 8
cos 8
(8 sin
8 (x y)
8 sin8
8/2/2019 Dan Um Jay A
219/560
Note that
B(x, y)
=
max{4 | cos 8 (x + y) |, 4 | sin 8 (x y1, for any (x, y) D.
Example
8/2/2019 Dan Um Jay A
220/560
Consider the nonlinear system
x = 34 ln
1 + 12(x + y)2
y =3
4
ln1 +1
2
(x
y)2
(i) Find a region D 2, such that th
x 1 = 3 ln
1 + 1(x + y )
8/2/2019 Dan Um Jay A
221/560
is guaranteed to converge to the root(x, y) = (0, 0) for any (x, y)
D.
Solution. Plz try.
8/2/2019 Dan Um Jay A
222/560
C H A PT ER System of Line
Equations
8/2/2019 Dan Um Jay A
223/560
Contents
1. Gauss Elimination Method2. Inverse of a Matrix
3. Condition Numbers and Errors
ill-Conditioned System
4. Iterative methods
System of Linear Equations
8/2/2019 Dan Um Jay A
224/560
Circuit analysis (Mesh and node equatio
Numerical solution of differential equation
Difference Method).
Numerical solution of integral equation
Element Method, Method of Moments).
nnnnnn
nn
nn
bxaxaxa
bxaxaxabxaxaxa
=+++
=+++=+++
2211
22222121
11212111 11 12 1
21 22 2
1 2
n
n
n n nn
a a a a a a
a a a
Notations
8/2/2019 Dan Um Jay A
225/560
The above system of equations can be w
matrix from as Ax = b.
A is called the coefficient matrix, A = (aij
b is called the constant vector= (b1,b2,
[A|b] is called the augmented matrix. If , then a multiplier is defined as
1
1 2 3 4i
i
a
m i
011 a
Consistency (Solvability)
8/2/2019 Dan Um Jay A
226/560
The linear system of equations Ax=bsolution, or said to be consistent if and on
Rank{A}=Rank{A|b}
A system is inconsistent when
Rank{A}
8/2/2019 Dan Um Jay A
227/560
The following operations applied to the au
matrix [A|b], yields an equivalent linear sys
Interchanges: Any two rows interchanged.
Scaling: Multiplying a row by a
constant.
Replacement: The row can be replace
sum of that row and a nonzero multip
An inconsistent example
8/2/2019 Dan Um Jay A
228/560
=
5
4
42
21
2
1
x
x
00
21Rank{A}=1
R k{A|b} 2
ERO:Multiply the first row with-2 and add to the second row
421
The systemof equations
are not
Uniqueness of solutions
8/2/2019 Dan Um Jay A
229/560
The system has a unique solution if a
only if Rank{A} = Rank{A|b} = n,
n is the order of the system.
Such systems are called full-rank
systems.
Full-rank systems Example
8/2/2019 Dan Um Jay A
230/560
If Rank{A}=n
Det{A} 0 A is nonsingular so inverand has a Unique solution.
=
2
4
11
21
2
1
x
x
Rank deficient matrices
8/2/2019 Dan Um Jay A
231/560
If Rank{A}=m
8/2/2019 Dan Um Jay A
232/560
A small change in the entries of matri
causes a large deviation in the solutio
1
2
1 2 3
0.48 0.99 1.47
x
x
=
1
2
1 2 3
0 49 0 99 1 47
x
x
=
=
1
1
2
1
x
x
=
0
3
2
1
x
x
Ill-conditioned continued.....
8/2/2019 Dan Um Jay A
233/560
A linear system of
equations is said to
be ill-conditioned
if the coefficient
matrix tends to be
singular
Solution Techniques
8/2/2019 Dan Um Jay A
234/560
Direct solution methods
Finds a solution in a finite number of opetransforming the system into an equivalent sy
is easier to solve. Diagonal, upper or lower triangular systems
to solve.
Number of operations is a function of system
Iterative solution methods
Computes succesive approximations of thvector for a given A and b, starting from an i
x0
Direct solution Methods
8/2/2019 Dan Um Jay A
235/560
Gaussian Elimination By using ERO, matrix A is transformed
upper triangular matrix (all element
diagonal are 0).
Back substitution is used to solve th
triangular system.
=
iiiniii
ni
b
b
x
x
aaa
aaa
11
1
1111
ERO
iinii
ni
x
x
aa
aaa
~~0
11111
Back Substitution
8/2/2019 Dan Um Jay A
236/560
After the forward elimination phase, the m
been transformed into upper triangularfo
Equation njust involves xn and so we ceasily.
Equation n-1 just involves xn-1 and xn, and
already knowxn we can findxn-1.
Working our way backwards through the e
we can findxn
, xn-1
,, x1
.
First step of elimination
8/2/2019 Dan Um Jay A
237/560
=
=
=
2()2(3
)2(2
2(
3
)2(
33
)2(
32
2(
2
)2(
23
)2(
22
1(
1
)1(
13
)1(
12
)1(
11
)1(11
)1(11
)1(
11
)1(
311,3
)1(
11
)1(
211,2
0
0
0
/
/
/
n
n
n
nn aaa
aaa
aaa
aaaa
aam
aam
aam
)1()1(
3
)1(
2
)1(
1
)1(
3
)1(
33
)1(
32
)1(
31
)1(
2
)1(
23
)1(
22
)1(
21
)1(
1
)1(
13
)1(
12
)1(
11
nnnnn
n
n
n
aaaa
aaaa
aaaa
aaaa
Pivotal element
Second step of elimination
8/2/2019 Dan Um Jay A
238/560
= )3(3
)3(
33
)2(
2
)2(
23
)2(
22
)1(
1
)1(
13
)1(
12
)1(
11
)2()2(
)2(
22
)2(
322,3 00
0
/n
n
n
aa
aaaaaaa
aam
)2()2(
3
)2(
2
)2(
3
)2(
33
)2(
32
)2(
2
)2(
23
)2(
22
)1(
1
)1(
13
)1(
12
)1(
11
0
0
0
nnnn
n
n
n
aaa
aaa
aaa
aaaa
Pivotal element
8/2/2019 Dan Um Jay A
239/560
Like this (n-1) steps would be done so
final matrix would be such that upper tr
gives the solutions and the numbers
diagonal entries are the multipliers.
Back substitution algorithm
8/2/2019 Dan Um Jay A
240/560
=
1
3
2
1
)(
)1(
1
)1(
11
)3(
3
)3(
33
)2(
2
)2(
23
)2(
22
)1(
1
)1(
13
)1(
12
)1(
11
0000
000
00
0
n
n
n
nn
n
nn
n
nn
n
n
n
b
x
x
x
x
x
a
aa
aa
aaa
aaaa
[ ]
1211
1
)()(
1
1
)1(
1)1(
11
1)(
)(
==
nnixabx
xaba
xa
bx
n
ii
n
n
nn
n
nn
nn
nn
nn
n
nn
Pivoting
8/2/2019 Dan Um Jay A
241/560
Computer uses finite-precisionarithmet
A small error is introduced in each a
operation, error propagates. When the pivotal element is very sm
multipliers will be large.
Adding numbers of widely differring magnlead to loss of significance.
To reduce error, row or column interchanges are made to maximizethe m
of the pivot element
Simple pivoting
8/2/2019 Dan Um Jay A
242/560
Suppose at the first step then this numnot be used to eliminate other numbers ofcolumns of other rows.
To avoid this in first step if use this aelement at first step.
Similarly if at second step, chose a
number from 2nd
column and from 2nd
row onw Avoiding zero for pivotal element is kn
simple pivoting.
In Simple pivoting also, due to finite digit aresults are in large errors
,011 =a
,021 a
,022 =a
Example: Without Pivoting
8/2/2019 Dan Um Jay A
243/560
=
22
.6
210.114.24
281.5133.1
2
1
x
x
=
7.113000.0
281.5133.1
2
1
x
x
1
2
0.9956
1 001
x
x
=
31.21
133.1
14.2421
==m
4-digit arithmetic
Loss of signifi
Example: With Partial Pivoting
8/2/2019 Dan Um Jay A
244/560
=
414.6
93.22
281.5133.1
210.114.24
2
1
x
x
=
338.5000.0
210.114.24
2
1
x
x
000.1
1
x
04693.014.24
133.121 ==m
Partial Pivoting
8/2/2019 Dan Um Jay A
245/560
In some cases, akk(k) is very close to zero. W
divide it with a bigger number, roundin
arises which causes the damage of exact s
Using such an element as the pivot ele
result in gross errors in the further calcu
the matrix.
To check it, we introduce partial pivo
complete pivoting.
Partial Pivoting Strategy
8/2/2019 Dan Um Jay A
246/560
For in the GE process at stage
Let j be the row index at which the maxim
attained.
Ifj > k, then interchange k-th andj-th row in Now, and this h
preventing the growth of elements in A(k) o
varying size and causes to reduce the
,11 nk
.|,|)(
nikaMaxck
ikk =
| | 1, 1,...,ikm i k n = +
Pivoting procedures
8/2/2019 Dan Um Jay A
247/560
)()()(
)()()(
)()()(
3(3
)3(3
)3(3
)3(33
2(
2
)2(
2
)2(
2
)2(
23
)2(
22
1(
1
)1(
1
)1(
1
)1(
13
)1(
12
)1(
11
000
000
000
00
0
iii
i
jn
i
jj
i
ji
i
in
i
ij
i
ii
nji
nji
nji
aaa
aaa
aaaa
aaaaa
aaaaaa
Eliminatedpart
Partial Row pivoting
8/2/2019 Dan Um Jay A
248/560
)()()(
)()()(
)()()(
)3(3
)3(3
)3(3
)3(33
)2(
2
)2(
2
)2(
2
)2(
23
)2(
22
)1(
1
)1(
1
)1(
1
)1(
13
)1(
12
)1(
11
000
000
00
0
iii
i
jn
i
jj
i
ji
i
in
i
ij
i
ii
nji
nji
nji
aaa
aaa
aaaa
aaaaa
aaaaaa
Inth
Partial Column pivoting
8/2/2019 Dan Um Jay A
249/560
)()()(
)()()(
)()()(
)3(3
)3(3
)3(3
)3(33
)2(
2
)2(
2
)2(
2
)2(
23
)2(
22
)1(
1
)1(
1
)1(
1
)1(
13
)1(
12
)1(
11
000
000
00
0
iii
i
jn
i
jj
i
ji
i
in
i
ij
i
ii
nji
nji
nji
aaa
aaa
aaaa
aaaaa
aaaaaa
Complete pivoting
8/2/2019 Dan Um Jay A
250/560
)()()(
)()()(
)()()(
)3(3
)3(3
)3(3
)3(33
)2(
2
)2(
2
)2(
2
)2(
23
)2(
22
)1(
1
)1(
1
)1(
1
)1(
13
)1(
12
)1(
11
000
000
000
00
0
iii
i
jn
i
jj
i
ji
i
in
i
ij
i
ii
nji
nji
nji
aaa
aaa
aaaa
aaaaa
aaaaaa
Inth
Remarks
8/2/2019 Dan Um Jay A
251/560
There is a possible change in order of the
in complete pivoting. So, the order of unkno
be reversed after back substitution in GE m In practical problems, the error beh
essentially the same in both partial piv
complete pivoting. The complete pivoting take more operatio
partial. So, we prefer the latter in mos
algorithms
8/2/2019 Dan Um Jay A
252/560
Chapter-3
System of Linear Equations
Dr. P. Dhanumjaya
Example 1
8/2/2019 Dan Um Jay A
253/560
Consider the following system
0.003000x1 + 59.14x2 = 59.5.291x1 6.130x2 = 46.
has the exact solution x1 = 10.00 and
We now find the solution using Gaussmethod and four-digit arithmetic with
8/2/2019 Dan Um Jay A
254/560
We form an augmented matrix, we ob
0.003 59.14... 59.17
5.291 6.130... 46.78
The multipliers are
mi1 =ai1
a11, i = 2, 3, 4,
8/2/2019 Dan Um Jay A
255/560
We now replace
aij =
aij
mi1a1
j, i, j
= 2,
3and
ai1 = 0, i = 2, 3, 4,
We have
0.003 59.14
... 59.17
(1764) 104300 10440
8/2/2019 Dan Um Jay A
256/560
Using back substitution, we get
x2 1
.001
, and x1 1
This shows that a small value of pivotin the numerator can be dramatically
Partial Pivoting
8/2/2019 Dan Um Jay A
257/560
The augmented matrix is
0.003 59.14.
.. 59.175.291 6.130
... 46.78
We takemax{|a11|, |a21|} = 5.291.
We replace r1 by r2 we get
8/2/2019 Dan Um Jay A
258/560
The multiplier is
m21 =
a21
a11 = 0.000567.
We replace
a22 = a22 m21a12,
we get
8/2/2019 Dan Um Jay A
259/560
This gives
x2= 1
.0
, and x1= 10
.0
Example 2
8/2/2019 Dan Um Jay A
260/560
Consider the following system
30.00
x1+ 591400
x2= 59175.291x1 6.130x2 = 46.78
The augmented matrix is
30.0 591400... 591700
5.291 6.130... 46.78
8/2/2019 Dan Um Jay A
261/560
The maximum value in the first colum
The multiplier is
m21 =a21
a11=
5.291
30.0= 0.176
We replace
a22 = a22 m21a12.
8/2/2019 Dan Um Jay A
262/560
We get
30.0 591400.
.. 59170(0.1764) 104300
... 10440
Using back substitution, we getx2 1.001, and x1 10
This shows that in a row if the coeffic
Scaled PartialPivoting
8/2/2019 Dan Um Jay A
263/560
Scaling
di = max1jn |
aij|
, i= 1
,2
,3
,
then divide ith row by di for i = 1, 2, 3
Choose numerically largest element oconcerned column as pivotal element
This is known as scaled partial pivotin
Note 1
8/2/2019 Dan Um Jay A
264/560
Note that
scaling is useful only when size of
coefficients vary too much.
Example 3
8/2/2019 Dan Um Jay A
265/560
Using Gaussian Elimination with partscaling and back substitution, solve th
algebraic system
6.684x1 + 2.925x2 + 9.835x3 =
5.543
x1+ 5
.953
x2+ 88
.63
x3= 8.375x1 + 2.988x2 + 8.681x3 =
Solution
8/2/2019 Dan Um Jay A
266/560
The augmented matrix is
6.684 2.925 9.835
..
. 10.755.543 5.953 88.63
... 19.81
8.375 2.988 8.681... 24.72
The scaling factors are
d19
.835
, d288
.63
, d30
8/2/2019 Dan Um Jay A
267/560
After scaled partial pivoting, the resulis
0.96475 0.34419 1... 2.847
0.06254 0.06717 1... 0.223
0.67961 0.29741 1... 1.093
The multipliers are
21
a21
0 06483
8/2/2019 Dan Um Jay A
268/560
We replace a22, a32 by
a22 = a22 m21a12,
a32 = a32 m31a12.
The resulting matrix is
0.96475 0.34419 1... 2.
0 0.06312 0.93517... 0.
8/2/2019 Dan Um Jay A
269/560
The multiplier is
m32 =
a32
a22 = 0.87056,
and replace a33 by
a33 = a33 m32a23.
The resulting matrix is
8/2/2019 Dan Um Jay A
270/560
Using the back substitution, we find
x3 = 1.82586,
x2 = 26.4352,
x1 = 10.4903.
Example 4
8/2/2019 Dan Um Jay A
271/560
Using Gaussian Elimination with partscaling and back substitution, solve th
algebraic system
0.995x1 + 1.54x2 + 4.51x3 =
0.995
x1+ 2
.16
x2+ 1
.19
x3= 0.298x1 + 0.577x2 + 1.42x3 =
Solution. Plz try.
Crouts Method
8/2/2019 Dan Um Jay A
272/560
Given any
AX = b,
we express the matrix A as a decomplower triangular matrix L and an uppematrix U such that
L U = A.
The different factorizations may be vie
resulting from different choices for the
8/2/2019 Dan Um Jay A
273/560
The two most common choices for thentries are
lii = 1, for each i = 1, 2, 3,
uii = 1, for each i = 1, 2, 3,
If lii = 1, for each i = 1, 2, 3, , n themethod is known as the Doolittle decomethod.
Crouts Method
8/2/2019 Dan Um Jay A
274/560
Let A be an n n matrix. To obtain thdecomposition of A, we must determ
entries lij (i j) and uij (i < j) such
8/2/2019 Dan Um Jay A
275/560
l11
l21 l22
l31 l32 l33. . .
.... . .
ln1 ln2 ln3... lnn
1 u12 u13
1 u23
1. .
.
a11 a12 a13
8/2/2019 Dan Um Jay A
276/560
The product of the ith row of L (fori = 1, 2, 3, , n) with the first column
simply the element li1.
The value is equated to ai1, that is
li1 = ai1, i = 1, 2, 3, , n.
These equations determine the first c
8/2/2019 Dan Um Jay A
277/560
Now multiply the first row of L with thcolumn of U (for j = 2, 3, 4, , n) and
to a1j produces the equation
l11u1j = a1j,
this impliesu1j =
a1j
l11.
8/2/2019 Dan Um Jay A
278/560
This determining the first row of U.
Similarly, we compute one more coluone more row of U.
In particular, the kth column of L andof U are determined as
lik = aik k1
j 1
lijujk, i = k, k + 1, k
O h ffi i i h b
8/2/2019 Dan Um Jay A
279/560
Once the coefficient matrix has been to its L U decomposition, then
AX = b L U X = b.
Let U X = b1 then Lb1 = b. This gives
l11
l21 l22
l31 l32 l33
b11b12
Thi i
8/2/2019 Dan Um Jay A
280/560
This gives
b11 =b1
l11,
and
b1i = 1lii
bi i1k=1
likb1k , i = 2, 3, 4,
Now using back substitution for
t X
8/2/2019 Dan Um Jay A
281/560
we get X as
xn
= b1n
,
xj = b1
j n
k=j+1
ujkxk, j = n 1, n
Example 5
S l th f ll i t f ti
8/2/2019 Dan Um Jay A
282/560
Solve the following system of equatioCrouts method
x1 + 2x2 + x3 = 2,
2x1 + x2 10x3 = 4,
2x1 + 3x2 x3 = 2.Here, the coefficient matrix is
1 2 1
2
W A i t f ll
8/2/2019 Dan Um Jay A
283/560
We express A into as follows
l11 0 0l21 l22 0
l31 l32 l33
1 u12 u130 1 u230 0 1
=
1 2
2
Thi i
8/2/2019 Dan Um Jay A
284/560
This gives
l11 = 1, l21 = 2, l31 = 2,
and
u12 =2
l11= 2, u13 = 1.
Similarly, we find
l22 3 l32 1 l33 1
Th f
8/2/2019 Dan Um Jay A
285/560
Therefore,
L =
1 0 0
2 3 02 1 1
,
and
U =
1 2 1
0 1 4
0 0 1
.
W d t t
8/2/2019 Dan Um Jay A
286/560
We need to compute
1 0 02 3 0
2 1 1
b
1
1
b12b13
=
24
2
This gives
b11 = 2, b12 = 0, b
13 = 2.
Finall e com te
8/2/2019 Dan Um Jay A
287/560
Finally, we compute
1 2 10 1 4
0 0 1
x1x2
x3
=
20
2
This gives
x3 = 2, x2 = 8, x1 = 1
8/2/2019 Dan Um Jay A
288/560
Chapter-3
System of Linear Equations
Dr. P. Dhanumjaya
Example 1
Consider the following system
8/2/2019 Dan Um Jay A
289/560
Consider the following system
0.003000x1 + 59.14x2 = 59.
5.291x1 6.130x2 = 46.has the exact solution x1 = 10.00 and
We now find the solution using Gaussmethod and four-digit arithmetic with
We form an augmented matrix we ob
8/2/2019 Dan Um Jay A
290/560
We form an augmented matrix, we ob
0.003 59.14... 59.17
5.291 6.130 ... 46.78
The multipliers are
mi1 =ai1
a11, i = 2, 3, 4,
We now replace
8/2/2019 Dan Um Jay A
291/560
We now replace
aij = aij
mi1a1j, i, j = 2, 3
and
ai1 = 0, i = 2, 3, 4,
We have
0.003 59.14... 59.17
(1764) 104300 10440
Using back substitution we get
8/2/2019 Dan Um Jay A
292/560
Using back substitution, we get
x2
1.001, and x1
1
This shows that a small value of pivotin the numerator can be dramatically
Partial Pivoting
The augmented matrix is
8/2/2019 Dan Um Jay A
293/560
The augmented matrix is
0.003 59.14... 59.17
5.291 6.130 ... 46.78
We takemax{|a11|, |a21|} = 5.291.
W l 1 b 2 t
The multiplier is
8/2/2019 Dan Um Jay A
294/560
The multiplier is
m21 =
a21
a11 = 0.000567.
We replace
a22 = a22 m21a12,we get
This gives
8/2/2019 Dan Um Jay A
295/560
This gives
x2 = 1.0, and x1 = 10.0
Example 2
Consider the following system
8/2/2019 Dan Um Jay A
296/560
Consider the following system
30.00x1 + 591400x2 = 5917
5.291x1 6.130x2 = 46.78The augmented matrix is
30.0 591400... 591700
5.291 6.130... 46.78
The maximum value in the first colum
8/2/2019 Dan Um Jay A
297/560
The maximum value in the first colum
The multiplier is
m21 =a21
a11=
5.291
30.0= 0.176
We replace
a22 = a22 m21a12.
We get
8/2/2019 Dan Um Jay A
298/560
We get
30.0 591400... 59170
(0.1764) 104300 ... 10440
Using back substitution, we getx2 1.001, and x1 10
Thi h th t i if th ffi
Scaled PartialPivoting
Scaling
8/2/2019 Dan Um Jay A
299/560
Scaling
di = max1jn |
aij|, i = 1, 2, 3,
then divide ith row by di for i = 1, 2, 3
Choose numerically largest element oconcerned column as pivotal element
Thi i k l d ti l i ti
Note 1
Note that
8/2/2019 Dan Um Jay A
300/560
Note that
scaling is useful only when size of
coefficients vary too much.
Example 3
Using Gaussian Elimination with part
8/2/2019 Dan Um Jay A
301/560
Using Gaussian Elimination with partscaling and back substitution, solve th
algebraic system6.684x1 + 2.925x2 + 9.835x3 =
5.543x1 + 5.953x2 + 88.63x3 =
8.375x1 + 2.988x2 + 8.681x3 =
Solution
The augmented matrix is
8/2/2019 Dan Um Jay A
302/560
The augmented matrix is
6.684 2.925 9.835
..
. 10.755.543 5.953 88.63
... 19.81
8.375 2.988 8.681... 24.72
The scaling factors are
d1 = 9.835, d2 = 88.63, d3 = 0
After scaled partial pivoting the resul
8/2/2019 Dan Um Jay A
303/560
After scaled partial pivoting, the resulis
0.96475 0.34419 1... 2.847
0.06254 0.06717 1... 0.223
0.67961 0.29741 1... 1.093
The multipliers are
21
a21
0 06483
We replace a22 a32 by
8/2/2019 Dan Um Jay A
304/560
We replace a22, a32 by
a22 = a22
m21a12,
a32 = a32 m31a12.The resulting matrix is
0.96475 0.34419 1... 2.
0 0.06312 0.93517... 0.
The multiplier is
8/2/2019 Dan Um Jay A
305/560
The multiplier is
m32 =
a32
a22 = 0.87056,
and replace a33 by
a33 = a33 m32a23.The resulting matrix is
Using the back substitution we find
8/2/2019 Dan Um Jay A
306/560
Using the back substitution, we find
x3 = 1.82586,
x2 = 26.4352,x1 = 10.4903.
Example 4
Using Gaussian Elimination with part
8/2/2019 Dan Um Jay A
307/560
Using Gaussian Elimination with partscaling and back substitution, solve th
algebraic system0.995x1 + 1.54x2 + 4.51x3 =
0.995x1 + 2.16x2 + 1.19x3 =
0.298x1 + 0.577x2 + 1.42x3 =
Solution. Plz try.
Crouts Method
Given any
8/2/2019 Dan Um Jay A
308/560
Given any
AX = b,
we express the matrix A as a decomplower triangular matrix L and an uppematrix U such that
L U = A.
The different factorizations may be vie
l i f diff h i f h
The two most common choices for th
8/2/2019 Dan Um Jay A
309/560
The two most common choices for thentries are
lii = 1, for each i = 1, 2, 3, uii = 1, for each i = 1, 2, 3,
If lii = 1, for each i = 1, 2, 3, , n themethod is known as the Doolittle decomethod.
Crouts Method
Let A be an n n matrix. To obtain th
8/2/2019 Dan Um Jay A
310/560
Let A be an n n matrix. To obtain thdecomposition of A, we must determ
entries lij (i j) and uij (i < j) such
8/2/2019 Dan Um Jay A
311/560
l11
l21 l22
l31 l32 l33. . .
.... . .
ln1 ln2 ln3... lnn
1 u12 u13
1 u23
1. . .
a11 a12 a13
The product of the ith row of L (for
8/2/2019 Dan Um Jay A
312/560
The product of the ith row of L (fori = 1, 2, 3, , n) with the first column simply the element li1.
The value is equated to ai1, that is
li1 = ai1, i = 1, 2, 3, , n.These equations determine the first c
Now multiply the first row of L with th
8/2/2019 Dan Um Jay A
313/560
Now multiply the first row of L with thcolumn of U (for j = 2, 3, 4, , n) andto a1j produces the equation
l11u1j = a1j ,
this impliesu1j =
a1j
l11.
This determining the first row of U .
8/2/2019 Dan Um Jay A
314/560
This determining the first row of U.
Similarly, we compute one more coluone more row of U.
In particular, the kth column of L and
of U are determined as
lik = aikk1
j
lij ujk, i = k, k + 1, k
Once the coefficient matrix has been
8/2/2019 Dan Um Jay A
315/560
O ce t e coe c e t at as beeto its L U decomposition, then
AX = b L U X = b.Let U X = b1 then Lb1 = b. This gives
l11
l21 l22
l l l
b11b1
This gives
8/2/2019 Dan Um Jay A
316/560
g
b11 =b1
l11,
and
b1i =1
lii
bi
i1k=1
likb1k
, i = 2, 3, 4,
Now using back substitution for
8/2/2019 Dan Um Jay A
317/560
Example 5
Solve the following system of equatio
8/2/2019 Dan Um Jay A
318/560
g y qCrouts method
x1 + 2x2 + x3 = 2,
2x1 + x2 10x3 = 4,
2x1 + 3
x2
x3 = 2
.
Here, the coefficient matrix is
We express A into as follows
8/2/2019 Dan Um Jay A
319/560
p
l11
0 0
l21 l22 0
l31 l32 l33
1 u12
u13
0 1 u230 0 1
=
1
2
2
This gives
8/2/2019 Dan Um Jay A
320/560
g
l11 = 1, l21 = 2, l31 = 2,
and
u12 =2
l11= 2, u13 = 1.
Similarly, we find
Therefore,
8/2/2019 Dan Um Jay A
321/560
,
L =
1 0 0
2 3 02 1 1
,
and
U =
1 2 1
0 1 4
.
We need to compute
8/2/2019 Dan Um Jay A
322/560
p
1 0 0
2 3 02 1 1
b1
1b12b13
=
2
4
2
This gives
b11 = 2, b12 = 0, b
13 = 2.
Finally, we compute
8/2/2019 Dan Um Jay A
323/560
y p
1 2 1
0 1 4
0 0 1
x1
x2
x3
=
2
0
2
This gives
x3 = 2, x2 = 8, x1 = 1
Vector and MatrixNorms
To measure the errors introduced dur
8/2/2019 Dan Um Jay A
324/560
solution of a linear system, it will be n
have a means for quantifying the size
This is done using matrix norms.
Vector Norms
The function : n is called a
8/2/2019 Dan Um Jay A
325/560
norm if for all X, Y n and all following properties hold:X 0,
X
= 0, iff X = 0,
X = || X andX + Y X + Y.
There are infinitely many vector norm
8/2/2019 Dan Um Jay A
326/560
be constructed.
We will restrict our attention to the coused norms:
1 (l1 norm), 2 (l2 norm,
Let X n, the l1-norm of X denoted
8/2/2019 Dan Um Jay A
327/560
defined by
X1 =n
i=1|xi|.
The l2-norm or Euclidean norm of X
by X2, is defined by
X =
n
x21
2
.
The l-norm or maximum norm is de
8/2/2019 Dan Um Jay A
328/560
X is defined byX = max
1in|xi|.
Example
Consider the three vectors
8/2/2019 Dan Um Jay A
329/560
X1 = [1
2 3]T,
X2 = [ 2 0 1 2]T,X3 = [ 0 1 4 2 1]T
The maximum norm of each of these computed as follows:
The l2-norm or Euclidean norm of eac
8/2/2019 Dan Um Jay A
330/560
X1
2 = 12 + (
2)2 + 32 =
14
X22 =
22 + 02 + (1)2 + 22 =
X3
2 = 0
2 + 12 + (
4)2 + 22 +
= 22 4.69.
Matrix Norms
A matrix norm is a function : n
8/2/2019 Dan Um Jay A
331/560
for all A, B nn and all satisA 0,A = 0, iff A = 0,
A = || A,A + B A + B andAB A B .
Let A be an n n matrix then
8/2/2019 Dan Um Jay A
332/560
The l1-norm denoted by
A
1 and
A1 = max1jn
ni=1
|aij|.
The l-norm or maximum norm dA and defined by
n
The l2-norm or Spectral norm den
8/2/2019 Dan Um Jay A
333/560
A2 and defined by
A2 =
(AT A),
where is the spectral radius and
(A) = maxi
|i|.
Example
We calculate both the l2 and the l n
8/2/2019 Dan Um Jay A
334/560
matrices
A1 = 1 2
4 3
and
A2 =
1 0 2
0 1 1 .
A1 = max{|1| + | 2| |4| + |3
8/2/2019 Dan Um Jay A
335/560
A1 max{|1| + | 2|, |4| + |3and
A2 = max{|1| + |0| + |2|, |0| + |1|
| 1
|+
|1
|+
To determine the l2-norm, we first com
1 4 1 2 1
The eigenvalues are
8/2/2019 Dan Um Jay A
336/560
=
1
2
31 517 .Hence,
(AT1 A1) =12
31 + 517
and
Similarly, we compute
8/2/2019 Dan Um Jay A
337/560
(AT2 A2) = 6.24914,
and
A2
2 =
6.24914
2.4998
Error Estimates &Condition Numbe
We now address the following questio
8/2/2019 Dan Um Jay A
338/560
1. Howmuch error can we expect wh
system of linear equations using Gelimination.
2. How does the error depend on the
of the coefficient matrix and the rigside vector?.
3. When A and b are known only app
Error estimates
Suppose X is an approximate solutio
8/2/2019 Dan Um Jay A
339/560
linear system
AX = b,whose exact solution is the vector X.
In practice, the exact solution to the sunknown, so the error in X
However, the residual vector
8/2/2019 Dan Um Jay A
340/560
r = AX
b,
can be easily computed.
Note
The residual measures the amount bi t l ti f il t ti f th
8/2/2019 Dan Um Jay A
341/560
approximate solution fails to satisfy thsystem.
When r = 0, this implies X is the exa
8/2/2019 Dan Um Jay A
342/560
so
e = 0.It seems reasonable to expect that whis small,
e
will be small as well.
Unfortunately, this need not always b
Example
Consider the linear system
8/2/2019 Dan Um Jay A
343/560
1 2
0.99 1.99 x1
x2
= 1
has X = [1 1]T as its exact solution.
Let the vector X = [1 0]T.
8/2/2019 Dan Um Jay A
344/560
The error e = [2 1]Tand
e = 2.However, the residual associated with
Then
0 01
8/2/2019 Dan Um Jay A
345/560
r = 0.01.Thus, the error is 200 times larger tharesidual.
Condition Numbe
The condition number of a matrix A is
8/2/2019 Dan Um Jay A
346/560
K(
A) =
A
A1
.
Theorem
Let A be a nonsingular matrix, X be ai l i h li
8/2/2019 Dan Um Jay A
347/560
approximate solution to the linear sys
AX = b, r = AX b and e = X X.any natural matrix norm ,
1
Ar e A1
rand
Proof
We calculate
8/2/2019 Dan Um Jay A
348/560
r = AX
b = AX
AX,
= A(X X) =Equivalently,
e = A1r.
Now let be any natural matrix no
From r = Ae, we obtain
8/2/2019 Dan Um Jay A
349/560
r
=
Ae
A
e
,
or
e
rA
.
Thus
1 1
Suppose the X = 0 and b = 0.
8/2/2019 Dan Um Jay A
350/560
Taking the norm on both sides of AX
b = AX A X.
This implies 1
X Ab .
Similarly, from X = A1b, we obtain
8/2/2019 Dan Um Jay A
351/560
1
A1 b 1
X .Then
1A1 b
1X
Ab .(2)
This provide lower and upper boundsrelative error in an approximate soluti
http://-/?-http://-/?-8/2/2019 Dan Um Jay A
352/560
relative error in an approximate soluti
AX = b,
interms of the relative residual rb and
of A and its inverse.
When K(A) is small (i.e., A A1
This shows that the relative residual pgood measure for the error in the app
8/2/2019 Dan Um Jay A
353/560
good measure for the error in the appsolution.
If K(A) is large, the relative residual cvery poor indicator of the accuracy of
approximate solution.
Example
Let 1 2
8/2/2019 Dan Um Jay A
354/560
A =
1 2
0.99 1.99
and
A1
= 199 200
99 100
.
Then1
Therefore, the relative error in an appsolution to a system with A as its coe
8/2/2019 Dan Um Jay A
355/560
solution to a system with A as its coe
matrix can be as small as
1
1197 times oas 1197 times the relative residual.
Suppose the entries in A and b are knapproximately
8/2/2019 Dan Um Jay A
356/560
approximately.
How does the error in the computed sdepend on the errors in A and b?.
Consider the original system of equatAX = b.
Let
A A + E
8/2/2019 Dan Um Jay A
357/560
A = A + E,
where E is the error matrix. ThenA + E
X = b,
where X represents true solution of t
We now estimate the error X X.
Then
8/2/2019 Dan Um Jay A
358/560
X = A1b = A1(AX) = A1(A
= X AWe have
X X = A1EX.Taking norm on both sides, we obtain
HenceX X E
8/2/2019 Dan Um Jay A
359/560
X X
X K(A)
E
A.
This shows that the relative error in thcomputed solution can be as large as
error in the coefficients of A multipliedcondition number.
Iterative Methods
It is natural to ask why we would wanneed to develop iterative techniques
8/2/2019 Dan Um Jay A
360/560
need to develop iterative techniques.
For systems of small dimension, therneed. Direct techniques will perform vefficiently.
However, linear systems arising fromapplications will frequently be quite la
For systems with large, sparse coefficmatrices direct techniques are often
8/2/2019 Dan Um Jay A
361/560
matrices, direct techniques are often efficient than iterative techniques.
Basic iterative techniques for the systsystem of equations are analogous to
fixed-point techniques.
The original linear system
Find the n-vector X such that
8/2/2019 Dan Um Jay A
362/560
AX
b = 0,
is first converted to the fixed point pro
Find the n-vector X such that
X = BX + c,
Starting from some initial approximat
solution of the fixed point problem X (
8/2/2019 Dan Um Jay A
363/560
solution of the fixed point problem X(
sequence of vectors {X(k)
} is compuaccording to the ruleX(k+1) = BX(k) + c,
where the matrix B is called the iterat
Example
Consider the system of equations
8/2/2019 Dan Um Jay A
364/560
5x1 + x2 + 2x3 = 10,
3x1 + 9x2 + 4x3 = 14x1 + 2x2 7x3 = 33
Now, we express the given system ofAX = b in the equivalent form
X = BX + c,
8/2/2019 Dan Um Jay A
365/560
as
x1
x2
x3
=
0 15 2539 0
49
17
27 0
x1
x2
x3
+
Strictly DiagonalDominant Matrice
An n n matrix A is strictly diagonall(column) dominant if for each row (co
8/2/2019 Dan Um Jay A
366/560
(column) dominant if for each row (comagnitude of the diagonal element islarger than sum of the magnitudes ofelements on the row (column), that isi (j),
|aii| >
j=1, j=i
|aij|,
Example
Let 3 1 1
8/2/2019 Dan Um Jay A
367/560
A =
3 1 1
2 6 39 7 20
.
The matrix is strictly diagonally row d
Since
|3| = 3 > 2 = | 1| + |1
Example
Let 3 2 2
8/2/2019 Dan Um Jay A
368/560
A =
3 2 2
2 6 39 7 20
.
The matrix is not strictly diagonally do
Since|3| = 3 < 4 = | 2| + |2|.
Symmetric PositiDefinite Matrices
A matrix A is Symmetric Positive Defisymmetric and
8/2/2019 Dan Um Jay A
369/560
symmetric and
XTAX > 0,
for any nonzero vector X.
Example
Let 3 1 1
8/2/2019 Dan Um Jay A
370/560
A =
3 1 1
1 4 21 2 5
.
Let X = [x1 x2 x3]T. Then
XTAX = 3x21 + 4x22 + 5x
23 + 2
8/2/2019 Dan Um Jay A
371/560
X AX 3x1 + 4x2 + 5x3 + 2
+ 4x2x3 2x1x3,= x21 + x
22 + x
23 + (x1 +
+ (x1
x3)
2 + 2(x2 +
which is clearly greater than zero for anon-zero X. Therefore, A is symmetr
Example
Consider the matrix
8/2/2019 Dan Um Jay A
372/560
A =
3
2
1
2 3 21 2 3
,
is not symmetric positive definite. Sin
XTAX = 2[(x1 x2)2 + (x1 x3)2 +
Theorem
Let AX = b can be written as X = BXsome norm of B, B < 1, then X =
8/2/2019 Dan Um Jay A
373/560
, ,a unique solution.
Further, the sequence {X(m)} genera
X(m+1) = BX(m) + c,
starting with some initial X(0) will con
We consider the homogeneous systeLet X = 0 is the solution. Then
8/2/2019 Dan Um Jay A
374/560
X = BX B XSince B < 1 implies that X < Xcontradicts that X = 0.This implies that X = 0.
This shows that the nonhomogeneouequations
8/2/2019 Dan Um Jay A
375/560
q
X = BX + c,
has an unique solution.
Let the unique solution X = . Then = B + c.
Then
X(m+1) BX(m) B B X
8/2/2019 Dan Um Jay A
376/560
X(m+1)
= BX(m)
B = B X
Now taking norm on both sides, we o
X(m+1)
=
B X(m)
B X(m) B 2 X(m1)
Taking limit on both sides as m ,( +1)
8/2/2019 Dan Um Jay A
377/560
limm
X(m+1)
= 0.
Since B < 1, this implies that
limmX(m+1)
= .
Hence, the sequence {X(m)} converg
Jacobi Method
The Jacobi method to solve AX = b iconverges, if the coefficient matrix A
8/2/2019 Dan Um Jay A
378/560
gdiagonally row dominant.
Proof. The system of equations in ma
a11 a12 a1n
a21 a22 a2n. . . . .
x1x2...
=
We express in the iteration form: fromequation, we get
8/2/2019 Dan Um Jay A
379/560
q , g
x1 = 0x1 a12a11
x2 a13a11
x3 a1na11
x
From the second equation, we get
x2 = a21a22
x1 0x2 a23a22
x3 a2na22
Similarly from the ith equation, we ob
8/2/2019 Dan Um Jay A
380/560
xi = ai1
aii x1 ai2
aii x2 ai3
aii x3 a
a
Now, we express in matrix form
8/2/2019 Dan Um Jay A
381/560
X =
0 a12a11
a1na11
a21a22
0 a2na22
......
......
...
an1ann
an2ann
0
X
The Jacobi iteration scheme is
X(m+1) BX(m) +
8/2/2019 Dan Um Jay A
382/560
X(m+1) = BX(m) + c,
convergent provided B < 1. Now
B = max1inn
j=1, j=i |
aij
||aii| | 2 | | 3 | .. | ndominant eigenvalue is not repeated
8/2/2019 Dan Um Jay A
386/560
real. Let
Y0 = 1x1 + 2x2 + . . . + nx
Generate the sequence {Yk}, whereYk = AYk1 = A
kY0.
Let xi = (xi1, xi2, . . . , xin).
L d h f
8/2/2019 Dan Um Jay A
387/560
Let Yki denote the ith component of y
Yki = 1k1x1i + 2
k2x2i + . . . +
ConsiderYk+1,i
Yk,i=
1k+11 x1i + 2
k+12 x2i + . . . +
1k1x1i + 2
k2x2i + . . . +
We can take 1 = 0 and x1i = 0.
Y 1 + d (2 )k+1 + + d
8/2/2019 Dan Um Jay A
388/560
Yk+1,i
Yk,i = 1{
1 + d2(21
)k+1 + . . . + d
1 + d2(21 )k + . . . + d
where dp =pxpi1x1i
for p = 2, 3, . . . , n .
Therefore
Yk i
Therefore we can take 1 =Yk+1,iYk,i
for s
large k
8/2/2019 Dan Um Jay A
389/560
large k.
Also
Yk =
k1{
1
x1 +
2(
2
1)
kx2 +
. . .+
This shows that for large k, Yk is prop
However, if | 1 |< 1, || Yk || 0 as k
if | | || Y || k
8/2/2019 Dan Um Jay A
390/560
if | 1 |> 1, || Yk || as k .
To avoid this in actual computation we
Let0 = 1x1 + 2x2 + . . . + n
and Z1 = A0.
In general, Zk+1 = Ak andk+1 = maxr | (Zk+1)r |. Therefore,
8/2/2019 Dan Um Jay A
391/560
Zk+1 = Ak = 1k.k1 . . . . 1
Ak
(Zk+1)i(k)i
= (Ak+10)i(Ak0)i
1k+11 x1i+2
k+12 x2i+...+n
k+1n xni
Hence (Zk+1)i(k)i 1 as k , provide
and x1i = 0.
8/2/2019 Dan Um Jay A
392/560
and x1i 0.
Also Zk+1 = Ak =1
k.k1....1Ak+10
= k+11
k.k1....1{1x1+2x2(
21
)k+1+. . .+
Hence Zk+1 is proportional to x1 for la
Example
A =
2 2 1
4 8 1
8/2/2019 Dan Um Jay A
393/560
A =
4 8 1
1 2 0
Eigen values are = 1, 3, 6.
Suppose 0 = (1, 2, 4)T.
z1 = A0 = (10, 16, 5)T.
8/2/2019 Dan Um Jay A
394/560
( )1 = 16
.1 = (.625, 1, .3125)
T
z2 = A1 = (2.9375, 5.1875, 2.625)T.
2 = 5
.1875
.
2 = (.56626, 1, .506)T
z3 = A2 = (2.62652, 5.22896, 2.56626(z3)2
Case 2
If numerically largest eigen value 1 imultiplicity r > 1.S
8/2/2019 Dan Um Jay A
395/560
Suppose1 = . . . = r, | r |> | r+1 | . . . |
Case 3
If two eigen values of maximum modubut of opposite sign.
8/2/2019 Dan Um Jay A
396/560
Example
A =
1 2 2
2 1 2
8/2/2019 Dan Um Jay A
397/560
A =
2 1 2
2 2 1
Eigen values are = 3, 3, 3.
Eigenvalues and Eigenvect
Dr. P. Dhanumjaya
8/2/2019 Dan Um Jay A
398/560
Dr. P. Dhanumjaya
Jacobi Method
Jacobi method is used to find the eigeand eigenvectors of a real and symm
8/2/2019 Dan Um Jay A
399/560
A matrix A = [aij ] is called symmetric
Similarity Matrice
If A and B are n n matrices, then BA if there exists a non-singular matrix
8/2/2019 Dan Um Jay A
400/560
B = P AP1.
We note that P AP1 is called similari
transformation.
Similarity transformation retains the seigenvalues but different eigenvectors
Note 1
A real symmetric matrix is diagonaliz
If A is any real symmetric matrix then
8/2/2019 Dan Um Jay A
401/560
If A is any real symmetric matrix, then
exists an orthogonal matrix P such th
PTAP = D (diagonal matrix
or
Ap = P D.
ElementaryOrthogonal Matri
Let the matrix Onn such that
Opp = cos = Oqq
8/2/2019 Dan Um Jay A
402/560
Opp cos Oqq,
Opq = sin ,Oqp = sin 1 p < q
The remaining entries of O are samethe unit matrix is called an elementarorthogonal matrix.
The Jacobi method consists of obtainmatrix P as a product of a sequence elementary orthogonal matrices i e
8/2/2019 Dan Um Jay A
403/560
elementary orthogonal matrices. i.e.,
P = O1 O2 On.In Jacobi method take elementary ortmatrices O1, O2, , On and generatesequence {Ar} as follows:
Take
A0 = A,
8/2/2019 Dan Um Jay A
404/560
A A,
Ai = OTi Ai1Oi, i = 1, 2, 3,
until A(r) becomes almost diagonal.
Eigenvalues of A can be approximate
diagonal entries of A(r).
Suppose at some step, we want to anpq element p < q of the current matrix
element is numerically largest out
8/2/2019 Dan Um Jay A
405/560
pq element is numerically largest out
off-diagonal entries.
The required transformation matrix O
following form:
Opp = Oqq = cos ,
Let
B = OTAO.
8/2/2019 Dan Um Jay A
406/560
Then the elements of B are as followbpk = apk cos aqk sin bqk = apk sin + aqk sin ,
bip = aip cos aiq sin ,biq = aip sin + aiq cos ,
bik = aik,2 2
8/2/2019 Dan Um Jay A
407/560
bpp = app cos
2
+ aqq sin
2
2apq sbqq = app sin
2 + aqq cos2 + 2apq c
Note that B is also symmetric. Thereorder to make bpq = 0, we must have
1(a a )sin2 + a cos2
This implies that
t 22apq
if
8/2/2019 Dan Um Jay A
408/560
tan2 =pq
(aqq app), if aqq
=
and
=
4 , if aqq = app.
We now compute the values of sin a
Now, if
8/2/2019 Dan Um Jay A
409/560
, > 0, sin =
2
> 0, < 0, sin =
2
< 0, > 0, sin = +
2
if = 0, then sin = 12
.
2
8/2/2019 Dan Um Jay A
410/560
We can compute cos = 1 sin2
.
We note that since 4 < 4 ,
12
< sin 12
.
Example 1
Use Jacobi method to find all eigenvacorresponding eigenvectors of the giv
8/2/2019 Dan Um Jay A
411/560
A =
1 2 22 3
2
2
2 1
.
Eigenvalues a
8/2/2019 Dan Um Jay A
412/560
Eigenvectors
Eigenvalues & Eigenvectors
Find a scalar and its corresponding nonzero
for a given square matrix A such that
8/2/2019 Dan Um Jay A
413/560
for a given square matrix A such that
If A is n x n matrix, then the value , for
equation (1) has non trivial solution is
eigenvalue of A and corresponding solution v
called eigenvector of A.
Ax x=
Example 1
Eigenvector of a 2
2 Matrix The vector is an eigenve
1
2
=
x
8/2/2019 Dan Um Jay A
414/560
Corresponding to the eigenvalue
since
2
3 0
8 1A
=
Examples
8/2/2019 Dan Um Jay A
415/560
Characteristic polynomial
8/2/2019 Dan Um Jay A
416/560
Example
8/2/2019 Dan Um Jay A
417/560
Characteristic polynomial
Computing eigenvalues using characteristicpolynomial is not recommended because of
8/2/2019 Dan Um Jay A
418/560
polynomial is not recommended because of
work in computing coefficients of characteristic po
sensitivity of coefficients of characteristic polynom
work in solving for roots of characteristic polynom
Characteristic polynomial is powerful theoretbut usually not useful computationally
Example
8/2/2019 Dan Um Jay A
419/560
Power Iteration
8/2/2019 Dan Um Jay A
420/560
The Power Method
Start with some random vector v, ||v
8/2/2019 Dan Um Jay A
421/560
Iterate v=(Av)/||Av|| The eigenvector with largest eigenva
tends to dominate
How fast?
Linear convergence, slowed down by ceigenvalues
Power Method
Case 1:Let
(that is dominant eigenvalue is not repeated.)
S ppose A has only eigenvalue of largest
||...|||||| 321 n >
8/2/2019 Dan Um Jay A
422/560
Suppose A has 1
only eigenvalue of largest and that 1 is real.
Consider the iteration formula:
xk+1
= Axk
where we start with some initial x0, so that:
xk
= Ak x0
Then xk
converges to the eigenve
Example of Power Method
Consider the follow matrix A
8/2/2019 Dan Um Jay A
423/560
=
100
120
014
A
Assume an arbitrary vector x0 = { 1 1 1}T
Example of Power Method
Multiply the matrix by the matrix [A] by {x}
8/2/2019 Dan Um Jay A
424/560
=
1
3
5
1
1
1
100
120
014
Normalize the result of the product
=
.0
.0
1
5
1
3
5
Example of Power Method
=
06.41
6.4
=
1
6.4
6.0
1
120
014
8/2/2019 Dan Um Jay A
425/560
=
0435.04783.0
2174.4
0435.0217.0
1
100120
014
.02.0
2.02.0100
2174.4
Example of Power Method
=
2165.0
1134.4
1134.0
1
120
014
8/2/2019 Dan Um Jay A
426/560
0103.00183.0100
=
1134.4
0103.0
2165.0
1134.4
As you continue to multiple each successiv
Power method
The special advantage of the power met
8/2/2019 Dan Um Jay A
427/560
The special advantage of the power metis that the eigenvector corresponds to thdominant eigenvalue and is generated asame time.
The disadvantage is that the method onl
Power Method:
Case 2: If numerically largest eigenvalue is multiple.
Suppose1 2 1... , | | | | ... | |r r r n += = = >
8/2/2019 Dan Um Jay A
428/560
Example:
(i) x0= [2,1]T, (ii) x0 = [1,0]
T, (iii) x0 = [0,1]T
4 1
, eigenvalues are5,5with eigenvector1 6A
=
Power Method:
80
1
4
5
5
1,5,
4
5
1
2
61
14
,61-
14
A,1
2
111
010
=
==
=
=
=
=
=
xz
Axzx
8/2/2019 Dan Um Jay A
429/560
19913.4114
,79167.4,75002.479167.4
79167.01
61-14
79167.0
1,8.4,
8.3
8.4
8.0
1
61-
14
8.0454161-
333
222
==
=
=
==
=
=
xz
xz
Power Method:
Case 3: If two eigenvalues of maximal modulus aropposite signs.
8/2/2019 Dan Um Jay A
430/560
Example :TT )3,1(x,2;)1,1(x2,,
1-3
11A ====
=
=
=
=
=
=
=
=
1-
14
4-
4
4
0
1-3
11
1
044
0
1-
1
1-3
11
1-
1
2
10
x
xx
Both the vectors are oalternatively, which ind
Power Method
=
=
=
=
1
0,
1-
1,
1
0,
1-
13210 xxxx
8/2/2019 Dan Um Jay A
431/560
=
=
=
+
=+
3-
1
1
02
1-
12
11
102
1-12
12
12
xx
xx
=
=
1
12
2
2
1
1
1-3
11
Power Method
Case 4: Dominant eigenvalue is complex conjugate of each oth
Example: =
= )1,1,1(,
5-23
411-
2-23
A 0xT
8/2/2019 Dan Um Jay A
432/560
=
=
=
==
=
=
==
=
=
,0588.3,
0 615437
1
0.36537
0588.3
1 8824
3.0588
1.1176
A
1
0.0588
1
,25.4,
1
0.0588
1
25.4
4.25
0.25
4.25
A
0
10.75
,4,
0
10.75
4
0
43
A
32
221
110
xx
xx
xx
Power Method:
Suppose numerically largest eigenvalue andcorresponding eigenvectors has been computed
power method.
8/2/2019 Dan Um J
Recommended