D言語で連立一次方程式を解く(ヤコビの反復法)

例として
$9x+2y+z+u=20$
$2x+8y-2z+u=16$
$-x-2y+7z-2u=8$
$x-y-2z+6u=17$
を考える.
この方程式を上から順に対角線上の変数について解くと
$x=(20-2y-z-u)/9$
$y=(16-2x+2z-u)/8$
$z=(8+x+2y+2u)/7$
$u=(17-x+y+2z)/6$
となる.
$x,y,z,u$ に適当な値を入れて右辺を計算し、
得られた値を新たな $x,y,z,u$ として、計算を繰り返す.
漸化式で書くと
$x_{n+1}=(20-2y_n-z_n-u_n)/9$
$y_{n+1}=(16-2x_n+2z_n-u_n)/8$
$z_{n+1}=(8+x_n+2y_n+2u_n)/7$
$u_{n+1}=(17-x_n+y_n+2z_n)/6$

import std.stdio;
import std.math;

const int N = 4;

void main(string[] args)
{
    double[N][N] a = [[9,2,1,1],[2,8,-2,1],[-1,-2,7,-2],[1,-1,-2,6]]; 
    double[N]    b = [20,16,8,17];
    double[N]    c = [0,0,0,0];

    jacobi(a,b,c);
    writefln("X");
    for (int i = 0; i < N; i++)
        writef("%14.10f\t", c[i]);
    writefln("");
}
void jacobi(double[N][N] a, double[N] b, ref double[N] c)
{
    while (true)
    {
        double[N] s;
        bool finish = true;
        for (int i = 0; i < N; i++)
        {
            s[i] = b[i];
            for (int j = 0; j < N; j++)
                if (j != i)
                    s[i] -= a[i][j] * c[j];
            s[i] /= a[i][i];
            if (fabs(s[i] - c[i]) > 0.0000000001) finish = false;
        }
        for (int i = 0; i < N; i++)
        {
            c[i] = s[i];
            writef("%14.10f\t", c[i]);
        }
        writefln("");

        if (finish) return;
    }
}

  2.2222222222	  2.0000000000	  1.1428571429	  2.8333333333	
  1.3359788360	  1.3759920635	  2.8412698413	  3.1772486772	
  1.2477219283	  1.9791666667	  2.6346371882	  3.7870921517	
  1.0688819252	  1.8733422960	  2.9686056521	  3.8334531858	
  1.0501396189	  1.9957492835	  2.9260675555	  3.9569452792	
  1.0139431776	  1.9743638243	  2.9936469635	  3.9662907960	
  1.0101482880	  1.9991395970	  2.9850360597	  3.9912857623	
  1.0028221093	  1.9948112226	  2.9987141438	  3.9931772381	
  1.0020540192	  1.9998258539	  2.9969712901	  3.9982362335	
  1.0005711965	  1.9989497885	  2.9997397420	  3.9986190691	
  1.0004157346	  1.9999647527	  2.9993869874	  3.9996430127	
  1.0001156105	  1.9997874366	  2.9999473236	  3.9997204988	
  1.0000841449	  1.9999928659	  2.9998759259	  3.9999277456	
  1.0000233996	  1.9999569771	  2.9999893383	  3.9999434288	
  1.0000170310	  1.9999985561	  2.9999748873	  3.9999853757	
  1.0000047361	  1.9999912921	  2.9999978421	  3.9999885500	
  1.0000034471	  1.9999997077	  2.9999949172	  3.9999970400	
  1.0000009586	  1.9999982375	  2.9999995632	  3.9999976825	
  1.0000006977	  1.9999999408	  2.9999989712	  3.9999994009	
  1.0000001940	  1.9999996433	  2.9999999116	  3.9999995309	
  1.0000001412	  1.9999999880	  2.9999997918	  3.9999998787	
  1.0000000393	  1.9999999278	  2.9999999821	  3.9999999051	
  1.0000000286	  1.9999999976	  2.9999999579	  3.9999999755	
  1.0000000079	  1.9999999854	  2.9999999964	  3.9999999808	
  1.0000000058	  1.9999999995	  2.9999999915	  3.9999999950	
  1.0000000016	  1.9999999970	  2.9999999993	  3.9999999961	
  1.0000000012	  1.9999999999	  2.9999999983	  3.9999999990	
  1.0000000003	  1.9999999994	  2.9999999999	  3.9999999992	
  1.0000000002	  2.0000000000	  2.9999999997	  3.9999999998	
  1.0000000001	  1.9999999999	  3.0000000000	  3.9999999998	
  1.0000000000	  2.0000000000	  2.9999999999	  4.0000000000	
  1.0000000000	  2.0000000000	  3.0000000000	  4.0000000000	
X
  1.0000000000	  2.0000000000	  3.0000000000	  4.0000000000

ONLY DO WHAT ONLY YOU CAN DO

こけたら立ちなはれ立ったら歩きなはれ

D言語で連立一次方程式を解く(ヤコビの反復法)

参考文献