Gaussian Quadrature strategy for $\int\limits_0^s {f\left( x \right){e^{ – ax}}} dx$

approximationintegrationnumerical methods

Are there any suitable Gaussian quadrature strategy for for this type of integral

$\int\limits_0^s {f\left( x \right){e^{ – ax}}} dx$ ?

Where $s$ is a positive real number.

I know that in that case of the interval from $0$ to ${ + \infty }$ we could use a Laguerre-Gauss Quadrature like this $\int_0^{ + \infty } {{e^{ – ax}}} f(x)dx \approx \frac{1}{a}\sum\limits_{i = 1}^n {{w_i}} f\left( {\frac{{{x_i}}}{a}} \right)$

But for my case, I could not come up with any answer.

Please help me, thank you very much !

Best Answer

In principle you can use the Gram-Schmidt procedure to derive the orthogonal polynomials - which will result in modified Laguerre polynomials, however these quickly become unwieldy $$ p_0(x)=1\\ p_1(x) = x - 1+\frac{s e^{-s}}{1-e^{-s}}\\ p_2(x) = x^2 - Ap_1(x) - B $$ where $$ A = \frac{\frac{s^{3} e^{s}}{e^{\left(2 \, s\right)} - e^{s}} + \frac{2 \, s^{2} e^{s}}{e^{\left(2 \, s\right)} - e^{s}} + \frac{4 \, s e^{s}}{e^{\left(2 \, s\right)} - e^{s}} - \frac{2 \, s}{e^{\left(2 \, s\right)} - e^{s}} - \frac{2 \, s}{e^{s} - 1} + \frac{4 \, e^{s}}{e^{\left(2 \, s\right)} - e^{s}} - \frac{4 \, e^{s}}{e^{s} - 1} - \frac{4}{e^{\left(2 \, s\right)} - e^{s}} + \frac{4}{e^{s} - 1}}{\frac{s^{2} e^{\left(2 \, s\right)}}{e^{\left(3 \, s\right)} - 2 \, e^{\left(2 \, s\right)} + e^{s}} - \frac{s^{2}}{e^{\left(2 \, s\right)} - 2 \, e^{s} + 1} + \frac{e^{\left(2 \, s\right)}}{e^{\left(3 \, s\right)} - 2 \, e^{\left(2 \, s\right)} + e^{s}} - \frac{e^{\left(2 \, s\right)}}{e^{\left(2 \, s\right)} - 2 \, e^{s} + 1} - \frac{2 \, e^{s}}{e^{\left(3 \, s\right)} - 2 \, e^{\left(2 \, s\right)} + e^{s}} + \frac{2 \, e^{s}}{e^{\left(2 \, s\right)} - 2 \, e^{s} + 1} + \frac{1}{e^{\left(3 \, s\right)} - 2 \, e^{\left(2 \, s\right)} + e^{s}} - \frac{1}{e^{\left(2 \, s\right)} - 2 \, e^{s} + 1}}\\ B = \frac{s^{2} e^{-s} + 2s e^{-s}}{e^{-s} - 1}+2 $$

By fixing $s$ this process can be done semi-numerically with a computer algebra system like Sagemath or Mathematica. For instance, for $s=0.4$ the polynomials are $$ p_0(x)=1\\ p_1(x)= x-0.186702087312105\\ p_2(x)=x^2 - 0.389329291486829x + 0.0246035818232814 $$ which has roots $x_i\in [0.0793791632266308,0.309950128260198]$.

The next step is to solve the tower equations to determine the weights $w_i$ $$ \sum_{i=0}^2 w_ip_n(x_i)=\delta_{n,0}0.3296799539643607. $$ The number multiplying the delta is the norm of $p_0(x)$. This leads to the matrix equation $$ \left(\begin{array}{rr} 1 & 1 \\ -0.107322924085475 & 0.123248040948093 \end{array}\right)\left(\begin{array}{c} w_0\\ w_1 \end{array}\right) =\left(\begin{array}{c} 0.3296799539643607 \\ 0 \end{array}\right) $$ which has solution $w_i\in \left[0.1762251741456236, 0.153454779818737\right]$.

I've checked this by writing a small routine in Sagemath

def gqexp(f):
   xi = vector([0.0793791632266308,0.309950128260198])
   wi = vector([0.1762251741456236,0.153454779818737])
   fncall = vector([f(x) for x in xi])
   return expand(wi.dot_product(fncall))
Related Question