以前にも同じ質問をさせていただいたのですが、解決できなかったので再送します。
Gnome端末を使っています。
x1+2x2+3x3=6
3x1+10x2-4x3=-29
-2x1-4x2+x3=9
という連立方程式をGauss-seidel法で求める問題です。
下記のようにプログラミングすると、
ILL CONDITION i=0
Iteration x1 x2 x3
1 1.0000 1.1034 1.7126
2 -0.2241 0.7406 1.2793
3 0.1135 0.8353 1.3964
4 0.0233 0.8098 1.3651
5 0.0475 0.8166 1.3735
6 0.0410 0.8148 1.3713
7 0.0428 0.8153 1.3719
8 0.0423 0.8152 1.3717
9 0.0424 0.8152 1.3717
10 0.0424 0.8152 1.3717
11 0.0424 0.8152 1.3717
となり、x1=1、x=ー2、x3=3が得られません。
また、
3x1+x2+x3=10
x1+5x2+2x3=21
x1+2x2+5x3=30
を同じプログラミングで解こうとすると、
Iteration x1 x2 x3
1 1.0000 0.9524 0.9032
2 0.8144 0.8752 0.9145
3 0.8210 0.8738 0.9144
4 0.8212 0.8738 0.9144
5 0.8212 0.8738 0.9144
のようになり、x1=1、x=2、x3=5が得られません。
“♯include<stdio.h>
♯include<math.h>
♯define eps 1.e-5
int main(void)
{
double a[3][3],b[3],x[3];
double s,w,anorm,xnorm;
int i,j,k,n;
scanf("%d",&n);
for(i=0;i<n;i++){
for(j=0;j<n;j++){
scanf("%lf",&a[i][j]);
}
scanf("%lf",&b[i]);}
for(i=0;i<n;i++){
s=0.0;
for(j=0;j<n;j++){
if(i==j) continue;
s+=fabs(a[i][j]);
}
if(fabs(a[i][i])<=s){
printf("ILL CONDITION i=%d\n",i);
break;
}
}
for(i=0;i<n;i++){
x[i]=0.0;
}
printf("Iteration x1 x2 x3\n");
for(k=1;;k++){
anorm=0.0;
xnorm=0.0;
for(i=0;i<n;i++){
w=b[i];
for(j=0;j<n;j++){
if(j==i) continue;
w -=a[i][j]*x[j];
}
w /=b[i];
anorm+=fabs(x[i]-w);
xnorm+=fabs(w);
x[i]=w;
}
printf("%6d",k);
for(i=0;i<n;i++){
printf(" %.4f",x[i]);
}
printf("\n");
if(anorm/xnorm<eps)
break;
}
return 0;
}
”
御指摘いただきたいです。
-
気になる質問をクリップする
クリップした質問は、後からいつでもマイページで確認できます。
またクリップした質問に回答があった際、通知やメールを受け取ることができます。
クリップを取り消します
-
良い質問の評価を上げる
以下のような質問は評価を上げましょう
- 質問内容が明確
- 自分も答えを知りたい
- 質問者以外のユーザにも役立つ
評価が高い質問は、TOPページの「注目」タブのフィードに表示されやすくなります。
質問の評価を上げたことを取り消します
-
評価を下げられる数の上限に達しました
評価を下げることができません
- 1日5回まで評価を下げられます
- 1日に1ユーザに対して2回まで評価を下げられます
質問の評価を下げる
teratailでは下記のような質問を「具体的に困っていることがない質問」、「サイトポリシーに違反する質問」と定義し、推奨していません。
- プログラミングに関係のない質問
- やってほしいことだけを記載した丸投げの質問
- 問題・課題が含まれていない質問
- 意図的に内容が抹消された質問
- 過去に投稿した質問と同じ内容の質問
- 広告と受け取られるような投稿
評価が下がると、TOPページの「アクティブ」「注目」タブのフィードに表示されにくくなります。
質問の評価を下げたことを取り消します
この機能は開放されていません
評価を下げる条件を満たしてません
質問の評価を下げる機能の利用条件
この機能を利用するためには、以下の事項を行う必要があります。
- 質問回答など一定の行動
-
メールアドレスの認証
メールアドレスの認証
-
質問評価に関するヘルプページの閲覧
質問評価に関するヘルプページの閲覧
checkベストアンサー
+1
なんとなく質問者さんが理解されていないように感じるため、修正したコードの全文を掲載します。
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#define eps 1.e-5
int main(void)
{
double a[3][3],b[3],x[3];
double s,w,anorm,xnorm;
int i,j,k,n;
scanf("%d", &n);
for(i=0; i<n; i++){
for(j=0; j<n; j++){
scanf("%lf", &a[i][j]);
}
scanf("%lf", &b[i]);
}
for(i=0;i<n;i++){
s=0.0;
for(j=0;j<n;j++){
if(i==j) continue;
s+=fabs(a[i][j]);
}
if(fabs(a[i][i])<=s){
printf("ILL CONDITION i=%d\n",i);
exit(1);
}
}
for (i = 0; i < n; i++) {
for (j = 0; j < n; j++) {
printf("%lfx%d ", a[i][j], j+1);
}
printf("= %lf\n", b[i]);
}
for(i=0; i < n; i++){
x[i] = 0.0;
}
printf("Iteration x1 x2 x3\n");
for(k=1; ; k++){
anorm = 0.0;
xnorm = 0.0;
for(i = 0; i < n; i++){
w=b[i];
for(j=0; j<n; j++){
if(j==i) continue;
w -= a[i][j]*x[j];
}
// w /= b[i]; // ←ここが間違っている
w /= a[i][i];
anorm += fabs(x[i]-w);
xnorm += fabs(w);
x[i] = w;
}
printf("%6d", k);
for(i = 0; i < n; i++){
printf(" %.4f", x[i]);
}
printf("\n");
if(anorm/xnorm < eps) break;
}
return 0;
}
ガウス=ザイデル法の更新式をよく見なおしてください。
なお、ひとつ目の例は、
ILL CONDITION i=0
と出力されているように、ガウス=ザイデル法の収束条件(対角優位)を満たしていないため収束しません。
以上のように修正すれば、2つ目の例は収束します。
% ./a.out
3
3 1 1 10
1 5 2 21
1 2 5 30
Iteration x1 x2 x3
1 3.3333 3.5333 3.9200
2 0.8489 2.4622 4.8453
3 0.8975 2.0824 4.9876
4 0.9767 2.0096 5.0008
5 0.9965 2.0004 5.0005
6 0.9997 1.9998 5.0001
7 1.0000 1.9999 5.0000
8 1.0000 2.0000 5.0000
投稿
-
回答の評価を上げる
以下のような回答は評価を上げましょう
- 正しい回答
- わかりやすい回答
- ためになる回答
評価が高い回答ほどページの上位に表示されます。
-
回答の評価を下げる
下記のような回答は推奨されていません。
- 間違っている回答
- 質問の回答になっていない投稿
- スパムや攻撃的な表現を用いた投稿
評価を下げる際はその理由を明確に伝え、適切な回答に修正してもらいましょう。
0
インデントを整えたり、コメントを入れて、プログラムをわかりやすくしましょう。
バグが発見しやすくなります。
2番目の連立方程式は、解x1=1、x=2、x3=5が得られました。
1番目の連立方程式は、解がGauss-Seidel 法では発散するようです。
#include<stdio.h>
#include<math.h>
/* 計算誤差の許容値 */
#define eps (1.e-5)
int main(void)
{
double a[3][3], b[3], x[3];
double s, anorm, xnorm, sum;
int i, j, k;
int n; /* 連立方程式の大きさ */
double new_x; /* 新しい近似解 */
/* データ入力 */
scanf("%d", &n);
for (i = 0; i < n; i++) {
for (j = 0; j < n; j++) {
scanf("%lf", &a[i][j]);
}
scanf("%lf", &b[i]);
}
for (i = 0; i < n; i++) {
s = 0.0;
for (j = 0; j < n; j++) {
if (i == j)
continue;
s += fabs(a[i][j]);
}
if (fabs(a[i][i]) <= s) {
printf("ILL CONDITION i=%d\n", i);
/* この条件のとき、近似解は発散? */
break;
}
}
/* 近似解の初期値代入 */
for (i = 0; i < n; i++) {
x[i] = 0.0;
}
printf("Iteration x1 x2 x3\n");
for (k = 1; ; k++) {
anorm = 0.0;
xnorm = 0.0;
for (i = 0; i < n; i++) {
sum = 0.0;
for (j = 0; j < n; j++) {
if (j == i)
continue;
sum += a[i][j] * x[j];
}
/* 反復計算後の近似解 */
new_x = 1.0 / a[i][i] * (b[i] - sum);
/* 近似解の変化量を加算 */
anorm += fabs(x[i] - new_x);
/* 近似解の総和計算 */
xnorm += fabs(new_x);
/* 新しい近似解を代入 */
x[i] = new_x;
}
printf("%6d", k);
for (i = 0; i < n; i++) {
printf(" %.4f", x[i]);
}
printf("\n");
/* 計算終了条件 */
if (anorm / xnorm < eps)
break;
}
return 0;
}
投稿
-
回答の評価を上げる
以下のような回答は評価を上げましょう
- 正しい回答
- わかりやすい回答
- ためになる回答
評価が高い回答ほどページの上位に表示されます。
-
回答の評価を下げる
下記のような回答は推奨されていません。
- 間違っている回答
- 質問の回答になっていない投稿
- スパムや攻撃的な表現を用いた投稿
評価を下げる際はその理由を明確に伝え、適切な回答に修正してもらいましょう。
0
ちょこっとwikiでGauss-seidel法をみただけですがGauss-seidel法のアルゴリズムとコードを見比べておかしいと思える点をコメントします。
前回の質問で正しいコードを間違ったものに直してます。
w /= a[i][i]
が正しいですが、これをw /= b[i]
に修正してしまってます。コメントした方は「Gauss-seidel法については知らないがコードをみた感じ間違いではないか?」とおっしゃっているだけですので、その指摘が正しいかどうかは質問者さんご自身で判断すべきだったと思います。収束条件
Gauss-seidelで解が収束する充分条件は(必要条件かはわかりません)「係数行列の各行で非対角要素の絶対値の和が対角要素の絶対値よりも小さい場合」とwikiにあります。プログラムコードでもそれをチェックしてはいますが単にメッセージを表示するだけで実行は継続していますね?最初の方程式は1行目がこの条件を満たしていませんので"IL CONDITION"とメッセージが出ます。実行を継続してみるとおっしゃるとおり収束しないように見えます。エラーメッセージを出すだけでプログラムを中断しない論理となっているのはどういう意図でしょう?「充分条件は満たしていなくてもとにかくやってみる」ということならわかりますが・・・当然ながら本当に収束できないケースでは無限ループに陥ってしまいます。収束条件を満たさないケースで処理をあきらめたいなら「収束条件を満たさないことが判明した時点で」returnなりexitで実行を中断すべきでしょう。
なお、収束条件を満たすように式を変形するような配慮をプログラム上でしてやれば解が求まるとは思いますが、そのような論理を入れると最早Gauss-seidel法ではなくGaussの消去法になってしまうのかも知れません。
投稿
-
回答の評価を上げる
以下のような回答は評価を上げましょう
- 正しい回答
- わかりやすい回答
- ためになる回答
評価が高い回答ほどページの上位に表示されます。
-
回答の評価を下げる
下記のような回答は推奨されていません。
- 間違っている回答
- 質問の回答になっていない投稿
- スパムや攻撃的な表現を用いた投稿
評価を下げる際はその理由を明確に伝え、適切な回答に修正してもらいましょう。
15分調べてもわからないことは、teratailで質問しよう!
- ただいまの回答率 88.12%
- 質問をまとめることで、思考を整理して素早く解決
- テンプレート機能で、簡単に質問をまとめられる
質問への追記・修正、ベストアンサー選択の依頼
KSwordOfHaste
2017/01/17 16:13
コードを載せるときは```で囲み字下げしましょう。