連立方程式のための反復法（主にヤコビ法）を、差分法（主にラプラス方程式）に適用する場合

Ax=bなる行列形式の連立方程式を解くことを考える。
求めるべきは未知なるxベクトルであり、係数行列Aとbベクトルは既知。

反復法は未知なるxベクトルに適当な加工を施し、出力たるx'ベクトルにまた同じ加工を施す。xとx'にほとんど変化が見られなくなったら収束とみなす。気を付けるべきは……

・ちゃんと収束するか
・一番最初のxベクトルは何にすべきか
・収束をどの程度どのように判定するか

ちゃんと収束するか

・行列の固有値の絶対値が全て１より小さければ収束する。ただし固有値を調べるのは結構大変である。下手したら反復法が収束するかどうか調べるための固有値を求めるために反復法を使用することになる。
・対角優位ならば収束する。対角優位であるとは、行列の対角要素（ピボット要素）の絶対値が、同じ列の他の要素の絶対値を全部足したものより大であること。
・差分法や有限要素法の文脈ではとりあえず試してみれば良い。

一番最初のxベクトルは何にすべきか

・基本的にはどのようなベクトルから始めても良いが、手法によっては最適な初期値ベクトルが存在したり、適切でない初期値ベクトルが存在したりする。つまり場合による。

収束をどの程度どのように判定するか

・L2ノルム：xベクトルとx'ベクトルの長さの差を判定する
・L1ノルム：xベクトルとx'ベクトルの成分の絶対値の総和の差を見る
・xベクトルとx'ベクトルの成分の2乗和の差を見る
などなど。
また、収束の判定は反復1回毎に行う必要はない。別に3回に1回でも10回に1回でも良い。

ヤコビ法

固有値を求めるための方法に同名のヤコビ法があるが、それとは別。
以下はWikipediaに乗ってるアルゴリズムをC#で書いたもの

参考文献
山本 哲朗 (著)
数値解析入門 (サイエンスライブラリ現代数学への入門)
森 正武 (著), 古屋 茂 (編集), 赤 摂也 (編集), 一松 信 (編集)
数値解析 (共立数学講座)

        public static double[] liner_Jacobi(double[,] A, double[] b, double[] x, double e = 0.0001)
       {
           int A_row_count = A.GetLength(0);//row
           int A_col_count = A.GetLength(1);//col

           double[] x_next = new double[x.Length];

           for (int row_index = 0; row_index < A_row_count; row_index++)
           {
               double pivot = A[row_index, row_index];

               double bt = b[row_index];
               for (int col_index = 0; col_index < A_col_count; col_index++)
               {
                   if (row_index == col_index) { continue; }
                   bt -= A[row_index, col_index] * x[col_index];
               }
               x_next[row_index] = bt / pivot;
           }
           
           //収束判定
           if (Math.Abs(norm2(x_next) - norm2(x)) < e)
           {
               return x_next;
           }
           else
           {
               return liner_Jacobi(A, b, x_next, e);
           }
       }

行列形式

        public static double[] liner_Jacobi_LDU(double[,] A, double[] b, double[] x, double e = 0.0001)
       {
           int A_row_count = A.GetLength(0);//row
           int A_col_count = A.GetLength(1);//col

           var D = Copy(A);
           var LU = Copy(A);

           //make D LU
           for (int row_index = 0; row_index < A_row_count; row_index++)
           {
               for (int col_index = 0; col_index < A_col_count; col_index++)
               {
                   if (row_index == col_index)
                   {
                       //D-1
                       D[row_index, col_index] = 1 / A[row_index, col_index];
                   }
                   else
                   {
                       LU[row_index, col_index] = A[row_index, col_index];
                   }
               }
           }

           var x_next = MatrixMultiple(D, Sub(b, MatrixMultiple(LU, x)));

           //収束判定
           if (Math.Abs(norm2(x_next) - norm2(x)) < e)
           {
               return x_next;
           }
           else
           {
               return liner_Jacobi_LDU(A, b, x_next, e);
           }
       }

