ONLY DO WHAT ONLY YOU CAN DO

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

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

例として
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

#include <iostream>
#include <iomanip>
#include <math.h>
using namespace std;

void jacobi(double a[4][4], double b[4], double c[4]);

int main()
{
    double a[4][4] = {{9,2,1,1},{2,8,-2,1},{-1,-2,7,-2},{1,-1,-2,6}}; 
    double b[4]    = {20,16,8,17};
    double c[4]    = {0,0,0,0};

    jacobi(a,b,c);
    for (int i = 0;i<=3;i++)
        cout << setw(12) << fixed << setprecision(10) << c[i] << "\t";
    cout << endl;
   
    return 0;
}

void jacobi(double a[4][4], double b[4], double c[4])
{
    for (int n = 0;n<=100;n++)
    {
        double s[4];
        double e = 0;
        for (int i = 0;i<=3;i++)
        {
            s[i] = b[i];
            for (int j=0;  j<=(i-1);j++) s[i] -= a[i][j] * c[j];
            for (int j=i+1;j<=3;    j++) s[i] -= a[i][j] * c[j];
            s[i] /= a[i][i];
            e += fabs(s[i] - c[i]);
        }
        for (int i=0;i<=3;i++)
        {
            c[i] = s[i];
            cout << setw(12) << fixed << setprecision(10) << c[i] << "\t";
        }
        cout << endl;

        if (e <= 0.0000000001) return;
    }
    cout << "収束しない" << endl;
}
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	
1.0000000000	2.0000000000	3.0000000000	4.0000000000	
1.0000000000	2.0000000000	3.0000000000	4.0000000000	
参考文献