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

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

上一篇讲到了：

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

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

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

f_n = \int_{-\infty}^{\infty} f(x) \psi_n(x) dx

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

下面是个测试例子，计算 ψ5(x)\psi_5(x) 在 2.52.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;
}

结果如下：

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

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

埃尔米特函数的计算（C++）相关推荐

pandas编写自定义函数计算多个数据列的加和（sum）、使用groupby函数和apply函数聚合计算分组内多个数据列的加和
pandas编写自定义函数计算多个数据列的加和(sum).使用groupby函数和apply函数聚合计算分组内多个数据列的加和目录
R语言dist函数距离计算实战（欧几里得距离、曼哈顿距离）
R语言dist函数距离计算实战(欧几里得距离.曼哈顿距离) 目录 R语言dist函数距离计算实战(欧几里得距离.曼哈顿距离)
Smooth_L1_Loss函数的计算方式
Smooth_L1_Loss函数的计算方式从今天开始,阅读faster rcnn的相关代码,并记录我对faster rcnn中特别的层的理解.本篇主要是对smooth_L1_Loss层进行解读. ...
php手册数组函数,PHP - Manual手册 - 函数参考 - Array 数组函数 - array_diff计算数组的差集...
PHP - Manual手册 - 函数参考 - Array 数组函数 - array_diff计算数组的差集 array_diff (PHP 4 >= 4.0.1, PHP 5) array_d ...
DOM中setTimeout() 方法用于在指定的毫秒数后调用函数或计算表达式。
setTimeout() 方法用于在指定的毫秒数后调用函数或计算表达式. <html> <head> <script type="text/javascript ...
整数阶贝塞尔函数c语言,整数阶复宗量变形贝塞尔函数的计算.pdf
整数阶复宗量变形贝塞尔函数的计算.pdf 焦作工学院学报(自然科学版),第卷,第期, 年月 20 2 2001 3 ( ), , , JOurnaI Of JiaOzuO Institute O ...
用位组函数来计算每个月中用户访问网页的天数。
mysql中BIT_COUNT的统计使用下面的例子显示了如何使用位组函数来计算每个月中用户访问网页的天数. CREATE TABLE t1 ( year YEAR ( 4 ), mont ...
php日期相减函数,倒计时函数_计算两个时间相差值_PHP函数
**PHP倒计时函数.求两个日期时间之间相差的时间函数.计算时差函数_PHP函数笔记** ```php /** * 求两个日期时间之间相差的时间 * (针对1970年1月1日之后,求之前可以采用泰勒公 ...
JavaScript eval() 函数，计算某个字符串，并执行其中的的 JavaScript 代码。
JavaScript eval() 函数,计算某个字符串,并执行其中的的 JavaScript 代码. 适合用于计算器的计算,等. 例子: eval("x=10;y=20;document. ...

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

埃尔米特多项式的计算

埃尔米特函数的计算（C++）相关推荐

最新文章

热门文章