T.R.マッカーラ (著), 三浦 功 (著), 田尾 陽一 (著)
計算機のための数値計算法概論 (サイエンスライブラリ情報電算機 8)
p153
Add()はdouble[]とdouble[]の成分ごとの和

        public static double[] liner_Jacobi_A0(double[,] A, double[] b, double[] x, double e = 0.0001)
       {
           int A_row_count = A.GetLength(0);//row
           int A_col_count = A.GetLength(1);//col

           var _A = Copy(A);
           var _b = Copy(b);

           for (int row_index = 0; row_index < A_row_count; row_index++)
           {
               double pivot = A[row_index, row_index];

               _b[row_index] /= pivot;
               for (int col_index = 0; col_index < A_col_count; col_index++)
               {
                   if (row_index == col_index)
                   {
                       _A[row_index, col_index] = 0;
                   }
                   else
                   {
                       _A[row_index, col_index] /= -pivot;
                   }
               }
           }

           var x_next = Add(MatrixMultiple(_A, x), _b);

           if (Math.Abs(norm2(x_next) - norm2(x)) < e)
           {
               return x_next;
           }
           else
           {
               return liner_Jacobi_A0(A, b, x_next, e);
           }
       }

差分法に適用する場合

ただし言語はprocessing

以下はラプラス方程式の反復計算の1回分。
Next[][]やPrev[][]は2次の配列であるが、これは2次元画像というか、2次元格子を走査しているからである。これは全ての成分が未知である1次のxベクトルに相当する。つまりxベクトルを加工してx'ベクトルを作っている。

    for(int i = 1; i<this.Prev.length-1; i++)
   {
     for(int j = 1; j<this.Prev[0].length-1; j++)
     {                
       //ラプラス方程式のヤコビ法
       this.Next[i][j]=(this.Prev[i-1][j]+this.Prev[i+1][j]+this.Prev[i][j-1]+this.Prev[i][j+1])/4;       
     }
   }    

どのみち成分を走査するので、ついでに収束を判定しても良い。

    //誤差初期化
   e_max = 0;
   
   //上下左右の端を走査しない
   for(int i = 1; i<this.Prev.length-1; i++)
   {
     for(int j = 1; j<this.Prev[0].length-1; j++)
     {                
       //ヤコビ法
       this.Next[i][j]=(this.Prev[i-1][j]+this.Prev[i+1][j]+this.Prev[i][j-1]+this.Prev[i][j+1])/4;
       
       //誤差
       double e = (this.Next[i][j]-this.Prev[i][j])*(this.Next[i][j]-this.Prev[i][j]);
       if(e_max<e)
       {
         e_max = e;
       }
     }
   }   

まとめると以下のような感じ。入力は境界のみ入力された行列。

  int iterative_max_count = 100; 
 int iterative_count = 0;
 double e_threshold = 0.01;
 double e_max = 0;  
 double[][] Prev;
 double[][] Next;

 void jacobi_method(double[][] matrix)
 {      
   iterative_count = 0;
   
   copy_matrix(matrix, this.Prev);       
   copy_matrix(matrix, this.Next); 
   
   this.jacobi_method_once();    
   while(e_threshold<e_max)
   {
     if(iterative_max_count<=iterative_count){break;}
     
     this.jacobi_method_once();
     iterative_count++;
   }
   
   //copy
   copy_matrix(this.Next, this.region_matrix);
 } 

 void jacobi_method_once()
 {    
   //誤差初期化
   e_max = 0;
   
   //上下左右の端を走査しない
   for(int i = 1; i<this.Prev.length-1; i++)
   {
     for(int j = 1; j<this.Prev[0].length-1; j++)
     {                
       //ヤコビ法
       this.Next[i][j]=(this.Prev[i-1][j]+this.Prev[i+1][j]+this.Prev[i][j-1]+this.Prev[i][j+1])/4;
       
       //誤差
       double e = (this.Next[i][j]-this.Prev[i][j])*(this.Next[i][j]-this.Prev[i][j]);
       if(e_max<e)
       {
         e_max = e;
       }
     }
   }    
   copy_matrix(this.Next, this.Prev);
 }

拡散方程式などと形式が似ているため、同じ要領でNext[][]を描画していけば反復法が収束していくさまが多分見れる（試してない）