Large, sparse linear systems of equations Ax = b arise in many applications e.g.
FEM of linear elasticity
Fluid dynamics
Schrodinger equation
Many approaches can be used to find x on GPU
FFT
Conjugate gradients
Jacobi relaxation
Multigrid algorithm
A few Jacoby iterations are enough to get good looking visual effect in
2D fluid simulations but long wave relaxation is very slow on fine
grid. Therefore we need multigrid algorithm for accurate simulation of heart tissue contraction.
Yang Chenglin sent me kindly link on his
FFT-based Droplet Simulation.
Conjugate gradients solver will be very interesting too.
1D weighted Jacobi relaxation algorithm
Consider one-dimensional boundary value problem
u" = f(x), 0 < x < 1, u(0) = u(1) = 0,
with the finite difference we get
ui - 1 - 2ui + ui + 1 =
h2fi ,
u0 = un = 0, h = 1/n,
xi = i h, i = 0, 1, ... , n.
We can solve these equations with respect to ui for Jacobi relaxation [1]
uinew =
(ui-1old + ui+1old - h2fi )/2.
Jacobi iterations converge to solution of the problem.
If we use weighted value uiold +
ω(uinew - uiold)
we get weighted Jacobi relaxation (we put fi = 0 further)
uinew = (1 - ω)ui +
ω(ui-1 + ui+1 )/2, (*)
Press the Step button below to see relaxation process for fi = 0,
n = 64 and initial condition
ui = [c1 sin(πi/n) + c3 sin(3πi/n) + c5 sin(5πi/n)]/3.
it
c1
c3
c5
ω
Ymin
Ymax
Error log( max |ui | ) is plotted against iteration number (the red curve).
Note that at first fast oscillations decay quickly, then the main harmonic decays slowly at
constant rate.
It is not difficult to check by direct substitution that
ui(m) = sin πmi/n = sin km i is eigenvector of (*) with the eigenvalue
λm = (1 - ω) + ω cos km
= 1 - 2ω sin2πm/2n.
λm determines relaxation rate for weighted Jacobi iterations.
For ω = 1/3, 1/2, 2/3, 1 they are plotted below
(-1 ≤ λm ≤ 1 and 0 ≤ m ≤ n, see also [1]).
For ω = 2/3 (the red curve) and n/2 < m < n we get
|λm| ≤ 1/3 therefore high harmonics decay at least 3 times after every
iteration.
For small m/n we get approximately
λm = 1 - ω(πm/n)2/2,
so relaxation is very slow for m/n < 1/100 values. As since convergence is faster
for small n we need caurse grids to accelerate relaxation. We consider
2D Poisson equation next.
[1] William L. Briggs, Van Emden Henson, Steve F. McCormick
A Multigrid Tutorial
(look at "Tutorial Slides")