Jacobi method demo

Code below is to solve the static version of the equation in previous FDM code demo, with B.C. changed as homogeneous. The Jacobi method is used to replace the LAPACK solver which gives inexact method. The exact solution is

Using initial guess as

one can see that most of the computing cost is spent on smoothing lowest mode.

#include <stdio.h>
#include <math.h>
#include <iostream>
#include <fstream>
//#include <vector>
//#include <sunperf.h>
using namespace std;
//typedef std::vector<float> fvector_s;

//Maximum of an array.
float max(float a[], int n) {
 float v = a[0];
 for (int i=1; i<n; i++) {
 if (a[i] > v) {
 v = a[i];
 }
 }
 return v;
}

int main ()
{
 int n=500;
 float tmp, b[n], a_d[n], a_u[n-1], a_l[n-1],    \
 tol=1.e-7, x0[n], x1[n], x2[n], error=1.;
 float h=1./(n+1), pi=3.1416;
 int count=0;
 
 //initializing a
 for(int i=0; i<n; ++i){
 a_d[i]=2.0;
 }
 for(int i=0; i<n-1; i++){
 a_l[i]=-1.;
 a_u[i]=-1.;
 }
  
 //initializaing b[]
 for(int i=0; i<n; ++i){
 b[i]=h*h;
 }

 //initial guess
 for(int i=0; i<n; ++i){
 x0[i]=sin(2*pi*h*(i+1))+0.5*sin(10*pi*h*(i+1))+\
       0.3*sin(20*pi*h*(i+1))+0.1*sin(40*pi*h*(i+1));
 x1[i]=x0[i];
 }
 
 ofstream outfile;
 outfile.open("jac.out");
 do{
 count+=1;
 for(int i=0; i<n; ++i){
 if(i==0)
 x2[i]=(b[i]-a_u[i]*x1[i+1])/a_d[i];
 else if(i==n-1)
 x2[i]=(b[i]-a_l[i-1]*x1[i-1])/a_d[i];
 else
 x2[i]=(b[i]-a_u[i]*x1[i+1]-a_l[i-1]*x1[i-1])/a_d[i];
 }
 for(int i=0; i<n; i++){
 x1[i]=x2[i]-x1[i];
 }
 error=max(x1, n)/max(x0, n);
 for(int i=0; i<n; i++){
 x1[i]=fabs(x2[i]-(-(i*h)*(i*h)/2+i*h/2));
 }
 outfile << count << "\t" << error << endl;
 for(int i=0; i<n; ++i){
 x1[i]=x2[i];
 }
 }while(error>tol);
 cout << "Number of iterations:" << count << endl;
 for(int i=0; i<n; i++){
 outfile << count << "\t" << x2[i] << endl;
 }
 
 outfile.close();
 return 0;
}

Result of Jacobi method at different iteration

Alpha shape

Alpha shape is a generalized way determining the geometric objects by a set of points, of which convex hull is a special case. It can be applied to FEM for boundary recognition. Essentially this method identifies the boundary of a set of points by eliminating those triangles (generated by triangulation of the points) which can be put into empty circles whose diagram are defined by threshold alpha. It is intuitive for a more simple criteria: eliminating the triangles whose circumcircles contains no points in the defined set. Below is a simple demo by Mathematica which tries to identify the boundary of a subset of points within circle radius 5 from a set points randomly distributed. The red arc isolates the point set   inside, the first step is to eliminate the triangles whose circumcirlces contains no points in , the next step is to eliminate the triangles contain nodes other than those in . More delicate results can be reached by tuning alpha circle instead of using simple circumcirlces.

Bowyer/Watson algorithm: 3D

3D version implementation of Bowyer/Watson algorithm is almost the same as the 2D. Instead of eliminating lines that belong only to the triangles whose circumcirlce contains the newly added point, in 3D we eliminate the facets that belong only to the tetrahedra whose circumsphere contains the newly added point, then form new tetrahedra by appending new point to those facets’ nodes. Here is an example generated by Mathematica scripts I compiled. Below on the right is tetrahedronization of 200 random points, together with enclosing tetrahedron. On the left is he tetrahedrons generated by first two points.

DT of 200 points

2 points

DistMesh

DistMesh is no doubt the most interesting way meshing a domain, simple and elegant. There is a MATHEMATICA version of the original Matlab code, but only in 2D. Base on that, I wrote a script for mesh test. Hopefully continue working on this one would bring some application to my research. Here is one simple case of generating a non-uniform mesh within a circle, whose distance funtion is the easist to define. At the end of termination, most of the ill-shaped triangles are corrected.

Classical iterative algorithms

For this topic, many excellent textbooks are available. Most of the modern numerical PDE solvers rely on one or another iterative algorithm taking care of the linear system given by discretization ( I know very little of direct methods). In this scenario, one has linear system

and tries to come up with some iterative scheme

with

meaning the convergence to the “real” solution, or, to be exact, the solution by direct solver if we disregard the round-off error. By this setup, one can see the solution is the fixed point of operator . So the convergence of the algorithm is equivalent to two proposals: the spectral radius , and there exists operator norm induced by some norm of the linear space to which belongs that this norm of is less than . The convergence speed is measured by average rate of convergence of (for iterations) defined by

which comes from the fact that . In order to put the original problem into this fixed point setup, is splitted:

If is splitted according to entries’ positions as , one has the classical algorithms:

Jacobi:

Gauss-Seidel:

SOR:

The last one is in the component wise form in which is the Gauss-Seidel result, and the relaxation parameter is between zero and two (even though the name “overrelaxation” comes when it is greater than one). Each of those algorithms has direct(but different in implementation) analogical block version. All these methods are stationary, i.e. in each iteration is unchanged. If this is not the case, we have unstationary methods, in which after each iteration information is re-gathered in order to determine a new for next iteration. When is a SPD matrix, a large family of unstationary methods can be formed by the solving a minimization problem(minimal residual method (MINRES)). Classical CG method I wrote about before belongs to this family. A modern member of this family is GMRES by Saad and Schultz in 1986.