My Favorite Song: Hanamaru☆Sensation

My Favorite Song: Hanamaru☆Sensation

連立方程式の解法 反復法１ ヤコビ法

連立方程式の解法３　反復法１　ヤコビ法

次の３元連立方程式があるとする。

　　

①をx、②をy、③をzについて解くと次のようになる。

　　

ここで、として上の式の右辺に代入すると、

　　

となる。

このを用いて

　　

さらに、この値を用いて

このように繰り返し計算すると、

　　

と、やがて、３元連立方程式の解である(x,y,z)=(1,2,3)に収束してゆく。

このように反復計算をして連立方程式の解を求める方法をヤコビ法という。

より一般の

　　

つまり、

　　

の場合、次のように計算すればよい。

STEP1　１行からn行まで

　　

を計算する。

STEP2　i=1〜nまでのすべてのをの値に更新する

STEP3　所定の精度に達していなければSTEP1に戻る。達していたら、計算を打ち切る。

#include <stdio.h>
#include <math.h>
#define N 100


void Jacovi(double a[][N], double x[], double b[], int n) {
    double xn[n];
    double eps = 1.e-6,res, res_max;
    int i, j, k;
    int limit = 100;
   
    printf("反復回数 x(1)      x(2)      x(3)     最大残差\n");
    for (k = 1 ; k <= limit; k++) {
        res_max = 0;
        for (i=0; i < n; i++) {
            res = b[i];
            for (j=0; j < n; j++) {
                res = res - a[i][j]*x[j];
            }
            if (res_max < fabs(res)) res_max=fabs(res);
            xn[i]= x[i]+res/a[i][i];
        }
        if (k%5==0) { // ５回に１度、計算途中の結果を表示
            printf("%2d %10.6f %10.6f %10.6f %10.6f\n", k, xn[0],xn[1], xn[2], res_max);
        }
        for (i=0; i<n; i++)
            x[i]=xn[i]; //xの更新
        if (eps>res_max) break;
    }
}

main() {
    int i, n;
    double a[N][N] = {{3,1,1}, {1,3,1}, {1,1,3}};
    double x[N]={0.};
    double b[N]= {8., 10., 12.};
   
    n=3;
   
    Jacovi(a,x,b,n);

    for (i = 0; i < n; i++) {
        printf("x[%d] = %f\n", i+1,x[i]);
    }
}

いま述べたヤコビ法は、直接、連立方程式の解を求める掃き出し法やガウス消去法とは異なり、連立方程式が唯一の解を持っていたとしても、計算が収束しなかったり、発散する場合がある。

優対角条件、つまり、

などの条件が成立しないとき、ヤコビ法などの反復法は振動したり、発散する。

このことは

　　

の②と③を入れ換え

　　

をヤコビ法で解いてみると、このことが確かめられる。

My Favorite Song: Ochame Kinou

My Favorite Song: Ochame Kinou

Orchestra Ver.

平行平板間の流れ（平面ポアズイユ流れ）

平行平板間の流れ（平面ポアズイユ流れ）

右図のような２つの無限平行板間の流れは、次の方程式で表される。

ここで、pは圧力、uはx軸に平行な流速、さらに、Qは（体積）流量、μは粘性係数である。

なお、（１）の境界条件は

（１）は、流体力学の基礎方程式であるナビエ・ストークス方程式

の最も簡単なものの一つで、ナビエ・ストークスの解析解（厳密解）が得られるものの一つで、平均流速

を用いると、解は

である。

u、yを

と無次元化すると、（５）式は

となり、UはY=1/2のときに

という最大値をとる。

さらに、pとxを次のように無次元化すると

（１）、（２）式は

ここで、Reはレイノルズ数と言われる、流体力学で流れを支配する、重要な無次元数で

である。

実は、このReが大きくなると、ナビエストークス方程式の非線形項の非線形性が強くなり、流れは時間、空間的に不規則――不規則って何だ(・・?――な乱流という、数学・物理学的に厄介なものになるのだけれど・・・。

それで、今回、やりたいのは、差分法を用いて、（１０）と（１１）を連立方程式に直し、これを数値的に解こうということ。

（７）式を見るとわかるけれど、この流れの速度分布は、Y=1/2を軸とする放物線で表されており、Y=1/2で対称である。

したがって、

だから、

（１０）式の境界条件は、

ではなく、

としてもよい。

なぜ、境界条件を、（１４）ではなく、（１５）とするかというと、計算領域を0≦Y≦1から0≦Y≦1/2と小さくし、計算量をできるだけ少なくしたいから。

[0,1/2]をn等分し、

とおき、

とし、この点における無次元流速Uをとする。

差分法を用いると

となるので、（１０）式は

となり、

ここで、

とおけば、（１９）式は

になる。

Y=0の境界条件U=0より

問題は、Y=1/2における境界条件

をどのように扱うかだけれど、中心差分を用いると、

となるので、i=nのとき

とすればよいだろう。

（１１）式も対称性を利用すると、

となり、左辺の積分を台形公式で近似すれば、

となる。

よって、

というn+1元の連立方程式が得られる。

未知数は、のn個とdP/dXの１個で、合せてn+1個だから、連立方程式（２５）を解くことによって、この未知数を全て求めることができる。

（２５）は行列で書き換えると、

したがって、係数行列を

にすることにより、ガウス消去法や掃き出し法を用いて連立方程式（２６）を解くことができる！！

（２６）の左辺の行列の成分の圧倒的多数は０の疎行列なので、単純なガウス消去法などの直接法は無駄な計算が多く、効率的ではないのだけれど・・・。

n=10の場合の計算結果は次のとおり。

速度U、圧力勾配dP/dXの計算誤差は約0.25％に収まっていることがわかる。

実は、ここまでのプログラムを作ると、それをもとに、より複雑なNS方程式の特殊形

で表される無限平板間の流れを数知的に解くこともできる！！

なのですが、上の方程式は非線形な上に、計算点が少なくとも数百〜数千程度になるので、ガウス消去法などの直接法で、上の連立微分方程式から得られる連立方程式を解くなど狂気の沙汰！！

連立方程式

を効率よく解く数値計算法、つまり、アルゴリズムを作らないといけない。

ガウス消去法の計算量はn³程度になるという話を前にしたけれど、ガウス消去法をもとに、計算量をnオーダーに留める計算手法、アルゴリズムを作ることは可能！！

我はと思うヒトは、この計算手法を考えてみるといいにゃ。

それさえできれば、

を数値的に解くプログラムを作ることもできる！！