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

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

上一篇讲到了：

$$f(x) = \sum_{n=0}^{\infty} f_n \psi_n(x)$$
其中：
$$f_n = \int_{-\infty}^{\infty} f(x) \psi_n(x) dx$$

所以这个函数逼近问题其实就是一个积分问题。但是有个前提，我们要能算出$\psi_n(x)$。这一篇就来写写如何用 C语言来实现这个 $\psi_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');
    }
}

程序输出如下：

1
2x
-2+4x^2
-12x+8x^3
12-48x^2+16x^4
120x-160x^3+32x^5
-120+720x^2-480x^4+64x^6
-1680x+3360x^3-1344x^5+128x^7
1680-13440x^2+13440x^4-3584x^6+256x^8
30240x-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。

下面是个测试例子，计算 $\psi_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.492627
0.563193
-2.33998
-0.212029
17.7321
-30.7134

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

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