[noise]
[vocalized-noise] Welcome. In the [noise]
last few sessions, I try to [noise] introduce
[noise] Jacobi and Gauss-Seidel methods, which
are [noise] the basic iterative solvers [vocalized-noise].
And then, [vocalized-noise] we also discussed
about the [vocalized-noise] properties of
the matrix A [vocalized-noise] [noise] for
which this [noise] solvers should [noise]
work. And we have seen that the matrix [noise]
has to be diagonally dominant or [vocalized-noise]
[noise] irreducibly diagonally [vocalized-noise]
dominant [vocalized-noise] that means, diagonal
dominant means ah all the diagonal terms must
be greater than sum absolute value of all
the diagonal terms must be greater than [noise]
sum of absolute value of the [noise] [vocalized-noise]
off diagonals [noise] for each particular
row [vocalized-noise] [noise].
And [vocalized-noise] diagonal [vocalized-noise]
irreducibly [noise] diagonally dominance means,
the [vocalized-noise] that sum [noise] the
[noise] all ah [vocalized-noise] absolute
value of ah [vocalized-noise] diagonal term
[vocalized-noise] must be greater than [noise]
equal to sum of the [noise] absolute values
[vocalized-noise] of the off diagonal terms
for each particular row [noise]. And there
should be at least one row [noise] for which
it is not [vocalized-noise] greater than [noise]
equal to [noise] rather it is [noise] [vocalized-noise]
greater than the [noise] sum of [noise] the
[noise] off diagonal terms, [noise] that I
was absolute value [noise] of the diagonal
term [noise] is greater than the [noise] sum
of the [vocalized-noise] [noise] absolute
value of the [noise] off diagonals [noise]
[vocalized-noise].
And for this cases, [noise] Jacobi and Gauss-Seidel
works [noise]. And Jacobi and Gauss-Seidel
[noise] works for this cases, because [vocalized-noise]
we looked into the [noise] iteration matrix
G, [noise] iteration matrix [vocalized-noise]
G [vocalized-noise] will have a spectral radius
[noise] less than 1, [noise] [vocalized-noise]
which ensure we will ensure that [noise] ah
the matrix norm of [vocalized-noise] G [noise]
will be [noise] less than 1. And [vocalized-noise]
as we will [vocalized-noise] [noise] do a
number of iterations, [noise] [vocalized-noise]
the initial error [noise] which will be multiplied
[noise] by matrix norm [noise] of G to the
power k; [noise] k is a large number, which
[vocalized-noise] [noise] this will be a very
small number [noise]. Therefore, it should
converge [noise] the error should be [vocalized-noise]
very close to [noise] 0, and it should converge
[noise] to the right solution [noise].
And if it is [vocalized-noise] [noise] converges
[noise] that means, if the guess [vocalized-noise]
[noise] after now a large number of iteration,
[noise] the solution converges [vocalized-noise]
[noise] that means, the guess an [noise] updated
[noise] value has [vocalized-noise] [noise]
very small difference in between them. It
is [noise] [vocalized-noise] must [noise]
converge to the [vocalized-noise] right solution
[noise] that is also property of [noise] Jacobi
[noise] and Gauss-Seidel [noise] [vocalized-noise].
So, now what we will do, we will do some numerical
[noise] experiments using Jacobi [noise] and
Gauss-Seidel [noise]. We also introduce the
terms [vocalized-noise] convergence rate,
[noise] convergence factor, and I have seen
[noise] that [vocalized-noise] it [noise]
depends on the matrix G, [noise] it depends
how is [vocalized-noise] splitting the [noise]
main matrix say [noise] the only condition
for the [vocalized-noise] we are solving [noise]
A x is equal to b. The [noise] only condition
on the main matrix A is that, it has to be
[vocalized-noise] diagonally dominant or irreducibly
diagonally dominant [vocalized-noise].
Now, it depends how are we [noise] [vocalized-noise]
splitting the matrix A [noise] to get the
matrix G. And [vocalized-noise] based on [noise]
that, we will get [noise] [vocalized-noise]
the ah eigenvalues [vocalized-noise] of G
the [noise] largest eigenvalue is the spectral
[noise] radius, [vocalized-noise] which is
also [noise] a measure of the [vocalized-noise]
[noise] convergence factor [noise] ah [noise]
as well as convergence rate [vocalized-noise]
of the [noise] [vocalized-noise] iterative
scheme.
So, [vocalized-noise] [noise] and [vocalized-noise]
[noise] we will tell us [noise] how fast the
iterations will converge. [vocalized-noise]
[noise] So, what [noise] will do, [noise]
we do [vocalized-noise] is called ah [noise]
numerical experiments [noise] that means,
we will [vocalized-noise] take [vocalized-noise]
few matrices [noise] [vocalized-noise] and
[noise] do experiment [noise] with these matrices
[noise]. So, we will try to [vocalized-noise]
ran Jacobi [noise] and Gauss-Seidel [noise]
on this matrices [vocalized-noise].
As well as we will see, how [vocalized-noise]
are there [vocalized-noise] [noise] [vocalized-noise]
how are there [noise] [vocalized-noise] G
matrices, [noise] what are the [vocalized-noise]
[noise] spectral radius of G matrix, what
is the convergence factor [noise] or convergence
rate, [noise] [vocalized-noise] so that [noise]
we can say a priory [noise] looking into the
G matrix, we should be able to tell that [vocalized-noise]
[noise] whether they will converge fast, [noise]
whether they will at all converge [noise]
or not etcetera [vocalized-noise] [noise].
So, for [noise] [vocalized-noise] this numerical
[noise] experiments, [noise] [vocalized-noise]
we have chosen ah [noise] [vocalized-noise]
particular matrix, [noise] which is the matrix
generated [noise] as finite difference [noise]
[vocalized-noise] approximation of d 2 T d
x square is equal to [noise] 0. We have discussed
[vocalized-noise] [noise] about how to convert
[noise] a differential equation [noise] into
a [noise] matrix equation, [noise] exactly
[noise] data has been [noise] done here [vocalized-noise].
And the boundary condition T [vocalized-noise]
x is 0 to 1, [noise] boundary condition is
T is equal to 0, it is 0, T is [noise] equal
to 1 that [vocalized-noise] [noise] value
is 1 [vocalized-noise]. And [noise] we now,
[vocalized-noise] so this will [noise] give
essentially a tridiagonal system, [noise]
we varied the number of [vocalized-noise]
intervals and got large number of points,
and got [vocalized-noise] larger matrix [noise]
systems [noise] [vocalized-noise].
So, a [noise] tridiagonal matrix [noise] is
obtained [vocalized-noise] [noise]. Matrix
size is varied [noise] by num[ber] [vocalized-noise]
varying the number of grid points or number
of intervals, [vocalized-noise] [noise] and
we got different [noise] matrices ah [vocalized-noise]
[noise]. But, the matrices [noise] will essentially
look like [noise] [vocalized-noise] a tridiagonal
matrix [noise] for this is an example of a
[vocalized-noise] [noise] the 8 by 8 [noise]
matrix [noise] [vocalized-noise]. However,
here we consider [noise] 10 by 10, [noise]
20 by 20, [noise] 40 by 40 [noise] matrices
[noise] [vocalized-noise].
And [noise] then with [noise] this matrix,
[noise] we [vocalized-noise] ran [vocalized-noise]
ran a Gauss-Seidel [noise] or Jacobi [vocalized-noise]
program. ah Later, we will also discuss how
to [vocalized-noise] [noise] write this program
[noise]. So, we have to [noise] implement
a [noise] Gauss-Seidel or Jacobi into a [vocalized-noise]
computer program, we [noise] ran a program
[vocalized-noise].
And noted, [noise] what are the number of
iterations [noise] to solve it [noise] [vocalized-noise].
Also noted into the [noise] [vocalized-noise]
G matrix, [noise] which comes out of [vocalized-noise]
so, this [vocalized-noise] this is basically
[noise] [vocalized-noise] how the [noise]
A matrix [noise] looks like here [noise] [vocalized-noise].
But, this is 8 by 8, [noise] we have taken
several [noise] matrices [vocalized-noise]
[noise]. So, seen [noise] how is the G matrix,
and then [vocalized-noise] try to find out
the eigenvalues [noise] of the G matrix from
which we can get the spectral radius, and
[noise] [vocalized-noise] see [noise] how
is the convergence property [noise] [vocalized-noise].
So, matrices are solved [noise] using [noise]
Jacobi and Gauss-Seidel method, [vocalized-noise]
because this is [noise] doing experiment [noise]
with the matrices and the method [noise] [vocalized-noise]
And but, this is not experiment in hardware,
we call it [noise] numerical experiment [noise]
[vocalized-noise]. The spectral radius [noise]
of matrix G, convergence rate [noise] [vocalized-noise]
convergence rate is [noise] [vocalized-noise]
ah log minus [noise] logarithm of ah [noise]
of the [noise] convergence factor [noise]
of the spectral radius, [noise] and the number
of iteration, [vocalized-noise] and [noise]
noted for different sizes of the matrices
[noise] [vocalized-noise].
So, for [noise] 10 by 10 matrix [noise] [vocalized-noise]
for [noise] Jacobi iteration, the spectral
radius is 0.94632; [vocalized-noise] [noise]
for Gauss-Seidel, [noise] the spectral radius
[noise] is 0.89533 [noise]. So, [noise] for
Gauss-Seidel, [noise] the spectral radius
is [noise] smaller [vocalized-noise] [noise].
The convergence [noise] rate is 0.05517 for
[noise] Jacobi, and for Gauss-Seidel [noise]
is 0.11034 [vocalized-noise]. So, convergence
rate is [vocalized-noise] almost [noise] double
[noise] in case of [noise] Gauss-Seidel than
Jacobi. And this is due to the fact, that
the [noise] spectral radius is smaller [vocalized-noise].
So, if the convergence is [vocalized-noise]
double, [noise] the rate will be double, [noise]
speed will be [noise] faster [noise]. So,
the number of [vocalized-noise] iterations
will be small, [noise] faster it will converge
[noise]. And we can see that [noise] [vocalized-noise]
this is [noise] almost have, [noise] if we
need [noise] 200 iterations [noise] for Jacobi
[noise] for Gauss-Seidel, [noise] we need
106 iterations [noise]. So, there is a relation
between convergence rate and number of [vocalized-noise]
iterations [vocalized-noise]. Roughly if we
multiply convergence [noise] rate with [noise]
number of iterations, [noise] there [noise]
[vocalized-noise] there [noise] that is roughly
[noise] constant, [noise] at least in this
case [noise]. In other cases, [noise] we will
also see [noise] [vocalized-noise] at least
for same size of matrix [vocalized-noise]
[noise].
So, [vocalized-noise] once we have a method
[vocalized-noise] [noise] in a Gauss-Seidel,
[noise] we have a laser spectral [noise] radius,
therefore the convergence will be faster [noise]
the number of iterations [noise] [vocalized-noise]
so, rho G is [noise] [vocalized-noise] also
convergence [noise] factor right [noise] rho
G is also [noise] convergence factor [noise]
factor [noise] general [noise] general convergence
factor [noise] [vocalized-noise]. So, if [noise]
this is smaller, [noise] the iterations [noise]
will [vocalized-noise] we need [noise] less
number of iterations, the [noise] scheme will
be faster [noise] [vocalized-noise].
Now, we go to 20 into 20, [noise] we can see
similar thing [noise] [vocalized-noise]. The
spectral radius [noise] [vocalized-noise]
of Jacobi is 0.98645, [noise] cell in not
lot more iterations [noise] spectral radius
is very close to 1 [vocalized-noise]. The
convergence [vocalized-noise] rate is 0.101364
[noise]. The Gauss-Seidel [noise] spectral
[noise] radius is 0.973, [vocalized-noise]
convergence [noise] rate is 0.027 [vocalized-noise].
And the number of iterations are [noise] in
almost close to [noise] half of each other
709 [noise] in Jacobi, [noise] and [noise]
374 in [noise] [vocalized-noise] ah Gauss-Seidel
[noise] [vocalized-noise].
And now, [vocalized-noise] [noise] we [noise]
can also see [noise] the [vocalized-noise]
also [noise] you can see that [noise] [vocalized-noise]
for convergence rate [noise] in 10 into 10
Jacobi [noise] is 0.055, [vocalized-noise]
it added in [vocalized-noise] ah Gauss-Seidel
[noise] is 0.013 [noise] [vocalized-noise]
nearly four times, [noise] convergence rate
is [vocalized-noise] reduced by nearly four
times. And [noise] iteration number is [noise]
increased [noise] by 23.5 times [vocalized-noise].
So, there is a [vocalized-noise] relation
between convergence [noise] rate and iteration
number for similar type [noise] of matrices
[vocalized-noise].
If we increase 40 into 40, [noise] the Jacobi
[noise] converse spectral radius is very high
[noise] 0.99661, [vocalized-noise] convergence
[noise] radius [noise] is 0.0034 [noise].
Gauss-Seidel 0.99324, [vocalized-noise] [noise]
convergence [noise] radius is 0.0068, [noise]
[vocalized-noise] because the iteration [noise]
number is very high [vocalized-noise]. I think
this is [noise] [vocalized-noise] swiped,
[noise] this is swiped [noise] [vocalized-noise]
this will be [noise] well ah [vocalized-noise]
change [noise] that in the [noise] later slide,
[noise] [vocalized-noise] this is 2435 [noise]
this is 1291 [noise] this is [noise] (Refer
Time: 09:02) while [noise] making the [noise]
table [noise] [vocalized-noise].
So, [noise] because the iteration [vocalized-noise]
spectral radius is very high very close to
1, [vocalized-noise] we need a [vocalized-noise]
[noise] very large number of iterations [noise]
[vocalized-noise] to 2435 [vocalized-noise]
[noise] iterations here [vocalized-noise].
And in Gauss-Seidel, [noise] we need [noise]
[vocalized-noise] it is also high, [noise]
we need 1200 [noise] iterations [vocalized-noise].
So, [noise] as the spectral radius [noise]
is increasing, we are also needing [noise]
large number of iterations [vocalized-noise].
However, Jacobi consistently give [noise]
[vocalized-noise] larger spectral radius than
Gauss-Seidel, [noise] and it consistently
[noise] takes more number of iterations than
[noise] Gauss-Seidel. And the convergence
rate is reversed [noise] in both the cases
[noise] [vocalized-noise].
So, the [noise] observations are, [noise]
[vocalized-noise] largest [noise] eigenvalues
[noise] or the spectral radius increases [noise]
with matrix size [noise] [vocalized-noise].
Matrix 20 into 20, [vocalized-noise] eigenvalue
[noise] largest [noise] eigenvalue is more
[noise] [vocalized-noise]. And [noise] spectral
radius is and [vocalized-noise] because [noise]
largest eigenvalue [noise] increases, [noise]
spectral radius [noise] ah also increases,
and number of iterations [noise] are more
[vocalized-noise]. Spectral radius [noise]
is higher for Jacobi, [noise] then [vocalized-noise]
than [noise] Gauss-Seidel [noise] [vocalized-noise].
And we [vocalized-noise] [noise] when [noise]
proposing [noise] Gauss-Seidel, almost intuitively
[noise] you propose that [noise] if we can
use some [noise] of the iterative values some
of the updated [noise] values, the iteration
should [noise] converge faster [vocalized-noise]
[noise]. So, this is an almost [vocalized-noise]
[noise] intuitive [vocalized-noise] all [vocalized-noise]
actually [noise] intuitive [vocalized-noise]
proposition, [noise] however we can see that
[noise] it is actually giving faster [noise]
convergence [vocalized-noise]. And why is
it giving [noise] faster convergence, [noise]
because in the [noise] [vocalized-noise] iteration
matrix level, [noise] the spectral radius
has been reduced [noise] well [vocalized-noise]
[noise].
The ah [noise] 3rd observation is [noise]
that number of [noise] iterations [noise]
required for [noise] convergence [noise] is
strongly related with [noise] convergence
rate [noise] [vocalized-noise]. Almost if
we have [vocalized-noise] ah [vocalized-noise]
[noise] like here [noise] convergence rate
is almost twice [noise] for Jacobi and Gauss-Seidel,
[vocalized-noise] and the number of iterations
is again [noise] almost [vocalized-noise]
[noise] half [noise] in Gauss-Seidel [vocalized-noise]
ah compared to Jacobi [noise].
So, larger [vocalized-noise] its [noise] convergence
rate and [noise] number of iterations are
related [noise] with each other [noise]. For
a same matrix, [noise] Jacobi gives higher
spectral radius, [noise] therefore [vocalized-noise]
slower convergence rate [noise] compared to
Gauss-Seidel [vocalized-noise]. And as we
will increase matrix size, [noise] the [vocalized-noise]
largest eigenvalue increases, [noise] convergence
rate decreases [vocalized-noise].
This is actually very important [noise] the
first observation, because [vocalized-noise]
ah we will come [noise] [vocalized-noise]
ah maybe at the [noise] very end of this [noise]
course, when [noise] we will discuss about
[noise] [vocalized-noise] what is called [noise]
multi grid, [noise] we will look into [vocalized-noise]
ah large matrix size [noise] [vocalized-noise]
ah issues [noise] with [noise] larger [noise]
matrix size [noise] [vocalized-noise].
But, if we are [noise] increasing [noise]
larger matrix size, [noise] we are taking
[vocalized-noise] lot of [vocalized-noise]
grid points. You were taking lot of [noise]
grid points, you are doing [vocalized-noise]
less [vocalized-noise] numerical errors [noise]
in terms of the [noise] deceleration error
[vocalized-noise] [noise]. However, [noise]
the number of iterations increase, and the
computational cost increase [noise]. Why the
number of iterations increase, [vocalized-noise]
because the eigenvalue of the [vocalized-noise]
largest eigenvalue of the [noise] G matrix
[noise] or spectral radius [noise] of the
G matrix [noise] increases, [noise] if we
take large matrix [noise]. So, [vocalized-noise]
large matrix [noise] many number of grid points
[noise] will give us [noise] accurate result,
[vocalized-noise] but it will be a slow process,
[vocalized-noise] and the computational cost
will be high [noise] [vocalized-noise].
So, now [noise] let us [noise] look into the
[noise] convergence in detail [noise]. Let
us [noise] run a Gauss-Seidel [noise] code
[noise] for a diagonally [noise] dominant
[noise] a small matrix [noise] ah 3 by 3 [noise]
diagonally dominant matrix [noise] [vocalized-noise].
We look [noise] into the convergence in terms
of [noise] L infinity [noise] norm of the
residual b minus A x, [vocalized-noise] and
we will also [noise] look into the [noise]
difference [noise] of x k plus 1 minus x k
[noise] [vocalized-noise].
The residual initially [noise] start with
[noise] some guess x, [vocalized-noise] x
is equal [noise] to 0 [noise]. First iteration,
residual [noise] is 1.24, [noise] and [noise]
x k plus 1 [noise] minus x k [noise] is 0.45
[vocalized-noise]. The next case, [noise]
we use the updated x, [noise] [vocalized-noise]
the residual [noise] becomes [noise] 0.42,
[noise] and the difference between x [noise]
[vocalized-noise] x k plus [vocalized-noise]
ah 1 minus x k [vocalized-noise] becomes [noise]
0.248.
So, as [noise] this difference [noise] as
this difference decreases [noise] we can see
[noise] as this difference decreases, [noise]
the residual [noise] also decreases [noise]
[vocalized-noise]. And this difference decreases
[noise] asymptotically [noise] at the end,
[noise] x k plus 1 [noise] minus x k [noise]
is 10 to the power minus 9 [noise] and the
residual is [vocalized-noise] [noise] ah 10
to the power minus 8 [noise].
We will say that [noise] [vocalized-noise]
the difference [noise] this uh [noise] this
has converged [vocalized-noise] [noise]. So,
A x minus [noise] ah the difference between
[vocalized-noise] b minus A x [noise] which
should be actually [noise] 0 is also [noise]
practically 0, [noise] because we are solving
[noise] [vocalized-noise] A x [noise] is equal
to b [noise] [vocalized-noise]. So, this should
also [noise] come to 0, [noise] [vocalized-noise]
as well as this also [noise] [vocalized-noise]
should go to 0 [noise] for convergence and
accuracy [noise] [vocalized-noise]. So, as
a [noise] solution becomes [noise] accurate,
the ah [vocalized-noise] difference goes to
0 [noise] [vocalized-noise].
Now, [noise] then make [noise] few [noise]
[vocalized-noise] observations from [noise]
that [vocalized-noise]. That both the parameters
[noise] x k plus 1 [noise] minus x k [noise]
x k plus 1 minus x k [vocalized-noise] and
residual, [noise] [noise] both these parameters
[noise] monotonically reduced [noise] to 0
[noise]. These are the L infinity norm [noise]
that is the [vocalized-noise] largest value
of the [vocalized-noise] both [noise] both
these are vectors, [noise] so largest [noise]
component of the [noise] vector [vocalized-noise]
[noise]. Both these [noise] parameters [noise]
monotonically [vocalized-noise] reduce to
0, [noise] there is no oscillation [noise]
[vocalized-noise]. It is [noise] it is never
increasing, [noise] it is continuously [noise]
reducing to 0.
So, [vocalized-noise] for Gauss-Seidel or
Jacobi [noise] or all [vocalized-noise] something
like these [noise] direct basic [noise] iterative
methods, [noise] [vocalized-noise] the errors
[noise] or the [vocalized-noise] residuals
must not [vocalized-noise] oscillate, they
must [noise] monotonically go to 0 [noise].
If uh [noise] if there is something different,
if it is [noise] oscillating, [noise] then
there is [vocalized-noise] some problem in
the implementation [noise]. Why, because every
time this residual [noise] is being multiplied
by G to [noise] the power k, [noise] which
is [vocalized-noise] reducing the [noise]
[vocalized-noise] norm of the vector [vocalized-noise]
[noise].
The difference x k plus 1 minus x k is [noise]
correlated with r is equal to b minus [noise]
ah [vocalized-noise] A x k [noise] b minus
A x k [noise]. As [noise] the residual and
[noise] the difference is correlated in a
sense, if [vocalized-noise] this is reducing,
[vocalized-noise] [noise] this is also reducing
[vocalized-noise]. Like as the smaller [noise]
the value of x k plus 1 minus x k is [noise]
being like, this [noise] is smaller than [noise]
this value [vocalized-noise]. Therefore, the
residual is also smaller [noise] than this
value, so they are [noise] related [noise]
[vocalized-noise].
Interestingly, [noise] the residual as well
as this difference, the [vocalized-noise]
rate of slowing down [noise] is not [vocalized-noise]
reduces with time [noise] [vocalized-noise].
For example, here [vocalized-noise] this difference
[vocalized-noise] reduces [noise] by 0.21,
[noise] in the next step this [vocalized-noise]
this difference reduces [noise] by [noise]
0.17 something like that [vocalized-noise].
Here this difference reduces by zero [noise]
point ah [noise] 3 into 10 to the power [noise]
minus 7 [vocalized-noise] [noise]. Here this
difference reduces [noise] by something order
of 10 to the power minus 9 [noise] [vocalized-noise].
So, [noise] [vocalized-noise] the rate at
which the differences are reducing is also
slowing down with number [noise] of iterations
[vocalized-noise].
So, if I [noise] say that the [noise] that
this is the x star [noise] and I start up
[noise] started with x 0, [noise] and this
[noise] is the [vocalized-noise] number of
[noise] iterations, [noise] this is x star
this this is x 0 [vocalized-noise]. This [noise]
kind of [noise] asymptotically approaches
[noise] it, the rate [noise] slows down [noise]
with time [noise]. It is an important [noise]
[vocalized-noise] observation here [vocalized-noise].
However, so [noise] this is [noise] what were
we gating using [noise] Gauss-Seidel [noise]
or this uh [noise] this this particular case
is [noise] [vocalized-noise] I think for ah
[noise] [vocalized-noise] [noise] (Refer Time:
16:28) this particular [noise] case is for
a [noise] Gauss-Seidel solver, [noise] these
for a [noise] Gauss-Seidel solver. [noise]
For Jacobi [noise] it will give [noise] [vocalized-noise]
similar result, [noise] but we will need [vocalized-noise]
a more number [noise] of iterations [vocalized-noise]
the [vocalized-noise] the rate at which [noise]
the residual falls down is [vocalized-noise]
smaller [noise] because G has a [noise] larger
[vocalized-noise] at the [noise] spectral
radius [vocalized-noise]. However, [noise]
so [noise] what do you see that [noise] they
are [noise] following certain pattern [noise]
[vocalized-noise] and this [vocalized-noise]
this pattern is determined [noise] by how
is G [noise] what are the eigenvalues of [noise]
G [vocalized-noise] ah [vocalized-noise] being
a [vocalized-noise] largest eigenvalue of
G [noise] in [noise] Gauss-Seidel [noise]
or Jacobi [noise] method [noise] [vocalized-noise].
Now, our question [noise] becomes [noise]
that [noise] can we [noise] improve it, [noise]
[vocalized-noise] can we [noise] [vocalized-noise]
increase the spectral [noise] radius [noise]
through some [noise] splitting will not be
[vocalized-noise] [noise] intuitively [noise]
we found out 1 splitting which is [noise]
Gauss-Seidel [vocalized-noise] can we have
[noise] some other splitting through which
[noise] we can improve the [vocalized-noise]
[noise] convergence rate [noise] [vocalized-noise].
How to improve [noise] convergence? Jacobi
or Gauss-Seidel [noise] gives [noise] x k
[noise] plus 1 [noise] is equal to [noise]
G x k plus f [vocalized-noise] [noise] we
will think of some [noise] and this G comes
out [noise] splitting of a is equal to a [noise]
is equal to [noise] a minus n which is a minus
d [noise] minus [vocalized-noise] d minus
[noise] n [noise] minus f [vocalized-noise]
or d minus n minus f something like that and
then [vocalized-noise] ah [noise] forming
ah [vocalized-noise] m inverse [noise] n as
G [noise] [vocalized-noise].
Now, we will [vocalized-noise] try to [vocalized-noise]
A is equal to m minus n [noise] [vocalized-noise]
m inverse [noise] n is equal to G [noise]
[vocalized-noise] that is the idea [noise]
we followed here [vocalized-noise]. Now, if
we try to improve [noise] convergence [noise]
that means, if we try to [vocalized-noise]
take m and n [vocalized-noise] in some different
way, so that the [vocalized-noise] G has a
[vocalized-noise] smaller spectral [noise]
radius [noise] can we do that [vocalized-noise].
So, again [noise] we will start [noise] with
[noise] some intuitive idea here [noise] [vocalized-noise]
and we will see [noise] the effect [noise]
in the [noise] [vocalized-noise] spectral
radius [noise] later [noise] [vocalized-noise].
So, [noise] [vocalized-noise] what we do at
every [noise] iteration step is [noise] that
x k plus 1 [noise] minus x k [noise] is equal
to [vocalized-noise] G x k [noise] minus x
k [noise] plus f [vocalized-noise].
So, [vocalized-noise] [noise] after every
iteration [noise] step [noise] [vocalized-noise]
x k plus 1 [noise] this is the iteration [noise]
step is modified [noise] this is the error
[vocalized-noise] be difference between [vocalized-noise]
consecutive [noise] ah guesses [vocalized-noise]
updates [noise]. The x k [noise] plus 1 is
[noise] modified by this term [noise] [vocalized-noise]
G minus I [noise] x k plus f [noise] [vocalized-noise].
Now, if we go to my [noise] [vocalized-noise]
previous [noise] slide [noise] [vocalized-noise]
what was looking here [vocalized-noise] is
that [noise] that this is [noise] [vocalized-noise]
this is G [noise] [vocalized-noise] minus
I [noise] x k [noise] plus f. [noise] [vocalized-noise]
Why is it taking so much [noise] time so [vocalized-noise]
convergence because [vocalized-noise] this
value is [noise] reducing [noise] with [noise]
each iteration [noise] this value is not a
large [noise] value this is a very [noise]
small value [vocalized-noise].
Had this been of [noise] [vocalized-noise]
larger value [noise] we could have [noise]
[vocalized-noise] achieve the solution [noise]
we could have conversed [noise] to the right
solution [noise] in a fast manner? [noise]
[vocalized-noise] Or [noise] ah [vocalized-noise]
I showed [noise] that [noise] this is [noise]
approaching [noise] this ah the actual solution
[noise] x star [noise] [vocalized-noise] x
0 in [noise] something like this [vocalized-noise].
Instead if we can use some [noise] [vocalized-noise]
faster [noise] solution or some faster way
to [noise] converge it [vocalized-noise] [noise].
So, what we will [noise] do here [noise] [vocalized-noise]
we will look into [noise] this difference
[noise] [vocalized-noise] and if you can [noise]
[vocalized-noise] increase this difference
because the [noise] it every time [noise]
we are [vocalized-noise] trying to approach
the [noise] final solution by [vocalized-noise]
[noise] adding some value [noise] with x k
[noise] and getting x k [noise] plus 1 [vocalized-noise].
So, if [vocalized-noise] what are we adding
here [noise] is G minus I [noise] x k plus
f [noise] in Gauss-Seidel [noise] and Jacobi.
[noise] [vocalized-noise] If we can [noise]
increase this [noise] or if we can scale this,
[noise] [vocalized-noise] so [noise] this
[vocalized-noise] this is the difference that
x [noise] covers in a particular iteration
[noise] till it [noise] convergence [noise]
to the exact solution. [noise] [vocalized-noise]
If x k plus 1 [noise] minus x k [noise] can
be increased at each iteration steps [noise]
number of iterations will be reduced [noise]
[vocalized-noise]. This goes through a very
[vocalized-noise] [noise] very simple [noise]
linear logic [noise] [vocalized-noise]. If
we can [vocalized-noise] increase this value
[noise] then x k plus 1 [noise] will be further
updated [noise]. So, it converges [noise]
to the right result [noise] there can be [vocalized-noise]
[noise] some [noise] ah [vocalized-noise]
cons also it [noise] might not converge it
may diverge. It is [vocalized-noise] increased
so much, there it is [noise] never coming
back [vocalized-noise].
Let us see [noise] what happens here. [noise]
[vocalized-noise] And this method [noise]
is called a successive over [noise] relaxation
[noise] that you over relaxed the solution
[vocalized-noise] you are updating the solution
[noise] using [vocalized-noise] some ah [vocalized-noise]
some [noise] relation [vocalized-noise] [noise]
from the [vocalized-noise] [noise] using gauge
value [noise] you are updating the [vocalized-noise]
[noise] solution [vocalized-noise]. The increase
the [vocalized-noise] [noise] amount of update
[noise] [vocalized-noise] update it more [noise]
ah [noise] over relaxation updating [noise]
[vocalized-noise] update [noise] [vocalized-noise].
So, Jacobi [noise] Gauss-Seidel step we will
[noise] look like this and over relaxation
step [vocalized-noise] will be like this value
is multiplied into [vocalized-noise] omega
G I [noise] this is this is under a bracket
[noise] [vocalized-noise] omega into [noise]
G minus [noise] [vocalized-noise] i x k [noise]
[vocalized-noise] and omega greater than 1
[noise] [vocalized-noise] [noise]. So, at
each iteration [noise] step update x as [noise]
x k plus 1 [noise] is equal to x k [noise]
[vocalized-noise] omega G minus i x k [noise]
sorry [noise] plus f [vocalized-noise] [noise]
this will be plus f [noise] [noise] (Refer
Time: 21:59) [vocalized-noise].
The SOR [vocalized-noise] [noise] if we [noise]
[vocalized-noise] write it down [noise] [vocalized-noise]
by D [vocalized-noise] D and [vocalized-noise]
I [vocalized-noise] SOR will be d minus omega
e [noise] x k plus 1 is [noise] equal to omega
[vocalized-noise] f plus 1 [noise] minus w
d [noise] [vocalized-noise] x k plus [noise]
omega b [vocalized-noise] [noise]. And the
iteration matrix [noise] [vocalized-noise]
or x k plus 1 [noise] is equal to D minus
omega inverse [noise] w F plus 1 minus or
1 [noise] D [noise] x k plus D minus omega
inverse [noise] omega b [vocalized-noise].
So, it is [noise] nothing but [noise] [vocalized-noise]
they [noise] follow this step [noise] x k
[noise] plus [noise] 1 is equal to x k [noise]
G [vocalized-noise] G minus I [noise] x k
plus x [vocalized-noise].
So, [vocalized-noise] find out the updated
value [noise] what is the difference between
the updated value and gauges [noise] [vocalized-noise]
multiplied [noise] by some factor omega which
is greater than 1 [vocalized-noise]. And then
add this with the gauge value this is your
[noise] new update value. [vocalized-noise]
[noise] And this is equivalent to [noise]
writing the matrix [noise] like this [vocalized-noise].
So, the G [noise] matrix iteration matrix
is D [noise] minus omega E inverse [noise]
[vocalized-noise] omega F plus [vocalized-noise]
I minus w omega D [vocalized-noise].
In case [noise] omega is equal to [noise]
1, [noise] this is [noise] not over relaxation
this is becomes [vocalized-noise] to the [noise]
Jacobi or Gauss-Seidel [noise] step [noise].
We can check that [noise] this D minus [noise]
E inverse [vocalized-noise] f plus D [noise]
this becomes the [noise] Jacobi [noise] step
in case [noise] omega is equal to 1. [vocalized-noise]
[noise] This is [noise] this is [vocalized-noise]
sorry Gauss-Seidel step [vocalized-noise]
over relaxation [noise] is applied over [noise]
Gauss-Seidel [vocalized-noise].
In case omega is equal to 1, [noise] [vocalized-noise]
this [noise] goes to [noise] Gauss-Seidel.
[noise] [vocalized-noise] [noise] What happens
[noise] in case omega [noise] is equal to
0 [noise] and check that [vocalized-noise]
[noise]. So, this is [noise] D inverse [noise]
f [vocalized-noise] [noise] plus D [noise]
right and [noise] f plus [noise] D is again
[noise] [vocalized-noise] ah [noise] [vocalized-noise]
D is equal to [noise] d minus c minus f is
equal to [noise] [vocalized-noise] a [noise].
So, [vocalized-noise] [noise] iteration [vocalized-noise]
does not progress [noise] it stops there [noise]
there is no [noise] update [noise] [vocalized-noise].
So, what should be [vocalized-noise] at omega
is equal to [vocalized-noise] 0 [vocalized-noise]
1, it is Gauss-Seidel; [noise] omega less
than 1 it is under relax [noise] the solution
[noise] is [noise] [vocalized-noise] become
slower [noise]. So, omega must be greater
than 1 [noise] what should be the [noise]
[vocalized-noise] right value [noise] of omega
for [noise] which we will. [noise] [vocalized-noise]
[noise].
For convergence [noise] we need [noise] [vocalized-noise]
[noise] the spectral radius [noise] of G less
than 1. [noise] [vocalized-noise] The theorem
if A is symmetric [noise] with positive diagonal
[noise] element an omega [noise] is in between
[noise] 0 to 2 [noise] successive [vocalized-noise]
over relaxation [noise] will converge for
[noise] any x 0, [noise] if A is positive
definite, [noise] [vocalized-noise] this is
a [noise] theorem like that [vocalized-noise]
[noise]. So, if omega [noise] greater than
2 [noise] successive [vocalized-noise] efficient
[vocalized-noise] will [noise] will not converge
[noise] the theorem [noise] is not satisfied
[noise]. So, it will [noise] diverge. [noise]
[vocalized-noise]
If omega is equal to 1 [noise] successive
[noise] [vocalized-noise] it is [noise] same
as basic Gauss-Seidel [noise] or Jacobi method
[noise]. If it is applied over Gauss-Seidel
[vocalized-noise] it is [vocalized-noise]
usually it is applied [noise] over Gauss-Seidel
[noise] because [noise] Gauss-Seidel is [noise]
better method than [noise] Jacobi [vocalized-noise].
And we want to further improve convergence
[vocalized-noise]. So, we imply [noise] successive
[vocalized-noise] discussion over [vocalized-noise]
the [noise] Gauss-Seidel [noise] [vocalized-noise].
If omega is [noise] less than 2, [noise] [vocalized-noise]
successive [vocalized-noise] omega is less
than 1 [noise] I am sorry, [noise] [vocalized-noise]
omega is less than 1 [noise] [vocalized-noise]
successive relaxation [vocalized-noise] is
actually [noise] under relaxing the [noise]
iterations are increasing [noise] the number
of iterations [noise] [vocalized-noise]. So,
the question is that [vocalized-noise] [noise]
omega [noise] should be between 1 to 2 [noise]
[vocalized-noise] for [noise] better convergence
[noise] [vocalized-noise] omega [noise] for
[noise] better convergence [noise] [vocalized-noise].
So, the question [noise] is [vocalized-noise]
what should be the value of [noise] the [noise]
omega [noise]. What is the optimum value [noise]
of omega? [noise] And [noise] in that [vocalized-noise]
in [noise] which case [noise] [vocalized-noise]
and what is the optimum value of omega [noise]
[vocalized-noise] when rho G is least [noise]
will have best [noise] convergence rate [noise]
[vocalized-noise]. So, to either look into
the G matrix [noise] and c [noise] with omega
c when it is having the [vocalized-noise]
[noise] least [vocalized-noise] ah [vocalized-noise]
over relaxation factor [vocalized-noise] least
[vocalized-noise] spectral radius.
So, the [vocalized-noise] so we said that
[noise] [vocalized-noise] that is the best
omega [noise] [vocalized-noise] omega is in
between 1 to 2 [noise] in 2 omega will [noise]
start diverging in [vocalized-noise] 1 it
is same as [noise] Gauss-Seidel [noise]. So,
[vocalized-noise] best omega [noise] should
be in between 1 to 2 [noise] [vocalized-noise].
However, [noise] there [vocalized-noise] is
[noise] 1 more theorem [noise] [vocalized-noise]
which [noise] helps us to determine [noise]
what should be [noise] the value of [noise]
optimum [noise] omega [noise] [vocalized-noise].
If A [noise] be a consistently [noise] ordered
matrix [noise] such that [vocalized-noise]
[noise] the diagonal element is non-zero [noise]
[vocalized-noise] ah there is a definition
[noise] of consistently [noise] ordered matrix
will [noise] not come here, but all these
matrices [noise] we deal with [noise] is usually
consistently ordered [noise] matrix [vocalized-noise]
[noise]. And omega [noise] is non-zero [vocalized-noise]
and [vocalized-noise] if [noise] [vocalized-noise]
lambda is a non-zero eigen [noise] value of
the SOR [noise] iteration matrix c, [vocalized-noise]
[noise] there is a scalar mu [noise] which
satisfies [noise] lambda plus [vocalized-noise]
omega minus 1 [noise] whole square is equal
to [noise] lambda square omega square mu square
[vocalized-noise] then mu is the [noise] eigenvalue
[noise] of the iteration matrix [noise] ah
Jacobi iteration matrix [vocalized-noise]
B [noise] [vocalized-noise].
So, Jacobi iteration matrix eigenvalue [vocalized-noise]
[noise] and the [noise] S 1 [noise] matrix
eigenvalue [vocalized-noise] are related by
the [noise] term [noise] omega [noise] [vocalized-noise].
So, [noise] and using this [noise] we can
find out for which omega [vocalized-noise]
[noise] if we know [noise] already the iteration
matrix [noise] we can [vocalized-noise] Jacobi
iteration [noise] [vocalized-noise] actual
matrix [noise] we can find out the [noise]
Jacobi iteration [noise] matrix eigenvalue
[noise] [vocalized-noise]. And find out [vocalized-noise]
for which omega [noise] we should have [noise]
[vocalized-noise] the [noise] least [noise]
value of lambda [vocalized-noise] or [vocalized-noise]
that will be the [noise] [vocalized-noise]
or that is [vocalized-noise] what eigenvalue
has least [noise] smallest spectral radius
[vocalized-noise]. Then we will say that [noise]
this is the [vocalized-noise] ah optimum omega
[noise] [vocalized-noise].
The converse is also [noise] true [noise]
if mu is an [noise] eigenvalue [noise] of
Jacobi iteration [noise] matrix Ps [noise]
[vocalized-noise] and a scalar [vocalized-noise]
lambda is there which satisfies [noise] lambda
plus 1 [noise] the omega minus 1 [noise] whole
square [noise] is equal to lambda [noise]
[vocalized-noise] omega square [noise] mu
square [noise]. Then lambda is an eigenvalue
[noise] of SOR [noise] iteration [noise] matrix
[noise] is the proofs [noise] are complex
[noise] we are not been [noise] [vocalized-noise]
discussing [noise] here [vocalized-noise].
Using the [noise] above relation [noise] we
can find out the [noise] optimum successful
[noise] relaxation factor [noise] [vocalized-noise]
omega of 2 [vocalized-noise] which is [noise]
obtained as [noise] [vocalized-noise] omega
of this 2 by [noise] root over 1 minus [vocalized-noise]
ah 1 plus root over [noise] 1 minus rho b
[noise] whole square [vocalized-noise] which
is the [vocalized-noise] this is the [noise]
[vocalized-noise] spectral radius [vocalized-noise]
of the Jacobi [noise] iteration matrix [noise]
B [noise] from the same problem [noise] [vocalized-noise].
Now, [noise] what we will do [noise] we will
[noise] do [noise] continue the numerical
experiment [noise] with [noise] what we did
earlier [noise] with different [noise] value
of omega [noise] for Gauss-Seidel [noise]
[vocalized-noise]. And what we can see is
that the case for [noise] which [vocalized-noise]
[noise] Jacobi will [noise] 20 into 20 matrix,
[noise] Jacobi needed [noise] ah 709 steps,
[noise] Gauss-Seidel needed [noise] 374 step
[noise] [vocalized-noise]. For varying omega,
[noise] omega 1.2, [noise] it is 258; [noise]
omega 1.8 it is [vocalized-noise] 133; [vocalized-noise]
omega 1.5 [noise] [vocalized-noise] on [vocalized-noise]
it is [noise] 83 steps [noise] [vocalized-noise].
From 374, [noise] it came down [noise] to
83 steps [noise] [vocalized-noise] again omega
started [noise] increasing here [vocalized-noise].
So, [noise] omega optimum [vocalized-noise]
[noise] should be in between 1.5 to 1.9 here.
[noise] [vocalized-noise]
And what is happening [noise] in the spectral
radius [noise] [vocalized-noise] at 1.8 [vocalized-noise].
The cases we observed [vocalized-noise] its
giving least [vocalized-noise] [noise] spectral
radius [noise] [vocalized-noise] like for
[noise] Jacobi [noise] the spectral radius
[noise] was 0.98 [noise] [vocalized-noise]
Gauss-Siedel the spectral [noise] radius is
[noise] 0.97 [noise] [vocalized-noise] omega
1.8 [noise] spectral radius is 0.8 and [noise]
that is why [noise] convergence is faster
[vocalized-noise] [noise].
We look into the larger matrix [noise] 40
into 40 [noise] matrix [vocalized-noise] Gauss-Seidel
is [noise] 0.99 [vocalized-noise] [noise]
322 [noise] [vocalized-noise]. And ah [vocalized-noise]
Jacobi is [noise] 0.99661 [noise] the number
of iterations are [vocalized-noise] 2400 [noise]
or [vocalized-noise] [noise] 1291. [noise]
[vocalized-noise] Omega [noise] 1.8 [noise]
that [vocalized-noise] spectral radius is
[noise] 0.92 [vocalized-noise] and the [noise]
number of iteration [noise] is [noise] 1 6
2 [noise] is minimum [vocalized-noise].
So, again if we look [noise] varying omega
for a [vocalized-noise] 10 by 10 matrix [noise]
is the least [noise] comes at [vocalized-noise]
[noise] omega [vocalized-noise] ah [noise]
1.5 [noise] [vocalized-noise] for ah 20 by
20 the least value comes at [noise] omega
83 [vocalized-noise]. Here at [vocalized-noise]
omega 1 [noise] [vocalized-noise] omega 1.8
[noise] again [vocalized-noise] [noise]. So,
[vocalized-noise] there is an optimum [vocalized-noise]
SOR factor [noise] [vocalized-noise] for which
is [noise] number of iterations are least.
[noise] [vocalized-noise] And this we can
calculate using [noise] the [noise] formula
we have given earlier. [noise] [vocalized-noise]
That for [vocalized-noise] 10 [noise] it is
1.51 [noise] where we will get the least [vocalized-noise]
spectral radius for [vocalized-noise] 20,
[noise] it is 1.7; [noise] [vocalized-noise]
for 40 it is 1.84 [noise]. And these are the
values [noise] for which [noise] we get [noise]
least spectral radius [noise]. After that
[noise] [vocalized-noise] the spectral radius
[noise] will increase as well as number of
iterations will increase. [noise] [vocalized-noise]
However, [vocalized-noise] we can reduce the
spectral radius significantly [vocalized-noise]
using successive [noise] over relaxation [vocalized-noise]
and that is how [vocalized-noise] we can increase
[vocalized-noise] reduce the number of iterations
also [vocalized-noise] substantially. [noise]
[vocalized-noise]
So, [noise] we try to [vocalized-noise] cover
[noise] [vocalized-noise] few important basic
iterative [noise] methods [noise]. We will
look into [vocalized-noise] [noise] more complex
iterative methods [noise] which gives [vocalized-noise]
faster solution [noise]. And also see in the
[vocalized-noise] iterative methods [noise]
which [noise] work for non diagonally [noise]
dominant matrices [noise] in later classes
[noise]. Also [vocalized-noise] another important
thing which we will look [vocalized-noise]
is [noise] how to write [noise] computer programs
[noise] using this iterative [noise] methods
[noise] [vocalized-noise] the problems [noise]
actually [noise] which have been used for
the [noise] numerical experiments [noise]
in this purpose how to [noise] [vocalized-noise]
how can we [noise] reproduce these [noise]
programs [vocalized-noise] through these [noise]
methods. [noise]
Thank you. [noise] [vocalized-noise]
