erjiedaoshu.cpp

用C++编的多元函数求二阶导数

#include "iostream.h"
#include "math.h"
void hessian(double (*pf)(double *x),int n,double *point,double *H);
double fun(double *x);

void main(void)
{
    int n=2;
    double x[2]={1,1};
double H[4];
hessian(fun,n,x,H);
for (int i=0;i<=3;i++)
    cout<<H[i]<<"  ";
    cout<<"\n";
}

void hessian(double (*pf)(double *x),int n,double *point,double *H)
{
    int i,j;
    double h=1e-3;
    double *temp=new double[n];
    for(i=1;i<=n;i++)
        temp[i-1]=point[i-1];
        for(i=1;i<=n;i++)
            {
                temp[i-1]+=h/2;
                H[(i-1)*n+i-1]=16*pf(temp);
                temp[i-1]-=h;
                H[(i-1)*n+i-1]+=16*pf(temp);
                temp[i-1]+=(3*h/2);
                H[(i-1)*n+i-1]-=pf(temp);
                temp[i-1]-=(2*h);
                H[(i-1)*n+i-1]-=pf(temp);
                H[(i-1)*n+i-1]-=(30*pf(point));
                H[(i-1)*n+i-1]/=(3*h*h);
                temp[i-1]=point[i-1];
            }
            for(i=1;i<=n;i++)
                for(j=i+1;j<=n;j++)
                    {
                        temp[i-1]+=h/2;
                        temp[j-1]+=h/2;
                        H[(i-1)*n+j-1]=16*pf(temp);
                        temp[i-1]-=h;
                        temp[j-1]-=h;
                        H[(i-1)*n+j-1]+=16*pf(temp);
                        temp[i-1]+=(3*h/2);
                        temp[j-1]+=(3*h/2);
                        H[(i-1)*n+j-1]-=pf(temp);
                        temp[i-1]-=(2*h);
                        temp[j-1]-=(2*h);
                        H[(i-1)*n+j-1]-=pf(temp);
                        H[(i-1)*n+j-1]-=(30*pf(point));
                        H[(i-1)*n+j-1]/=(3*h*h);
                        H[(i-1)*n+j-1]-=H[(i-1)*n+i-1];
                        H[(i-1)*n+j-1]-=H[(j-1)*n+j-1];
                        H[(i-1)*n+j-1]/=2;
                        temp[i-1]=point[i-1];
                        temp[j-1]=point[j-1];
                    }
                    for(i=1;i<=n;i++)
                        for(j=i+1;j<=n;j++)
                            H[(j-1)*n+i-1]=H[(i-1)*n+j-1];
                            delete[] temp;
                        }
                        
                        double fun(double *x)
                        {
                            return x[0]*exp(x[0]-x[1]*x[1]);
                        }