連立方程式の解法3 反復法1 [ネコ騙し数学]
連立方程式の解法3 反復法1
ヤコビ法
次の3元連立方程式があるとする。
①をx、②をy、③をzについて解くと次のようになる。
ここで、として上の式の右辺に代入すると、
となる。
このを用いて
さらに、この値を用いて
このように繰り返し計算すると、
と、やがて、3元連立方程式の解である(x,y,z)=(1,2,3)に収束してゆく。
このように反復計算をして連立方程式の解を求める方法をヤコビ法という。
より一般の
つまり、
の場合、次のように計算すればよい。
STEP1 1行からn行まで
を計算する。
STEP2 i=1〜nまでのすべてのをの値に更新する
STEP3 所定の精度に達していなければSTEP1に戻る。達していたら、計算を打ち切る。
いま述べたヤコビ法は、直接、連立方程式の解を求める掃き出し法やガウス消去法とは異なり、連立方程式が唯一の解を持っていたとしても、計算が収束しなかったり、発散する場合がある。
優対角条件、つまり、
などの条件が成立しないとき、ヤコビ法などの反復法は振動したり、発散する。
このことは
の②と③を入れ換え
をヤコビ法で解かせると、このことが確かめられる。
C言語によるサンプルプログラムを次に示す。
#include <stdio.h>
#include <math.h>
#define N 100
void Jacovi(double a[][N], double x[], double b[], int n) {
double eps = 1.e-6, sum=0, res=0;
double xn[N];
int i, j, k;
int limit = 100;
printf("反復回数 x(1) x(2) x(3) 残差\n");
for (k = 1 ; k <= limit; k++) {
res = 0;
for (i=0; i < n; i++) {
sum = b[i];
for (j=0; j < n; j++) {
sum = sum - a[i][j]*x[j];
}
res = fabs(sum);
xn[i]= x[i]+sum/a[i][i];
}
if (k%5==0) { // 5回に1度、計算途中の結果を表示
printf("%2d %10.6f %10.6f %10.6f %10.6f\n", k, xn[0],xn[1], xn[2], res);
}
for (i=0; i<n; i++) {
x[i]=xn[i];
}
if (eps>res) 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]);
}
}
#include <math.h>
#define N 100
void Jacovi(double a[][N], double x[], double b[], int n) {
double eps = 1.e-6, sum=0, res=0;
double xn[N];
int i, j, k;
int limit = 100;
printf("反復回数 x(1) x(2) x(3) 残差\n");
for (k = 1 ; k <= limit; k++) {
res = 0;
for (i=0; i < n; i++) {
sum = b[i];
for (j=0; j < n; j++) {
sum = sum - a[i][j]*x[j];
}
res = fabs(sum);
xn[i]= x[i]+sum/a[i][i];
}
if (k%5==0) { // 5回に1度、計算途中の結果を表示
printf("%2d %10.6f %10.6f %10.6f %10.6f\n", k, xn[0],xn[1], xn[2], res);
}
for (i=0; i<n; i++) {
x[i]=xn[i];
}
if (eps>res) 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]);
}
}
タグ:数値解析
2017-10-29 13:03
nice!(0)
コメント(0)
コメント 0