Why this kind transformation of Gauss-Hermite quadrature is wrong?

Question

I am confused with numerical integration on Gauss-Hermite quadrature method. Here is an example of what I want to calculate: $$\int_{-\infty}^{+\infty}(x-a)^2e^{-x^2}dx$$ I can directly use Gauss-Hermite quadrature for $f(x)=(x-a)^2$, for the numerical integration, $I=\sum_i f(x_i)w_i=\sum_i(x_i-a)^2w_i$, I can also use this transformation $y=x-a$, then the integration becomes: $$\int_{-\infty}^{+\infty}y^2e^{-(y+a)^2}dy$$ if I use scipy.integrate.quad to calculate these two quadrature, the result are the same, but if I transform the second integration in the following form: $$\int_{-\infty}^{+\infty}(y^2e^{-a^2-2ay})e^{-y^2}dy$$ if I treat $y^2e^{-a^2-2ay}$ as $f(y)$, then the summation should be $\sum_if(y_i)w_i=\sum_i(y_i^2e^{-2ay_i-y_i^2})w_i$, but the result is wrong, so why can't I do this transformation?

here is the python code and what I get from it:

#!/usr/bin/env python
from math import *
import numpy as np
import numpy.polynomial.hermite as nh
import scipy.integrate as si
x,w=nh.hermgauss(28)
print sum(np.power(x-10,2)*w)/sqrt(pi)
print sum(np.power(x,2)*np.exp(-20*x-100)*w)/sqrt(pi)
print si.quad(lambda x:np.power(x,2)*exp(-np.power(x+10,2)),-np.inf,+np.inf)[0]/sqrt(pi)
print si.quad(lambda x:np.power(x-10,2)*exp(-np.power(x,2)),-np.inf,+np.inf)[0]/sqrt(pi)

here is the result

100.5
0.000187032762349
100.5
100.5

Answer 1

Not enough terms. You are only guaranteed equality with infinitely many terms, and I suspect the size of $e^{-100}$ may be a problem. Here is some code (run it after running your code): it works fine for small $a$, choking around $a=6$.

>>> for a in range(11):
...     print(a)
...     print sum(np.power(x-a,2)*w)/sqrt(pi)
...     print sum(np.power(x,2)*np.exp(-2*a*x-a*a)*w)/sqrt(pi)
... 
0
0.5
0.5
1
1.5
1.5
2
4.5
4.5
3
9.5
9.49999999993
4
16.5
16.4999860521
5
25.5
25.4687260564
6
36.5
33.5074807039
7
49.5
21.4184121374
8
64.5
2.98254069772
9
81.5
0.063331734501
10
100.5
0.000187032762349
>>> 

In addition, looking at the quadrature points,

>>> x
array([-6.59160544, -5.85701464, -5.24328537, -4.69075652, -4.17663674,
   -3.68913424, -3.22111208, -2.76779535, -2.32574984, -1.8923605 ,
   -1.46553726, -1.04353527, -0.62483672, -0.20806738,  0.20806738,
    0.62483672,  1.04353527,  1.46553726,  1.8923605 ,  2.32574984,
    2.76779535,  3.22111208,  3.68913424,  4.17663674,  4.69075652,
    5.24328537,  5.85701464,  6.59160544])

We see that we are only using values of $f$ at points $x_i$ such that $|x_i|<7$. For a function that is nearly 0 on this region, this behaviour is not too surprising.

With the weights,

w
array([  1.14013935e-19,   8.31593795e-16,   6.63943671e-13,
     1.47585317e-10,   1.32568250e-08,   5.85771972e-07,
     1.43455042e-05,   2.10618100e-04,   1.95733129e-03,
     1.19684232e-02,   4.95148893e-02,   1.41394610e-01,
     2.82561391e-01,   3.98604718e-01,   3.98604718e-01,
     2.82561391e-01,   1.41394610e-01,   4.95148893e-02,
     1.19684232e-02,   1.95733129e-03,   2.10618100e-04,
     1.43455042e-05,   5.85771972e-07,   1.32568250e-08,
     1.47585317e-10,   6.63943671e-13,   8.31593795e-16,
     1.14013935e-19])

we have the bound $$\left|\sum_{i=1}^{28} f(x_i) w_i \right| ≤ 28\left(\max_i w_i\right) \sup_{x∈[-6.6,-6.6]} |f(x)| ≤ 10 \sup_{x∈[-6.6,-6.6]} |f(x)| $$

Answer 2

I do not know if the following lines really answer your question, but this is too long for a comment. You do not need any approximation/quadrature formula for computing the integral over $\mathbb{R}$ of a polynomial multiplied by $e^{-x^2}$. Let

$$ I_n = \int_{-\infty}^{+\infty} x^n e^{-x^2}\,dx\tag{1} $$ If $n$ is odd the value of $I_n$ is clearly zero, while through the substitution $x=\sqrt{z}$ we have: $$ I_{2m}= 2\int_{0}^{+\infty}x^{2m} e^{-x^2}\,dx = \int_{0}^{+\infty}z^{m-\frac{1}{2}}e^{-z}\,dz = \Gamma\left(m+\frac{1}{2}\right)=\frac{(2m-1)!!}{2^m}\sqrt{\pi}\tag{2} $$ hence: $$ \int_{-\infty}^{+\infty}(x-a)^2 e^{-x^2}\,dx = I_2+a^2 I_0 = \color{red}{\left(a^2+\frac{1}{2}\right)\sqrt{\pi}}.\tag{3}$$