Root finding & Integration

←Previous

Homework assignment

Use the skeleton code in Integrator.py, Functions.py and Error.py to implement the Newton-Cotes, Trapezoidal and Simpson rules discussed in the lecture and the Bessel functions J_n(x) in a series and an integral representation.

The assignment consists of the following parts:

1. Implement the integration methods from the lecture with a fixed number of panels in integrate_Nsteps(self,x0,x1,N). To this end, fill in the following methods (in Integrator.py):

   • NCStep: \[I_f^{(a,b)}(N) = \Delta x\sum\limits_{i=0}^{N-1} f(x_i)\,,\]
   • TrapezStep (note the starting point for the summation!): \[I_f^{(a,b)}(N) = \Delta x\left[\frac{f(a)+f(b)}{2}+\sum\limits_{i=i}^{N-1} f(x_i) \right] = \frac{\Delta x}{2} \left[\vphantom{\sum\limits_{i=1}^N}f(a)+2f(x_1)+2f(x_2)+\dots + 2f(x_{N-1})+f(b)\right]\,,\]
   • SimpsonStep: \[I_f^{(a,b)}(N) = \frac{\Delta x}{3} \left[\vphantom{\sum\limits_{i=1}^N}f(a)+4f(x_1)+2f(x_2)+ 4f(x_3)+\dots + 2f(x_{N-2})+4f(x_{N-1})+f(b)\right]\,.\]

   where in all cases $x_i = a+i\Delta x$ and $\Delta x = (b-a)/N$.

   For testing purposes you may take as "kernel" the function $x^n$, implemented as xtoN in the Functions.py file.

2. Use these methods to implement an integration, which runs until a predefined accuracy goal has been reached, in the function
   integrate_Accuracy(self,x0,x1,Ninit,accuracy) in Integrator.py.

   Start with a fixed number of panels, N and store the result for this number of panels $I_f^{(a,b)}(N)$ and for twice the number of panels, $I_f^{(a,b)}(2N)$. Keep doubling the number of panels $N$, $N\to 2N$, until the error function (imported from Error.py) applied to two consecutive results yields a result smaller than the accuracy goal.

   Hint: You may want to reuse the method integrate_Nsteps(self,x0,x1,N).

3. Implement a numerical treatment of the Bessel functions $J_n(x)$ in two ways:

   • Through a series expansion, implemented as BesselJ(order,x,accuracy) in Functions.py: (typo fixed 2010-02-17 21:10)

     \[J_\nu(x) = \sum\limits_{m=0}^\infty \frac{(-1)^m}{m!\Gamma(m+\nu+1)}\, \left(\frac{x}{2}\right)^{\nu+2m}\,.\]

     Here you may assume that $\nu$ always is an integer number, such that $\Gamma(m+\nu+1)$ simplifies to \[\Gamma(m+\nu+1) = (m+\nu)!\,.\]

     Use the accuracy statement together with the error function to terminate the series after the first $k$ terms.

   • Through an integral representation ($n\in {\bf N}$), implemented as BesselJ_From_Int, similar to the original test function xtoN:

     \[\begin{eqnarray*} J_{2n}(x) &=& \frac{1}{\pi}\,\int\limits_{0}^{\pi}\, {\rm d}t\,\cos\left(x\sin t\right)\,\cos\left(2n\,t\right) \\ J_{2n+1}(x) &=& \frac{1}{\pi}\int\limits_{0}^{\pi}\, {\rm d}t\,\sin\left(x\sin t\right)\,\sin\left((2n+1)\,t\right)\,, \end{eqnarray*}\]

   Hint: These Bessel functions and other special functions are provided in scipy. You may use these implementations to check your results, but, of course, you cannot use them as part of the homework solution.

   Hint: The Gamma-function has been added, again by a series and an integral representation, to see how this should work for the case of the Bessel function. Note, however, that the convergence behaviour of the Γ-function is bad.