Numerical Algorithms for Visual

Computing II

Assignment 9

Jan Hendrik Dithmar

2031259

Andreas Steinel

2029984

Assignment 9  Dithmar (2031259), Steinel (2029984)

9.1 Getting Symmetric

We have to derive symmetric versions of Gauss-Seidel (Part a) and SOR (Part b):

a) Symmetric Gauss Seidel (A = L + D + U):

x^{m+1} = −(D + U)^{−1}L [\underbrace{−(D + L)^{−1}Ux^{m} + (D + L)^{−1}b}_{ordinary Gauss-Seidel}] + (D + U)^{−1}b
= (D + U)^{−1}L(D + L)^{−1}Ux^{m} − (D + U)^{−1}L(D + L)^{−1}b + (D + U)^{−1}b
= \underbrace{(D + U)^{−1}L(D + L)^{−1}U}_{M_{SGS}} x^{m} + \underbrace{((D + U)^{−1} − (D + U)^{−1}L(D + L)^{−1})}_{N_{SGS}} b

For the case that A spd, this boils down to:

x^{m+1} = (D + L)^{−1}L(D + L)^{−1}Lx^{m} + ((D + L)^{−1} − (D + L)^{−1}L(D + L)^{−1}) b
= ((D + L)^{−1}L)^{2} x_{m} + ((D + L)^{−1} − (D + L)^{−1}L(D + L)^{−1}) b

b) Symmetric SOR:

x^{m+1} = (D + ωU)^{−1}[(1 − ω)D − ωL] (\underbrace{(D + ωL)^{−1}[(1 − ω)D − ωU]x^{m} + ω(D + ωL)^{−1}b}_{ordinary SOR}) + ω(D + ωU)^{−1}b
= \underbrace{(D + ωU)^{−1}[(1 − ω)D − ωL](D + ωL)^{−1}[(1 − ω)D − ωU]}_{MSSOR} x^{m} + \underbrace{((D + ωU)^{−1}[(1 − ω)D − ωL]ω(D + ωL)^{−1} + ω(D + ωU)^{−1})}_{NSSOR} b

As a further proof, lets assume ω = 1:

x^{m+1} = \underbrace{(D + U)^{−1}L(D + L)^{−1}U}_{MSSOR,ω=1} x^{m} + \underbrace{((D + U)^{−1}(−L)(D + L)^{−1} + (D + U)^{−1})}_{NSSOR,ω=1} b
= M_{SGS}x^{m} + N_{SGS}b

So, we see that for ω = 1, the SSOR method boils down to the SGS method.

9.2 Preposterous in Hilbert’s space

In this exercise, we investigate the Hilbert Matrix with corresponding Hilber vector as given by:

a_{i,j} = \frac{1}{i + j − 1}
H = (a_{i,j})_{i,j=1}^{n}
h_{i} = ∑_{j=1}^{n} a_{i,j}

This problem is a standard problem for numerical analysis which can shows the limitations of numerical methods.

a. We have to compute the real (in the mathematical sense) solution of

Hx = h

which is for row i defined as follows:

∑_{j=1}^{n} a_{i,j} x_{i} = ∑_{j=1}^{n} a_{i,j}

So, one can see that the solution is fullfilled by x_{i} = 1. Therefore, x = 1.

b. One can easily see that the trace of the matrix is not significantly increasing by the addition of one dimension, e.g.,

tr(H_{n+1}) = ∑_{i=1}^{n+1} a_{ii}
= tr(H_{n}) + \frac{1}{2 (n + 1) − 1}

For n → ∞, this series converges, which can easily be shown with the ratio test (it is always positive, so we omit the absolute values):

\frac{a_{n+1}}{a_{n}} = \frac{\frac{1}{2(n+1)−1}}{\frac{1}{2n−1}}
= \frac{2n − 1}{2n + 1} < C

Therefore, tr(H) = ∑_{i=1}^{n} λ_{i} is bound to C which is also less than 1 and therefore it converges absolutely.

So, we know that the trace represents the sum of all eigenvalues. In worst case, we have the same eigenvalues and the smallest eigenvalue is therefore ^{\frac{∑_{i=1}^{n} a_{ii}}{n}}, which is also decreasing by increasing n (n ≥ a_{ii}∀i. If we have different eigenvalues, the smallest eigenvalue is even smaller and also becomes even smaller.

This leads to the fact that the smallest eigenvalue becomes even smaller with increasing dimension n. In addition, the condition number increases tremendously:

dimension       condition number
 2 ->                        19
 3 ->                       524
 5 ->                    476607
10 ->            16024717801185
15 ->        398078708382037888
25 ->       5774300303076995072
50 ->      92977564954409222144
100 ->     141035992700580937728

Due to this big condition numbers, an interative solver have problems in determining the true error (compare to chapter 17, especially Theorem 17.4).

c. The CG-error in the maximum norm in comparison to the correct solution 1:

dimension      absolute error
 2 -> 0.00000000000000288658
 3 -> 0.00000000000069633188
 5 -> 0.00356480305521222895
10 -> 0.00008775740479349814
15 -> 0.00007187035314881740
25 -> 0.00002137496845033837
50 -> 0.00001437576758900327
100 -> 0.00049619274308021488

d. The Preconditioned CG-error has been computed as

dimension      absolute error
 2 -> 0.00000000000000011102
 3 -> 0.00000000000003508305
 5 -> 0.00000000438746028397
10 -> 0.05328292484288166264
15 -> 0.05922950956662620392
25 -> 0.04165658976247055278
50 -> 0.02999276484401813825
100 -> 0.02835902228564224359