- Jacobi.JPG (47.68 KiB) 閲覧数: 9150 回
先ほどのJacobi法の計算で∞-ノルムの場所に関数を使用してプログラムを書きたいです。途中まで書いてみたのですがここから分からないです。教えてほしいです。
コード:
#include <stdio.h>
#include <math.h>
#define n 5 /* 行列サイズ */
#define Nmax 100 /* 最大反復回数 */
#define eps 1.0e-10 /* 許容誤差 */
/* 関数のプロトタイプ宣言 */
void matVec(double *Ax, double (*A)[n], double *x);
double norm2(double *x);
double normInf(double *x);
int main(){
int i, j, k;
int iter;
double sum;
/* 係数行列 */
double A[n][n] = {{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[n] = {9, 11, 11, 11, 11};
double x_old[n]; /* 更新前の近似解 */
double x_new[n]; /* 更新後の近似解 */
double err[n]; /* 誤差 */
double r[n]; /* 残差 */
/* 初期近似解の設定 */
for(i=0; i<n; i++){
x_old[i] = 0.0;
}
/* Jacobi法の反復開始 */
iter=0;
while(1){
/* 反復回数をカウント */
iter++;
/* 近似解の更新 */
for(i=0; i<n; i++){
/* 分からない */
}
/* 誤差(隣り合う近似解の差分)を計算 */
for(i=0; i<n; i++){
err[i] = x_new[i] - x_old[i];
}
/* 収束判定 */
if(normInf(err) < eps * normInf(x_new)){
break;
}
/* 最大反復回数の確認 */
if(iter == Nmax){
printf("収束しない(最大反復回数に到達しました)");
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]);
}
/* 残差2-ノルムの表示 */
matVec(r,A,x_new); /* 行列ベクトル積 A*x の結果を一時的に r に保存 */
for(i=0; i<n; i++){
r[i] = b[i] - r[i]; /* r=b-A*x を計算 */
}
printf("残差2ノルム ||b-Ax||_2 = %.2e\n", norm2(r));
return 0;
}
/* 行列ベクトル積を計算する関数 */
void matVec(double *Ax, double (*A)[n], double *x){
/* 分からない */
}
/* ベクトルの2-ノルムを計算する関数 */
double norm2(double *x){
/* 分からない */
}
/* ベクトルのinf-ノルムを計算する関数 */
double normInf(double *x){
/* 分からない */
}