jacobiのNewton法

質問

何回やっても収束しません。
初期値は共に1.0でε（誤差）は0.001ぐらいで回数を1000にしてます。
答えがありえないぐらい大きくなったりします。
どこがちがうでしょうか？わかる方お願いします。

#include<stdio.h>
#include<math.h>

double xn[10][100],ipu,jacobi,gyou[2],gyouretu[2][2],a[2];
int genkai,i,j;
void nyuryoku(void);
double f(double,double);
double g(double,double);
double fx(double,double);
double fy(double,double);
double gx(double,double);
double gy(double,double);
void keisan(void);
void syuturyoku(void);

int main(void)
{
nyuryoku();
keisan();
syuturyoku();
exit(0);
}

void nyuryoku(void)
{
for(j=0;j<2;j++){
printf("初期値x,yを入力してください");
scanf("%lf",&xn[j][0]);
}
printf("εを入力してください。");
scanf("%lf",&ipu);
printf("収束しないと判断する回数は？");
scanf("%d",&genkai);
}

void keisan(void)
{
for(i=0;i<genkai;i++){
gyouretu[0][0]=fx(xn[0][i],xn[1][i]);
gyouretu[0][1]=fy(xn[0][i],xn[1][i]);
gyouretu[1][0]=gx(xn[0][i],xn[1][i]);
gyouretu[1][1]=gy(xn[0][i],xn[1][i]);
a[0]=f(xn[0][i],xn[1][i]);
a[1]=g(xn[0][i],xn[1][i]);
gyou[0]=gyouretu[0][0]*a[0]+gyouretu[0][1]*a[1];
gyou[1]=gyouretu[1][0]*a[0]+gyouretu[1][1]*a[1];
jacobi=1/(gyouretu[1][1]*gyouretu[0][0]-gyouretu[0][1]*gyouretu[1][0]);
xn[0][i+1]=xn[0][i]-jacobi*gyou[0];
xn[1][i+1]=xn[1][i]-jacobi*gyou[1];

if(fabs(f(xn[0][i+1],xn[1][i+1]))<ipu && fabs(g(xn[0][i+1],xn[1][i+1]))<ipu)break;
}
}

double f(double x,double y)
{
double z;
z=x*x-y*y;
return z;
}

double g(double x,double y)
{
double z;
z=6*x*y-5;
return z;
}

double gy(double x,double y)
{
double z;
z=6*y;
return z;
}

double fy(double x,double y)
{
double z;
z=-2*y;
return z;
}
double gx(double x,double y)
{
double z;
z=6*x;
return z;
}

double fx(double x,double y)
{
double z;
z=2*x;
return z;
}

void syuturyoku(void)
{
for(j=0;j<=i;j++){
printf("X%d=%.15lf\n",j,xn[0][j]);
printf("y%d=%.15lf\n",j,xn[1][j]);
}
printf("ε=%.15lf\n",ipu);
}

回答No.1

asuncion · ベストアンサー

理論上、どういった値に収束するはずなのですか？

また、最大計算回数を1000にしているとのことですが、
配列定義が double xn[10][100] であるのに1000回も計算したら、
配列の定義範囲を超えてしまいます。
これはまずいです。

回答No.4

Interest

その後進展はありましたか？

Javaのカテゴリでも
http://okwave.jp/qa3034520.html
で同じ質問してたんですね。

逆行列の求め方ができていないのが原因ですから、それを直さない限りどこで質問しても一緒です。

回答No.3

Interest

こういうのみると、かっち～んときます(-_-#)

http://oshiete1.goo.ne.jp/qa3038640.html
（http://okwave.jp/qa3038640.html）
のほうが、まだプログラムの書き方がマシでしたよ。

今回の式でも前回とまったく同じことが考慮できていなくて、
ヤコビ行列の逆行列を求める計算が間違っています。
１．行列式＝０になったら逆行列が存在しないので、計算できない（ので回避策を打つ必要がある）こと
２．逆行列を求める際に、プラスマイナスを間違えていること
この２つが解決できなければ、何度プログラムを書いても間違えたままです。

ヤコビ行列、なのに
jacobi=1/(gyouretu[1][1]*gyouretu[0][0]-gyouretu[0][1]*gyouretu[1][0]);
って、これは 1/det ですよね。det = 0 ならどうする？

課題丸投げせずに、あきらめず何度も自力で挑戦しているところだけは、感心しておきます。

回答No.2

tatsu99

http://oshiete1.goo.ne.jp/qa3038640.html
の質問はどうなったのですか？