Introduction
It looks like I will be quite busy for a while so I won’t have much time for the blog. I might only post easy things like this one for a while… anyway, here it is, with only brief explanations and the JavaScript code. I didn’t feel the need and didn’t have the time and patience to write a C++ program for this topic. It’s just too simple and besides, it’s nice to see the code action in the page.
Short Description
Links
So, this post is about the relaxation method combined with a multigrid method applied on the Laplace equation.
The methods are easy and I suppose Wikipedia pages are a good start, but here there are several other links for help, I just googled them, you can find many more on the subject. First, a paper. A web page is here and another one here.
Very short theory description
You can find the Laplace equation (or more general, the Poisson equation) in various topics, for example originating from Gauss law if you use the electric potential instead of the electric field, or in the heat equation in a steady state when the time derivative drops out (also see: diffusion equation).
What is nice about this equation, apart from the superposition principle is that the solutions are harmonic functions which means that the value in a point is the average of the points around it (for a more rigorous explanation, please see the Wikipedia page). This allows us to use a relaxation method to solve it, although there are faster methods that one could find (for example, one could use Fourier transform). Despite this, the method is easy to understand and can be a start for other methods of solving the equation.
The method is very easy, first the equation is discretized using a finite difference method then one iteratively averages the points in the discretized space. This can be coupled with a multigrid method by starting out with a coarse grid first.
Code in action
That should be enough theory, here is the result in action, on a very simple model on a square with the boundary condition of two sides with the field with value 1 and the other two with value -1.
The code just iterates relaxation until the difference between the old solution and the new one is under a certain threshold, then the resolution is increased and the same is done again until it reaches a certain resolution, then it starts again.
It could be instructive to change the code to not use the multigrid method, but instead to start from the beginning with the smallest mesh size.
The Code
Here is the model class:
function SquareModel(Size)
{
this.Size = Size;
this.field = [];
this.IndexForSize = function(row, col, size) { return size*row+col; };
this.Index = function(row, col) { return this.IndexForSize(row, col, this.Size); };
this.Value = function(matrix, row, col) { return matrix[this.Index(row, col)]; };
this.Boundary = function(row, col) { return row == 0 || col == 0 || row == this.Size - 1 || col == this.Size - 1; };
this.Field = function(row, col) { return this.field[this.Index(row, col)]; };
this.SetBoundary = function() {
for (i = 1; i < this.Size - 1; ++i)
{
this.field[this.Index(0, i)] = -1;
this.field[this.Index(this.Size - 1, i)] = -1;
this.field[this.Index(i, this.Size - 1)] = 1;
this.field[this.Index(i, 0)] = 1;
}
};
this.Init = function() {
for (i = 0; i < this.Size; ++i)
for (j = 0; j < this.Size; ++j)
this.field[this.Index(i, j)] = 0;
this.SetBoundary();
}
this.Init();
}
There is not much to say about it apart of what I already described. I guess you could change it with some other model…
Here is the Relaxation class:
function Relaxation(relaxationModel)
{
this.StartSize = relaxationModel.Size;
this.iteration = 1;
this.newField = [];
this.model = relaxationModel;
this.SwapFields = function() {
tmp = this.model.field;
this.model.field = this.newField;
this.newField = tmp;
};
this.MakeGridSmaller = function() {
for (i = 0 ; i < this.model.Size; ++i)
for (j = 0 ; j < this.model.Size; ++j)
{
this.newField[this.model.IndexForSize(2*i, 2*j, this.model.Size*2)] = this.model.Field(i, j);
this.newField[this.model.IndexForSize(2*i+1, 2*j, this.model.Size*2)] = (this.model.Field(i, j) + (i<this.model.Size - 1 ? this.model.Field(i+1, j) : this.model.Field(i, j)))/2;
this.newField[this.model.IndexForSize(2*i, 2*j+1, this.model.Size*2)] = (this.model.Field(i, j) + (j<this.model.Size - 1 ? this.model.Field(i, j+1) : this.model.Field(i, j)))/2;
this.newField[this.model.IndexForSize(2*i+1, 2*j+1, this.model.Size*2)] = (this.model.Field(i, j) + (i<this.model.Size - 1 ? this.model.Field(i+1, j) : this.model.Field(i, j)) +
(j<this.model.Size - 1 ? this.model.Field(i, j+1) : this.model.Field(i, j)) +
(i<this.model.Size - 1 && j<this.model.Size - 1 ? this.model.Field(i+1, j+1) : this.model.Field(i, j)))/4;
}
this.SwapFields();
this.model.Size *= 2;
this.model.SetBoundary();
};
this.Reset = function() {
this.model.Size = this.StartSize;
this.model.Init();
this.iteration = 1;
};
this.Relax = function() {
change = 0;
for (i = 0 ; i < this.model.Size; ++i)
for (j = 0 ; j < this.model.Size; ++j)
if (!this.model.Boundary(i, j))
{
oldVal = this.model.Field(i, j);
newVal = (this.model.Field(i - 1, j) + this.model.Field(i, j - 1) + this.model.Field(i, j + 1) + this.model.Field(i + 1, j))/4.0;
this.model.field[this.model.Index(i, j)] = newVal;
dif = oldVal - newVal;
change += dif * dif;
}
++this.iteration;
return Math.sqrt(change);
};
}
The code uses the Gauss-Seidel method instead of the slower Jacobi method. As a note, the interpolation method I made in the MakeGridSmaller method is very crude, it shifts the values up and to the left than they should be, one can do much better and it should be done better in ‘real life’ code. I didn’t have patience to implement something better in this short time I allocated for this post.
In order to be fully functional, the code needs also functions to create the objects, successively apply the methods and display the results, and here it is:
function DisplayModel(canvas, model) {
function Color(val) {
r = 0;
g = 0;
b = 0;
if (val > 0)
{
r = Math.ceil(255. * val);
g = Math.floor(255. * (1. - val));
}
else
{
val *= -1;
g = Math.floor(255. * (1. - val));
b = Math.ceil(255. * val);
}
return "rgb(" + r.toString() + "," + g.toString() + "," + b.toString() + ")";
};
ctx = canvas.getContext("2d");
displaySize = canvas.width / model.Size;
rectSize = Math.ceil(displaySize);
for (i = 0 ; i < model.Size; ++i)
for (j = 0 ; j < model.Size; ++j)
{
ctx.fillStyle = Color(model.Field(i, j));
ctx.fillRect(j * displaySize, i * displaySize, rectSize, rectSize);
}
};
ParticularModel1 = new SquareModel(8);
Relaxation1 = new Relaxation(ParticularModel1);
canvas = document.getElementById("relaxationCanvas");
function Tick() {
DisplayModel(canvas, Relaxation1.model);
error = Relaxation1.Relax();
if (error / (Relaxation1.model.Size * Relaxation1.model.Size) < 0.0000001)
{
Relaxation1.MakeGridSmaller();
if (Relaxation1.model.Size == 256) Relaxation1.Reset();
}
}
setInterval(Tick, 10);
Everything is packed into an anonymous closure and that’s about it:
(function() {
})();
Conclusion
Maybe I was too brief on this subject but it’s weekend… as usual, if you find something wrong, please point it out.
There are many things one might try to improve such algorithm, for example the interpolation when switching to a finer mesh could be improved and using over-relaxation will speed up the convergence.
The post Relaxation Method first appeared on Computational Physics.
