• 締切済み

 与関数get_objの或る点におけるヘッセ行列を求めるプログラムを下

 与関数get_objの或る点におけるヘッセ行列を求めるプログラムを下記の通り組んだのですが、点の座標を346879 2 2等、オーダーが離れた値を与えると、ヘッセ行列の要素がすべて0となるなど、理論とかけ離れた形として出力されることがあります。  与関数get_objはどの次元も2次関数であることから、ヘッセ行列は対角成分が2、それ以外は0となるはずなのに、このような出力が得られる場合があり、困っております。 原因は何か、教えて下さると助かります。 また、本プログラムを書いていて思ったのですが、多変数関数のn階の微分を行うプログラム作成法のヒントを教えて頂けないでしょうか? #include<iostream> #include<cmath> #include<ctime> using namespace std; const int N=3;//次元数 const double h=10e-4;//微分係数の定義用 int i,j,k,n;//for文用 double get_obj(double *x) {return x[0]*x[0]+x[1]*x[1]+x[2]*x[2];}//関数をセット void set_start(double *x); void dif(int n,double *x,double *grad);//n次成分を微分 void dif2(int i,int j,double *x,double **grad2);//第i,j成分を微分(2階微分) void set_grad(double *x,double *grad);//grad計算 void set_grad2(double *x,double grad2[][N]); void setH(double grad2[][N],double H[][N]);//H==grad2だから、別に使わなくても良い class parameter { public: double x[N]; double grad[N]; double grad2[N][N];//2階微分係数(ヘッセ行列に等しい) double H[N][N];//ヘッセ行列 }; int main() { parameter ag; set_start(ag.x); set_grad2(ag.x,ag.grad2); setH(ag.grad2,ag.H); return 0; } void set_start(double *x) { for(n=0;n<N;n++){ cout<<"第"<<n<<"次元の初期値を入力"<<endl; cin>>x[n]; } } void dif2(int i,int j,double *x,double grad2[N][N]) { double bunshi1,bunshi2,bunshi3,bunshi4; x[i]+=h;x[j]+=h; bunshi1=get_obj(x);//f(x_i+h,x_j+h)←x_i,x_j以外の引数は省略 x[i]-=h; bunshi2=get_obj(x);//f(x_j+h) x[i]+=h;x[j]-=h; bunshi3=get_obj(x);//f(x_i+h) x[i]-=h; bunshi4=get_obj(x);//f() grad2[i][j]=(bunshi1+bunshi4-bunshi2-bunshi3)/(h*h);//2階微分定義式 cout<<"grad2["<<i<<"]["<<j<<"]="<<grad2[i][j]<<endl; } void dif(int n,double *x,double *grad) { double bunshi1; double bunshi2; x[n]+=h; bunshi1=get_obj(x); x[n]-=h; bunshi2=get_obj(x); grad[n]=(bunshi1-bunshi2)/h; } void set_grad2(double *x,double grad2[N][N]) { for(i=0;i<N;i++){ for(j=0;j<N;j++){ dif2(i,j,x,grad2); } } } void setH(double grad2[][N],double H[][N]) { for(i=0;i<N;i++){ for(j=0;j<N;j++){ H[i][j]=grad2[i][j]; } }

みんなの回答

  • Tacosan
  • ベストアンサー率23% (3656/15482)
回答No.2

導関数ももらえれば簡単なんだけどねぇ....

  • Tacosan
  • ベストアンサー率23% (3656/15482)
回答No.1

「すべて 0 になる」というのは疑問ですが, 「桁落ち」してるような気がします. 対策は... この形では不可能かも.

expiz
質問者

お礼

回答ありがとうございます。 「この形」にこだわらければ、どのような解決法があるのでしょうか?

関連するQ&A