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
u_{i - 1} - 2u_{i} + u_{i + 1} =
h^{2}f_{i} ,
u_{0} = u_{n} = 0, h = 1/n,
x_{i} = i h, i = 0, 1, ... , n.
We can solve these equations with respect to u_{i} for Jacobi relaxation [1]
u_{i}^{new} =
(u_{i-1}^{old} + u_{i+1}^{old} - h^{2}f_{i} )/2.
Jacobi iterations converge to solution of the problem.
If we use weighted value u_{i}^{old} +
ω(u_{i}^{new} - u_{i}^{old})
we get weighted Jacobi relaxation (we put f_{i} = 0 further)
u_{i}^{new} = (1 - ω)u_{i} +
ω(u_{i-1} + u_{i+1} )/2, (*)
Press the Step button below to see relaxation process for f_{i} = 0,
n = 64 and initial condition
u_{i} = [c1 sin(πi/n) + c3 sin(3πi/n) + c5 sin(5πi/n)]/3.
it
c1
c3
c5
ω
Ymin
Ymax

Error log( max |u_{i }| ) 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
u_{i}^{(m)} = sin πmi/n = sin k_{m }i is eigenvector of (*) with the eigenvalue
λ_{m} = (1 - ω) + ω cos k_{m}
= 1 - 2ω sin^{2}π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")