Solving Linear Systems
Linear systems arise in many practical applications and solving
linear systems is a key issue
in modern scientific computing.
In this context, we will focus on solving a tridiagonal linear system
arising from solving
numerical Ordinary Differential Equation (ODE).
ODE Example
We consider the following second order Ordinary Differential Equation
(ODE)
f''(x)=x, f(0) = f(1) = 0, x Î [0,1].
It can be shown easily that the solution of the above ODE is
f(x) = |
(x - 1) x (x + 1) ¾¾¾¾¾¾¾
6 |
Remark: Closed form solution is not generally available and
hence numerical methods are
popular tools for solving this kind of problem.
Here we consider solving the above ODE numerically.
Discretization of ODE
One may first divide the interval [0,1] into n equal parts
(discretization).
We then let
xi = |
i ¾ n |
for i = 0, 1, 2, ..., n |
We note that f(x0) = f(0) = f(xn) = f(1) = 0.
We then adopt the following approximation of f''(x) at the points xi:
f''(xi) = |
f(xi-1) - 2 f(xi) +
f(xi+1)
¾¾¾¾¾¾¾¾¾
1 / n2 |
for i = 1, 2, ..., n - 1. |
Remark: In this case we have n-1 interior grid points in (0,1).
From above we have
Af = |
æ ç ç ç ç è |
2 | -1 | 0 | 0 | 0 |
ö ÷ ÷ ÷ ÷ ø |
-1 | 2 | -1 | .. | : |
0 | -1 | : | : | 0 |
: | : | : | 2 | -1 |
0 | .. | .. | -1 | 2 |
|
æ ç ç ç ç è |
f(x1) |
ö ÷ ÷ ÷ ÷ ø |
f(x2) |
: |
: |
f(xn-1) |
|
= |
1 ¾ n2 |
æ ç ç ç ç è |
f"(x1) |
ö ÷ ÷ ÷ ÷ ø |
f"(x2) |
: |
: |
f"(xn-1) |
|
= b. |
By solving the above tridiagonal linear system, one may get an
approximation of the solution of
ODE at the points xi Î [0,1].
Remark: It is expected that as n ® ¥,
the approximation will tend to
the true solution.
But we will not discuss this issue in this context. Interested reader
can consult the book listed in References section.
There are many techniques for solving a linear system of equations.
Generally speaking, two classes of methods: direct methods and
iterative methods.
One example of direct methods is the LU factorization.
Given a non-singular linear system of equations Af = b,
one would like to factorize A = LU, where L is a lower triangular
matrix and U
is a upper triangular matrix.
If we can do so, the solution of the linear system Af = b can be
obtained by first solving Ly = b by forward substitution
and
then solving Lf = y by backward substitution.
We will discuss the Doolittle's LU factorization shortly.
In an iteration method, a sequence of approximate solutions {xk}
of Af = b is generated such that
For iterative methods, one aims at saving computational cost by
obtaining an approximate
solution.
This means that we will terminate the iteration when certain criteria
is met.
Stopping Criteria 1: The iteration is terminated when
||Axk - b|| £ t
where t is a pre-determined tolerance (e.g. t = 0.000001) and
||.|| is vector norm.
One possible choice of the vector norm is ||.|| is the Euclidean
norm
||x||2 = ( |
n S i=1 |
vi2 )1/2 |
Stopping Criteria 2: Another popular criteria is that the
iteration is terminated when
||xk - xk-1|| £ t.
A class of iterative methods can be derived by using the idea of
matrix splitting.
We will discuss this idea in detail.
|