与関数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成分を微分（２階微分）
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];//２階微分係数（ヘッセ行列に等しい）
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);//２階微分定義式
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];
}
}
お礼
回答ありがとうございます。 もしわかりましたら回答お願いします。