埃尔米特函数的计算（C++）

前两篇博客介绍了埃尔米特多项式和埃尔米特函数的基本性质。我研究这些的目的其实是为了解决一个函数逼近问题。也就是我要用埃尔米特函数去逼近另外一个函数。有了前两篇的铺垫，这个工作似乎挺简单。

上一篇讲到了：

f(x)=∑n=0∞fnψn(x)

其中：

fn=∫∞−∞f(x)ψn(x)dx

所以这个函数逼近问题其实就是一个积分问题。但是有个前提，我们要能算出ψn(x)。这一篇就来写写如何用 C语言来实现这个 ψn(x)。

埃尔米特多项式的计算

计算埃尔米特函数之前，我们先来计算埃尔米特多项式。埃尔米特多项式是普通多项式的特例，所以我们先来写个多项式计算的类。

class Polynomials{public:    Polynomials(unsigned int N = 0);    ~Polynomials() {delete[] m_coeff;}    double coefficient(unsigned int m) const;    void setCoefficient(unsigned int m, double c);    void setCoefficient(double coeff[]);    /**     * @brief eval 计算多项式在 x 处的值     * @param x     * @return 返回多项式在 x 处的值     */    double eval(double x) const;    /**     * @brief eval_d 计算多项式在 x 处的值和导数值     * @param x     * @param dp [out] 返回导数值     * @return 返回多项式在 x 处的值     */    double eval_d(double x, double &dp) const;    /**     * @brief eval_dd 计算多项式在 x 处的第 0 阶直到第 M 阶导数值     * @param x     * @param dp [out] 返回第 0 阶直到第 M 阶导数值     * @param M     */    void eval_dd(double x, double dp[], int M) const;    void trim(double epsAbs);    void debug(char s);protected:    unsigned int m_N;    double *m_coeff;};

具体的实现代码如下：

inline void zeros(double arr[], int N){    while(N-- != 0)    {        *arr++ = 0;    }}Polynomials::Polynomials(unsigned int N){    m_N = N;    m_coeff = new double [N + 1];    zeros(m_coeff, N+1);}double Polynomials::coefficient(unsigned int m) const{    return m_coeff[m];}double Polynomials::eval(double x) const{    return poly(m_coeff, m_N, x);}double Polynomials::eval_d(double x, double &dp) const{    return dpoly(m_coeff, m_N, x, dp);}void Polynomials::eval_dd(double x, double dp[], int M) const{    ddpoly(m_coeff, m_N, x, dp, M);}void Polynomials::setCoefficient(unsigned int m, double c){    if(m <= m_N)    {        m_coeff[m] = c;    }}void Polynomials::setCoefficient(double coeff[]){    for(unsigned int i = 0; i <= m_N; i++)    {        m_coeff[i] = coeff[i];    }}void Polynomials::trim(double epsAbs){    epsAbs = fabs(epsAbs);    for(unsigned int i = 0; i <= m_N; i++)    {        if(fabs(m_coeff[i]) < epsAbs)        {            m_coeff[i] = 0;        }    }}void Polynomials::debug(char s){    bool start = true;    if(m_coeff[0] != 0)    {        std::cout << m_coeff[0];        start = false;    }    if(m_N >= 1 && m_coeff[1] != 0)    {        if(start)        {            std::cout << m_coeff[1] << s;            start = false;        }        else        {            std::cout << "+" << m_coeff[1] << s;        }    }    for(unsigned int i = 2; i < m_N + 1; i++)    {        if(m_coeff[i] > 0)        {            if(start)            {                std::cout << m_coeff[i] << s << "^" << i;                start = false;            }            else            {                std::cout << "+" << m_coeff[i] << s << "^" << i;            }        }        else if(m_coeff[i] < 0)        {            std::cout << m_coeff[i] << s << "^" << i;            start = false;        }    }    std::cout << std::endl;}

上面的代码用到了 poly()、dpoly()、ddpoly()函数，这三个函数在我另一篇博客中有介绍：
多项式求值的几个简单算法（C语言）

有了这个基础代码只有就可以实现埃尔米特多项式了：

//物理学家的埃尔米特多项式class HermitePolynomials : public Polynomials{public:    HermitePolynomials(unsigned int N = 0);private:    void buildCoeff();};

代码实现如下，用到了埃尔米特多项式的递推公式：

