newton法のjacobi

質問

以下のプログラムを作ったんですが実行したときに答えが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.1

sakusaker7 · ベストアンサー

これはJavaのプログラムじゃなくてCでは?

とりあえずCコンパイラで実行ファイルを作って試してみましたが、
xの値が発散する一方のようです。

detの式がどういうものかは知りませんが、収束するようなもので
ないとニュートン法では解は求まらないのではないでしょうか?

初期値がよくないものでないか、f'(x) が本当に正しいのか
再確認してみてはいかがでしょう。