間違い？jacobiを用いたnewton法

質問

以下のプログラムを作ったんですが実行したときに答えがnanになります。どなたか間違いがわかる方よろしくおねがします。

const double e =1e-10;
enum {N=100};

const double x0 =1;
const double y0 =1;

double f(const double x,double y) {return x*x-y*y;}
double g(const double x,double y) {return 6*x*y-5;}
double dfdx(const double x,double y){return 2*x;}
double dfdy(const double x,double y){return -2*y;}
double dgdx(const double x,double y){return 6*y;}
double dgdy(const double x,double y){return 6*x;}
double det(double x,double y);

main()
{
double x[N+1],y[N+1];
int i;
double det2;

x[0]=x0;
y[0]=y0;

for(i=0;i<=20;i++)
{

det2=det(x[i],y[i]);

x[i+1]=x[i]-(dgdy(x[i],y[i])*f(x[i],y[i]))-(dfdy(x[i],y[i])*g(x[i],y[i]))/det2;

y[i+1]=y[i]-(dgdx(x[i],y[i])*f(x[i],y[i]))-(dfdx(x[i],y[i])*g(x[i],y[i]))/det2;

printf("x = %lf\n",x[i+1]);
printf("y = %lf\n",y[i+1]);
printf("\n");

}

exit(0);
}

double det(double a,double b)
{
double ans;
ans=1/(dfdx(a,b)*dgdy(a,b)-dfdy(a,b)*dgdx(a,b));

return ans;
}

回答No.3

Interest · ベストアンサー

ANo.1 = 2 です。
バラバラ書いてごめんなさい。

> x[i+1]=x[i]-(dgdy(x[i],y[i])*f(x[i],y[i]))-(dfdy(x[i],y[i])*g(x[i],y[i]))/det2;

これはいいけど、

> y[i+1]=y[i]-(dgdx(x[i],y[i])*f(x[i],y[i]))-(dfdx(x[i],y[i])*g(x[i],y[i]))/det2;

これ、プラスマイナスの付け方間違っていませんか？

というのは、ニュートン法を使うためにヤコビ行列の逆行列を使いますよね？　まずは、逆行列が存在しないケース（すなわち行列式＝０の場合）が抜けていませんか、というのがANo.1,2 の指摘。今回の指摘は、逆行列に f(x,y), g(x,y) を書けるのにプラスマイナスが間違っちゃってませんか、という指摘です。

プログラム云々の前に手で数式を展開して、最初の１、２回目程度手計算した結果とプログラムが吐き出す答えをつき合わしてみてはいかがでしょうか。

回答No.2

Interest

というか、その前に、

> double det(double a,double b)

という関数はきっと行列式を計算したいのだと思うけど、

> ans=1/(dfdx(a,b)*dgdy(a,b)-dfdy(a,b)*dgdx(a,b));

って、
det = dfdx(a,b)*dgdy(a,b)-dfdy(a,b)*dgdx(a,b);
じゃないですか？

回答No.1

Interest

> 実行したときに答えがnanになります。

というのは、ループの何回目からですか？
そのときのx,yの値は？

ans=1/(dfdx(a,b)*dgdy(a,b)-dfdy(a,b)*dgdx(a,b));
で、分母が０になる可能性は？
まずはここにエラー処理を入れてみませんか？

それから、

> const double e =1e-10;

というのはどこかで使ってますか？

> enum {N=100};

としながら、forループが２０回なのはなぜですか？