Numerical integration for integral with weight $e^{-x^2/2}$

356 Views Asked by At

The task is to designate $a, b, x$ and $y$ such that the $Q$ is precise numerical integration for polynomials of the degree $\leq3$ $$Q(f)=af(x)+bf(y)$$ $$\int_{0}^{\infty} f(x)e^{-x^2/2} dx$$ Of course i can create a set of equations by substituting $f(x)$ with $1, x, x^2$ and $x^3$ and calculate it, which by the way I have already done with a little help of wolframalpha, but I know that it can be done in a tricky way using interpolation polynomial, maybe it has something to do with Hermites polynomial? Will be glad to get any way to solve it other than calculating this set of equations! :)

Mistakably i calculated it wrong, the result is very very ugly!

2

There are 2 best solutions below

0
On BEST ANSWER

I don't see a simple relationship with Gauss-Hermite quadrature because the zeros of the Hermite polynomials are rational, but here the roots are not. The orthogonal (monic) polynomials are $$\begin{align}\pi_0(x)&=1\\ \pi_1(x)&=x-\sqrt\frac2{\pi}\\ \pi_2(x)&=x^2-\frac{\sqrt{\frac{\pi}2}}{\frac{\pi}2-1}x+\frac{2-\frac{\pi}2}{\frac{\pi}2-1}\end{align}$$ We can solve by the quadratic formula to get $\pi_2(x_1)=\pi_2(x_2)=0$ for the nodes, and since $$\int_0^{\infty}(\pi_1(x))^2e^{-x^2/2}dx=\sqrt{\frac{\pi}2}-\sqrt{\frac2{\pi}}$$ for monic polynomials the formula for the weights gives us $$w_i=\frac{\sqrt{\frac{\pi}2}-\sqrt{\frac2{\pi}}}{\left(2x_i-\frac{\sqrt{\frac{\pi}2}}{\frac{\pi}2-1}\right)\left(x_i-\sqrt\frac2{\pi}\right)}$$ Since the weights end up looking ugly, I just wrote up a program that conputes the nodes and weights and then checks the formula for polynomials of degree $\le3$.

    program hermite2
   use ISO_FORTRAN_ENV
   implicit none
   integer, parameter :: wp = REAL64
   real(wp) a10, a21, a20
   real(wp), parameter :: pi = 4*atan(1.0_wp)
   real(wp) x1, x2
   real(wp) w1, w2

   a10 = -sqrt(2/pi)
   a21 = -sqrt(pi/2)/(pi/2-1)
   a20 = (2-pi/2)/(pi/2-1)
   x1 = (-a21-sqrt(a21**2-4*a20))/2
   x2 = (-a21+sqrt(a21**2-4*a20))/2
   w1 = (a10-1/a10)/(2*x1+a21)/(x1+a10)
   w2 = (a10-1/a10)/(2*x2+a21)/(x2+a10)
   write(*,*) 'x1 = ', x1
   write(*,*) 'w1 = ', w1
   write(*,*) 'x2 = ', x2
   write(*,*) 'w2 = ', w2
   write(*,*) w1+w2,sqrt(pi/2)
   write(*,*) w1*x1+w2*x2,1.0_wp
   write(*,*) w1*x1**2+w2*x2**2,sqrt(pi/2)
   write(*,*) w1*x1**3+w2*x2**3,2.0_wp
end program hermite2

Output was:

 x1 =   0.424538328648333
 w1 =   0.905845053005361
 x2 =    1.77119082811243
 w2 =   0.347469084310139
   1.25331413731550        1.25331413731550
   1.00000000000000        1.00000000000000
   1.25331413731550        1.25331413731550
   2.00000000000000        2.00000000000000
2
On

See for example the tables cited in Section 2 of http://vixra.org/abs/1303.0013