Guide through MGC
Level of sophistication: this report should be understandable for students
with a batchelor's degree in the exact sciences.
Meanwhile this report, with its 64 pages, has reached the size of a booklet.
Contrary to most writeups of contemporary mathematics, its level of detail
is quite high. And the amount of referencing to sources, other than itself,
is quite low :-)
Newton-Rhapson Algorithm
The starting point of the report is the observation that a division can also
be performed by employing the Newton-Rhapson algorithm. Find the solution of
the equation 1/x = a . Giving the iterations: p(n+1) = 2.p(n) - a.p(n)^2 .
This can be generalized for matrices A , instead of numbers a . Let A = 1-M ,
then it can be proven for the n-th iterand that: P(n) = (1-M^m)(1-M)^(-1) ,
where m = 2^n . Thus establishing a relationship with Geometric Series.
The geometric series of matrices turns out to be equivalent with an
"incremental Jacobi" iterative solution method.
Another way of looking at the series is the so-called Euler expansion, which
involves products of (one plus matrix powers of two).
2 x 2 matrix
Application of the Euler expansion is demonstrated at hand of a 2 x 2
matrix. An exact solution is found almost immediately.
MultiGrid
At hand of an arbitrary one-dimensional space-like grid - corresponding with
a three-diagonal matrix - it is demonstrated that application of the Euler
expansion is equivalent with subsequent splitting of the grid in even and
odd-numbered courser grids. This in turn corresponds with a splitting of the
associated matrix in blocks; even and odd indices become uncoupled in this
manner. The notion of re-normalization is introduced. The procedure ends up
in the identity matrix. And the solution of the linear 1-D system is found.
Remark:
I don't know if this method is "true" MultiGrid, but the name is certainly not
a misnomer, Because multiple grids are actually involved. So I don't
mind either :-)
Direct Solver
The Newton Rhapson method with Euler expansion for the equation X^(-1) = A
is thus a Direct Solver. This is confirmed once more, because an equivalence
can be established with Gaussian elimination.
Thus, for the 1-D case, we have established the following equivalences:
Newton Rhapson - Geometric Series - (incremental Jacobi method) - Euler
Expansion - MultiGrid - Gaussian elimination
Persistent Properties
From now on, our one-dimensional grid will be restricted even more, namely
to a uniform grid. The notion of Persistent Properties is introduced:
properties which are independent of grid coursening or grid refinement. It is
found that the off-diagonal elements of the system matrix are only persistent
if they are negative or zero (: one of Patankar's Four Basic Rules). Moreover,
it is demonstrated that only a few types of matrices are persistent through
all multigrids. The sum of the off-diagonal elements being smaller, equal
or greater than one (a + b = 1, a + b < 1, a + b > 1) is a persistent
property too.
Quotient Function
The quotient of the off-diagonal matrix elements a/b is studied as a persistent
function of the grid spacing dx . The result is: a/b(dx) = exp(P.dx) . Here P
is an arbitrary constant (which is yet to be determined).
I felt not entirely comfortable with the "proof" in this section, though.
Conjecture: f(2.x) = f(x).f(x) , x continuous, if and only if f(x) = exp(P.x) .
Problem: prove or disprove the necessary part (sufficient being trivial).
This result has been achieved originally by sort of transfinite induction: it
is "proved" for grid spacings approaching zero - the domain of Calculus - and
then established for Numerical Analysis, while working backwards to finite
sized grid spacings. In terms of Chaos Theory: the theorem is proved for
the Analytical Attractor in a/b(0) = 1 and then very much extrapolated.
Therefore this section has been augmented by an alternative proof, where
results from 'Direct Solver' have been employed, applied to a non-
uniform 1-D mesh with a rather dedicated layout. This results in a slightly
different conjecture:
f(p+q) = f(p).f(q) if and only if f(x) = exp(P.x) .
The latter can be proved, though (much easier - if at all - than the former :-).
Some Stable Solutions
With the previous result at hand, it is not difficult to find all solutions of
the three-diagonal system, for the very special - but persistent - case that
a + b = 1 . Meaning that any temperature (i.e. arbitrary quantity) at the grid
is a weighted mean of the neighbouring quantities. Solutions are recognized as
those for 1-D Convection and Diffusion.
Product Function
The product of the off-diagonal matrix elements a.b is studied as a function
of doubling the grid spacing:
a.b := [ a.b / (2.a.b - 1) ]^2 .
While solving the
linear equations system, the denominator (2.a.b-1) is formed every time again.
Values of a.b which can become 1/2, in the end, are therefore called
dangerous.
But it's actually easier to work backwards, towards refined grids: starting with
the value a.b = 1/2 and then looking for values of a.b which can eventually
lead to a zero denominator.
It is established that the following quantity is persistent through multigrids:
1 - 4.a.b . This quantity is called a discriminant. It is shown that
the value a.b = 1/4 effectively splits the domain of interest into two distinct
pieces: a.b ≤ 1/4 and a.b ≥ 1/4 . Here the interval [1/4, ∞] is
identified as our dangerous interval. In contrast herewith, the interval [0,1/4]
should be identified as being safe: because it can be shown that these
values of a.b will never result in a zero denominator.
Trigonometric Connections
There exists most surprising relationships for the grid refinement function f ,
if it is written as f(dx) = 1/(a.b) . Grid coursening - that is doubling the
grid spacing - then results in f(2.dx) = [ f(dx) - 2 ]^2 . This is very much
analogous with doubling an angle with the cosine function: g(x) = 2 + cos(2.x) .
Because with 2 + 2.cos(2.x) = 2.(2.cos^2(x) - 1) + 2 = 4.cos^2(x) =
([2 + 2.cos(x)] - 2)^2 , we find: g(2.x) = [ g(x) - 2 ]^2 .
Thus we are
tempted to identify the functions f and g . The more precise argument in the
report reveals that the following formula represents all subsequent iterations
in the dangerous domain:
1/(a.b) = 2 + 2.cos(k.pi/2^N) , where k = 0,1,2, .. ,2^N .
A multitude of conclusions can be drawn by employing the above formula.
It is demonstrated, for example, that solving linear equations in the
dangerous domain is comparable with unstable equilibrium in physics.
It leads to unpredictable, chaotic (!) behaviour of the iterands,
which nevertheless can be described with high precision.
Apart from this link to Chaos Theory, there is another interesting
link to binary-reflected Gray-codes.
The Hyperbolic Connection
Formulas for doubling and halving an "angle" are not only obeyed by the
trigonometric cosine function. They are also valid for its counterpart, which
is associated with an orthogonal hyperbola instead of a circle, the so-called
hyperbolic cosine. Let g(x) = 2 + 2.cosh(x) , then the following formula is
found for the safe domain of interest: 1/(a.b) = 2 + 2.cosh(Q.dx) . Which is
thus valid for 0 ≤ a.b ≤ 1/4 . With this hyperbolic connection, it seems
that everything fits in place. Explicit formulas for the off-diagonal matrix
elements are found. And the exact solution of any uniform tri-diagonal system
can be written down.
Two 'sci.math' posters
give a more concise overview of the above: |
|
Governing Equation
Imagine an immensely refined uniform grid. For such a grid, the solution of our
linear equations, corresponding with a matrix of infinite size, should become
exactly
the same as the outcome of an analytical model. The latter must be found with
help of calculus. It turns out that the linear equations, indeed, are
equivalent with a second order linear differential equation (ODE):
A.T'' + B.T' + C.T = 0 . Where T = unknown (temperature) , A,B,C = constants.
Moreover, the constants A,B,C must be such that the discriminant of the
characteristic (quadratic) equation is positive or zero. And there
is no danger to confuse this discriminant with the one defined for the linear
system of equations, because both discriminants turn out to be the "same".
The rest of the report consists of the investigation of special cases
which deserve some further attention, a re-statement of the main results,
a summary of Persistent Properties and a couple of Appendices.