HermitePolynomials::HermitePolynomials(unsigned int N)    :Polynomials(N){    switch (N)    {    case 0:        setCoefficient(0, 1.0);        break;    case 1:        setCoefficient(0, 0.0);        setCoefficient(1, 2.0);    default:        buildCoeff();        break;    }}void HermitePolynomials::buildCoeff(){    double *p1 = new double [m_N + 1];    double *p2 = new double [m_N + 1];    double *p = new double[m_N + 1];    for(unsigned int i = 0; i <= m_N; i++)    {        p[i] = 0.0;        p1[i] = 0.0;        p2[i] = 0.0;    }    p[1] = 2.0;    p1[0] = 1.0;    for(unsigned int n = 2; n <= m_N; n++)    {        for(unsigned int j = 0; j <= n; j++)        {            p2[j] = p1[j];            p1[j] = p[j];        }        //H_n(x) = 2x H_{n-1}(x) - 2(n-1)H_{n-2}(x)        p[0] = -2.0 * (n - 1) * p2[0];        for(unsigned int j = 1; j <= n; j++)        {            p[j] = 2 * p1[j - 1] - 2.0 * (n - 1) * p2[j];        }    }    setCoefficient(p);    trim(0.5);    delete[] p;    delete[] p1;    delete[] p2;}

下面是个测试代码：

void test1(){    HermitePolynomials *p[11];    for(int i = 0; i <= 10; i ++)    {        p[i] = new HermitePolynomials(i);        p[i]->debug('x');    }}

程序输出如下：

12x-2+4x^2-12x+8x^312-48x^2+16x^4120x-160x^3+32x^5-120+720x^2-480x^4+64x^6-1680x+3360x^3-1344x^5+128x^71680-13440x^2+13440x^4-3584x^6+256x^830240x-80640x^3+48384x^5-9216x^7+512x^9-30240+302400x^2-403200x^4+161280x^6-23040x^8+1024x^10

经验算，这些值都是正确的。

最后是在埃尔米特多项式的基础上构造埃尔米特函数：

class HermiteFunction{public:    HermiteFunction(unsigned int N);    ~HermiteFunction();    /**     * @brief eval_exp 计算带权重、归一化的函数的值     * @param x     * @return 函数值     */    double eval_exp(double x) const;    /**     * @brief eval_d_exp 计算函数值和导数值     * @param x     * @param dp [out] 导数值     * @return 函数值     */    double eval_d_exp(double x, double &dp) const;    /**     * @brief eval_d5_exp 计算第 0 阶（也就是函数值）到第 5 阶导数值     * @param x     * @param dp [out] dp[0] .. dp[5] 分别存放第 0 阶（也就是函数值）到第 5 阶导数值     */    void eval_d5_exp(double x, double dp[6]) const;private:    HermitePolynomials *m_poly;    double m_norm;};

其中 eval_d5_exp() 函数是用来计算埃尔米特函数的前 5 阶导数的。因为我做的问题中需要求这 5 阶导数，所以就实现了这么个函数。普通的应用应该用不到。

下面是代码实现：

inline double factorial(int n){    double f = 1;    for(int i = 2; i <= n; i++)    {        f = f * i;    }    return f;}inline double pow_int(double x, int n){    double ret = 1;    while(n)    {        if(n & 1) ret = ret * x;        x = x * x;        n = (unsigned int)n >> 1;    }    return ret;}HermiteFunction::HermiteFunction(unsigned int N){   m_poly = new HermitePolynomials(N);   m_norm = SQRT_PI * pow_int(2.0, N) * factorial(N);   m_norm = 1.0 / sqrt(m_norm);}HermiteFunction::~HermiteFunction(){    delete m_poly;}double HermiteFunction::eval_exp(double x) const{    return m_norm * m_poly->eval(x) * exp(- x * x / 2.0);}double HermiteFunction::eval_d_exp(double x, double &dp) const{    double y = m_poly->eval_d(x, dp);    double e = exp(- x * x / 2.0);    dp = m_norm * (dp - y * x) * exp(- x * x / 2.0);    return m_norm * y * e;}void HermiteFunction::eval_d5_exp(double x, double dp[6]) const{    double dd[6];    zeros(dd, 6);    double x2 = x * x;    double x3 = x2 * x;    double x4 = x2 * x2;    double x5 = x2 * x3;    double e = exp(- x2 / 2.0);    m_poly->eval_dd(x, dd, 5);    dp[0] = m_norm * dd[0] * e;    dp[1] = m_norm * (- x * dd[0] + dd[1]) * e;    dp[2] = m_norm * ((x2 - 1) * dd[0] - 2 * x * dd[1] + dd[2]) * e;    dp[3] = m_norm * (-x * (x2 - 3) * dd[0] + (3 * x2 - 3) * dd[1] -3 * x * dd[2] + dd[3]) * e;    dp[4] = m_norm * ((3 - 6 * x2 + x4) * dd[0] + (12 * x - 4 * x3) * dd[1] + (-6 + 6 * x2) * dd[2] - 4 * x * dd[3] + dd[4]) * e;    dp[5] = m_norm * ((-15 * x + 10 * x3 - x5) * dd[0]                     + (15 - 30 * x2 + 5 * x4) * dd[1]                     + (30 * x -10 * x3) * dd[2]                     + (-10 + 10 * x2) * dd[3] -5 * x * dd[4] + dd[5]) * e;}

这里面又涉及到些公式推导，手工推导很麻烦，我是用 Mathematica 计算的。另外，Mathematica 中有HermiteH[] 函数。我的代码验证工作也大量的使用了 Mathematica。

下面是个测试例子，计算 ψ5(x) 在 2.5 处的函数值和直到 5 阶导数的值。

void test2(){    HermiteFunction p5(5);   double dp[6];   p5.eval_d5_exp(2.5, dp);   std::cout << dp[0] << std::endl;   std::cout << dp[1] << std::endl;   std::cout << dp[2] << std::endl;   std::cout << dp[3] << std::endl;   std::cout << dp[4] << std::endl;   std::cout << dp[5] << std::endl;}

结果如下：

0.4926270.563193-2.33998-0.21202917.7321-30.7134

这个结果与 Mathematica 的计算结果比对过，是正确的。

好了，到此，埃尔米特函数的计算问题就基本完成了。代码写的比较仓卒，没有特意的优化。希望我的代码对那些对数值计算感兴趣的同学能有所帮助。