Relaxation Method

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)
       {        
           // two borders to -1, two 1           
           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();
}

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)

{

// two borders to -1, two 1

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() {
       // make a smaller grid by interpolating the values from the source
       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;

             // the above is not quite correct for the rightmost column and last row, but it doesn't matter, SetBoundary should set them to the correct values, anyway
          }

      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);
   };
}

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() {

// make a smaller grid by interpolating the values from the source

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;

// the above is not quite correct for the rightmost column and last row, but it doesn't matter, SetBoundary should set them to the correct values, anyway

}

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);

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() {

// all the code from above is here

})();

(function() {

// all the code from above is here

})();

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.

CodeProject

Relaxation Methods for Partial Diﬀerential Equations: Applications to Electrostatics by David G. Robertson ↩
Poisson’s Equation and Relaxation Methods part of a Computational Physics lecture ↩
Partial differential equations Numerical Methods, Natural Sciences Tripos 1B by Stuart Dalziel ↩

Introduction

Short Description

Links

Very short theory description

Code in action

The Code

Conclusion

Like this:

Related

Leave a Reply (for privacy concerns, please read the Privacy Policy Page accessible from the top menu)Cancel reply

Introduction

Short Description

Links

Very short theory description

Code in action

The Code

Conclusion

Share this:

Like this:

Related

Leave a Reply (for privacy concerns, please read the Privacy Policy Page accessible from the top menu)Cancel reply