数値解析 Jacobi法

[うーちゃん]

Jacobi.JPG (47.68 KiB) 閲覧数: 2772 回

数値解析のプログラムが分からなくて困ってます。ある程度は完成させたのですがここから先進まなくて困ってます。分かる方お願いします。
[img]





#include <stdio.h>
#include <math.h>
int main(){
int i, j, k, n;
int iter, Nmax;
double eps, sum, norm_err, norm_x;

/* 行列サイズ */
n = 5;

/* 係数行列 */
double A[5][5] = {{6, 1, 1, 1, 0},
{1, 7, 1, 1, 1},
{0, 1, 8, 1, 1},
{0, 0, 1, 9, 1},
{0, 0, 0, 1, 10}};

/* 右辺項 */
double b[5] = {9, 11, 11, 11, 11};

double x_old[5]; /* 更新前の近似解 */
double x_new[5]; /* 更新後の近似解 */
double err[5]; /* 誤差 */

/* 最大反復回数 */
Nmax = 100;

/* 許容誤差 */
eps = 1.0e-10;

/* 初期近似解の設定 */
for(i=0; i<n; i++){
x_old = 0.0;
}

/* Jacobi法の反復開始 */
iter=0;
while(1){

/* 反復回数をカウント */
iter++;

/* 近似解の更新 */
for(i=0; i<n; i++){



/* 分からない*/
}

/* 誤差(隣り合う近似解の差分）を計算 */
for(i=0; i<n; i++){
err = x_new - x_old;
}

/* 誤差の∞ノルムを計算 */
norm_err = 0.0;
for(i=0; i<n; i++){
if(fabs(err) > norm_err){
norm_err = fabs(err);
}
}

/* 近似解の∞ノルムを計算 */
norm_x = 0.0;
for(i=0; i<n; i++){

/* 分からない*/

}

/* 収束判定 */
if(norm_err < eps * norm_x){
break;
}

/* 最大反復回数の確認 */
if(iter == Nmax){
printf("収束しない（最大反復回数に到達しました）");
break;
}

/* 近似解の入れ替え */
for(i=0; i<n; i++){
x_old = x_new;
}
}

/* 反復回数の表示 */
printf("反復回数：%d\n\n",iter);

/* 近似解の表示 */
printf("近似解 x=\n");
for(i=0; i<n; i++){
printf("%.16g\n", x_new);
}

return 0;
}

[かずま]

コード:

codeタグを使わないと、[i] がイタリック体タグと解釈されて削除され、
それ以後、半角文字がイタリック体になるので注意してください。

これでよいのでしょうか？

コード:

#include <stdio.h>  // printf
#include <math.h>   // fabs

int main()
{
	int i, j;
	double norm_err, norm_x;
	const int n = 5;         // 行列サイズ
	double A[5][5] = {       // 係数行列
		{ 6, 1, 1, 1, 0 },
		{ 1, 7, 1, 1, 1 },
		{ 0, 1, 8, 1, 1 },
		{ 0, 0, 1, 9, 1 },
		{ 0, 0, 0, 1, 10},
	};
	double b[5] = { 9, 11, 11, 11, 11 };   // 右辺項
	double x_old[5] = { 0.0 };     // 更新前の近似解
	double x_new[5];               // 更新後の近似解
	double err[5];                 // 誤差
	const int Nmax = 100;          // 最大反復回数
	const double eps = 1.0e-10;    // 許容誤差
	int iter = 0;                  // Jacobi法の反復開始

	while (1) {
		iter++;                    // 反復回数をカウント
		for (i = 0; i < n; i++) {  // 近似解の更新
			x_new[i] = b[i];
			for (j = 0; j < n; j++)
				if (j != i) x_new[i] -= A[i][j] * x_old[i];
			x_new[i] /= A[i][i];
		}
#if 0  // 0 を 1 に変えると、途中経過の表示
		printf("%3d:", iter);
		for (i = 0; i < n; i++) printf("%14.10f", x_new[i]);
		putchar('\n');
#endif
		// 誤差(隣り合う近似解の差分）を計算
		for (i = 0; i < n; i++) err[i] = x_new[i] - x_old[i];
		norm_err = 0.0;  // 誤差の∞ノルムを計算
		for (i = 0; i < n; i++)
			if (fabs(err[i]) > norm_err) norm_err = fabs(err[i]);
		norm_x = 0.0;   // 近似解の∞ノルムを計算
		for (i = 0; i < n; i++)
			if (fabs(x_old[i]) > norm_x) norm_x = fabs(x_old[i]);
		if (norm_err < eps * norm_x) break;    // 収束判定
		if (iter == Nmax) {         // 最大反復回数の確認
			puts("収束しない（最大反復回数に到達しました）");
			break;
		}
		for (i = 0; i < n; i++) x_old[i] = x_new[i];  // 近似解の更新
	}
	printf("反復回数：%d\n\n",iter);
	printf("近似解 x=\n");
	for (i = 0; i < n; i++) printf("%.16g\n", x_new[i]);
	return 0;
}