I have been trying to port and understand the code used in the paper Real-Time Fluid Dynamics for Games. One function in the paper, the method to diffuse one of the vector fields is said to use the Gauss-Siedel relaxation to approximate the solutions.
Here is the function:
void diffuse ( int N, int b, float * x, float * x0, float diff, float dt )
{
int i, j, k;
float a=dt*diff*N*N;
for ( k=0 ; k<20 ; k++ ) {
for ( i=1 ; i<=N ; i++ ) {
for ( j=1 ; j<=N ; j++ ) {
x[IX(i,j)] = (x0[IX(i,j)] + a*(x[IX(i-1,j)]+x[IX(i+1,j)]+
x[IX(i,j-1)]+x[IX(i,j+1)]))/(1+4*a);
}
}
set_bnd ( N, b, x );
}
}
I understand what the set_bnd(N, b, x)
call is doing (handling interactions with the edges of the simulation), and I know that the outermost for-loop is doing the iteration 20 times, but how does this relate to the formula for the gauss-siedel method? Also, what is the term a
doing here, and why is the sum of the surrounding cells in the matrix multiplied by $\frac{a}{1+4a}$?