User:Egm6322.s14.team2/report3

Report 3 PEA 2 Spring 2014 Team 2

Problems are given in lecture notes

R3.1:

Problem Statement

sec.67a, Pb-67.1 (equilibrium of line defects, Legendre diff. eq., Legendre polynomials)

Provide a detailed justification for the following relations:

$S(u){\text{behaves like }}(u-u_{i})^{-1}$

(Eq 1-1)

and thus,

$P^{''}(u)\to 0{\text{ as }}u\to u_{i}$

(Eq 1-1)

$c=N(N-1)$

(Eq 1-2)

$P(u)$ satisfies the differential equation,

$(1-u^{2})P^{''}+N(N-1)P=0$

(Eq 1-3)

$(1-u^{2})Q^{''}-4uQ^{'}+{N(N-1)-2}Q=0$

(Eq 1-4)

Find

(a): (b):

Solution

On our honor, we did this assignment on our own, without looking at the solutions in previous semesters or other online solutions.

Given at equilibrium,

$\sum _{j=1,j{\not =}i}^{N}{\dfrac {1}{u_{i}-u_{j}}}=0$

(Eq 1-5)

(Eq 1-5) is the definition of S(u), a summation of N terms.

$S(u)=\sum _{i=1}^{N}{\dfrac {1}{u-u_{i}}}$

(Eq 1-6)

Substituting (Eq 1-5)for the N-1 terms (excluding j=i) in (Eq 1-6) and evaluating the behavior of the function as u approaches u_i results in (Eq 1-7).

$S(u)\sim {\dfrac {1}{u-u_{i}}}$

(Eq 1-7)

P(u) is defined in equation (Eq 1-8) and the derivative is displayed in (Eq 1-9).

$P(u)=\prod _{j=1}^{N}(u-u_{j})$

(Eq 1-8)

$P'(u)=\sum _{i=1}^{N}\prod _{j=1,j\neq {i}}^{N}(u-u_{j})$

(Eq 1-9)

Substituting (Eq 1-8) and (Eq 1-9) into (Eq 1-6), (Eq 1-10) was established.

$S(u)={\dfrac {P'(u)}{P(u)}}$

(Eq 1-10)

Using the chain rule the following equation can be defined using (Eq 1-10).

${\dfrac {P''(u)}{P'(u)}}=S(u)+{\dfrac {S'(u)}{S(u)}}$

(Eq 1-11)

(Eq 1-12) is the derivative of S as $u\to u_{i}$.

$S'(u)\sim {\dfrac {1}{(u-u_{i})^{2}}}$

(Eq 1-12)

Substituting (Eq 1-8) and (Eq 1-9) into (Eq 1-12), multiplying the last term by (u-u_i)/(u-u_i), and evaluating the behavior of the system as $u\to u_{i}$ yields P(u)/P'(u) approaches 0 as $u\to u_{i}$.

${\dfrac {P''(u)}{P'(u)}}\sim {{\dfrac {1}{u-u_{i}}}-{\dfrac {1}{u-u_{i}}}}$

(Eq 1-13)

Evaluating the product in (Eq 1-8), a power series with constant coefficients that are products of the u_i positions can be expressed as follows:

$P(u)=\sum _{n=0}^{N}c_{n}u^{n}$

(Eq 1-14)

And the first and second derivatives can also be expressed as a power series as follows:

$P'(u)=\sum _{n=1}^{N}nc_{n}u^{n-1}$

(Eq 1-15)

$P''(u)=\sum _{n=2}^{N}n(n-1)c_{n}u^{n-2}$

(Eq 1-16)

Given that P(u) is proportional to ${\frac {P(u)}{u^{2}-1}}$,

$P''(u)=c(u^{2}-1)^{-1}P(u)$

(Eq 1-17)

Substituting the power series expressions of the solutions into the above equation:

$\sum _{n=2}^{N}c_{n}n(n-1)u^{n-2}={\frac {c\sum _{n=0}^{N}c_{n}u^{n}}{(u+1)(u-1)}}$

(Eq 1-18)

Considering just the final term in the summation:

$c_{N}N(N-1)=cc_{N}$

(Eq 1-19)

The above equation simplifies to c=N(N-1). Substituting the coefficient, c, back into (Eq 1-17)

$(1-u^{2})P''(u)+N(N-1)P(u)$

(Eq 1-20)

Now, using the transformation for P as a function of Q, and the first and second derivatives:

$P(u)=(u^{2}-1)Q(u)$

(Eq 1-21)

$P'(u)=2uQ(u)+(u^{2}-1)Q'(u)$

(Eq 1-22)

$P''(u)=2Q(u)+4uQ'(u)+(u^{2}-1)Q''(u)$

(Eq 1-23)

Substituting the transformation into the differential equation in (Eq 1-20)

$(1-u^{2})Q''(u)-4uQ(u)+(N(N-1)-2)Q(u)=0$

(Eq 1-24)

R3.2: Hermite Differential Equation Form of Line Defect Differential Equation

Problem Statement

sec.67b, Pb-67.2 (equilibrium of line defects, Hermite diff. eq., Hermite polynomials)

Equation (Eq 2-1) is given:

$P^{''}(u)-uP^{'}(u)+NP(u)=0$

(Eq 2-1)

Find a suitable change of variable to show that (Eq 2-1) ((6) p.68-5) can be written as the Hermite differential equation of the form:

$y^{''}(x)-2xy^{'}(x)+2ny(x)=0$

(Eq 2-2)

For the 2 cases of $N=7,8$ particles, find the non-dimensional equilibrium positions and plot them.

Find

(a): A suitable change of variable to show that (Eq 2-1) can be written as the Hermite differential equation of the form of (Eq 2-2)

Solution

On our honor, we did this assignment on our own, without looking at the solutions in previous semesters or other online solutions.

Part (a): Find Suitable Change of Variable

For this problem we seek to determine a suitable change of variables that will give us (Eq 2-2) from (Eq 2-1). The variable of interest is changing from u to x in this problem. We need to find some sort of relation between the two variables that will gives us what we want. We can observe the differences between the two forms of the differential equations and conclude that a factor of 2 on the zeroth and first derivative terms is the only difference.

Plot Solution for 2 Cases

Here we are asked to plot the solution to the Hermite differential equation for the cases when $N=7,8$.

Case for $N=7$

Here the Hermite differential equation takes the following form.

$y^{''}(x)-2xy^{'}(x)+14y(x)=0$

(Eq 2-2)

Case for $N=8$

Here the Hermite differential equation takes the following form.

$y^{''}(x)-2xy^{'}(x)+16y(x)=0$

(Eq 2-2)

R3.3:

Problem Statement

sec.70a, Pb-70.1 (Kepler’s equation, acceleration vector in planar motion using Euler formula)

Show,

$\mathbf {a} ={\frac {d}{dt}}\,\mathbf {v} =[{\ddot {r}}-r({\dot {\theta }})^{2}]\mathbf {e} _{r}+[{\frac {1}{r}}\,{\frac {d}{dt}}\,(r^{2}\,{\dot {\theta }})]\mathbf {e} _{\theta }$

(Eq 2-1)

Find

(a): (b):

Solution

On our honor, we did this assignment on our own, without looking at the solutions in previous semesters or other online solutions.

Using polar coordinates, the unit vectors that form a basis in the orbital plane and their first order time derivatives are:

$\mathbf {e_{r}} =e^{i{\theta }}$

(Eq 2-2)

$\mathbf {e_{\theta }} =ie^{i{\theta }}$

(Eq 2-3)

$\mathbf {\dot {e_{r}}} =i{\dot {\theta }}e^{i{\theta }}={\dot {\theta }}\mathbf {e_{\theta }}$

(Eq 2-4)

$\mathbf {\dot {e_{\theta }}} =-{\dot {\theta }}e^{i{\theta }}=-{\dot {\theta }}\mathbf {e_{r}}$

(Eq 2-5)

The position vector, r, in this basis is:

$\mathbf {r} =r\mathbf {e_{r}}$

(Eq 2-6)

The time derivative of the position vector is the velocity:

$\mathbf {v} ={\frac {d}{dt}}\mathbf {r} ={\dot {r}}\mathbf {e_{r}} +r{\dot {\theta }}\mathbf {e_{\theta }}$

(Eq 2-7)

The second derivative with respect to time is the acceleration:

$\mathbf {a} ={\frac {d}{dt}}\mathbf {v} ={\ddot {r}}\mathbf {e_{r}} +{\dot {r}}{\dot {\theta }}\mathbf {e_{\theta }} +{\dot {r}}{\dot {\theta }}\mathbf {e_{\theta }} +r{\ddot {\theta }}\mathbf {e_{\theta }} -r{\dot {\theta }}^{2}\mathbf {e_{r}}$

(Eq 2-8)

After combining like terms and factoring a ${\frac {1}{r}}$ out of the $e_{\theta }$ term, the equation above simplifies to:

$\mathbf {a} =({\ddot {r}}-r{\dot {\theta }}^{2})\mathbf {e_{r}} +{\frac {1}{r}}(2r{\dot {r}}{\dot {\theta }}+r^{2}{\ddot {\theta }})\mathbf {e_{\theta }}$

(Eq 2-9)

Recognizing the derivative of $r^{2}{\dot {\theta }}$ and performing the integration, Kepler's equation takes the form:

$\mathbf {a} =({\ddot {r}}-r{\dot {\theta }}^{2})\mathbf {e_{r}} +{\frac {1}{r}}{\frac {d}{dt}}(r^{2}{\dot {\theta }})\mathbf {e_{\theta }}$

(Eq 2-10)

R3.4: Solutions to First Oder Linear ODE with Variable Coefficients using Euler Integrating Factor

Problem Statement

A general First Order Linear Differential Equation with Varying Coefficients (L1-ODE-VC) has the following form:

$a_{1}(x)y'+a_{0}(x)y=b(x)$

(Eq 3-1)

Find

(a)

$a_{1}(x)=1,a_{0}(x)=x,b(x)=2x+3$

(Eq 3-2)

(b)

$a_{1}(x)=x^{2}+1,a_{0}(x)=x,b(x)=2x$

(Eq 3-3)

Solution

On our honor, we did this assignment on our own, without looking at the solutions in previous semesters or other online solutions.

Part a

The equation in 1 can be written as follows:

$y'+xy=2x+3$

(Eq 3-1)

The general solution for the Euler Integrating Factor Method is:

$y(x)={\frac {1}{h(x)}}\left[\int ^{x}{h(s)b(s)ds+k_{2}}\right]$

(Eq 3-2)

And $h(x)$ is:

$h(x)=e^{\left[\int ^{x}a_{0}(s)ds+k_{1}\right]}$

(Eq 3-3)

($k_{1}$ is arbitrary, it scales the integrating factor by a constant but has no effect on the final solution. Here we set it to 0 for ease of computation.)

Plugging $a_{0}$ into $h(x)$ and integrating gives:

$h(x)=e^{\frac {x^{2}}{2}}$

(Eq 3-4)

Plugging $h(x)$ into the general solution we get:

$y(x)={\frac {1}{e^{x^{2}/2}}}\int ^{x}{\left[e^{s^{2}/2}(2s+3)ds+k_{2}\right]}$

(Eq 3-5)

Integrating this gives the solution:

$y(x)={\frac {1}{e^{\frac {x^{2}}{2}}}}\left[(3{\sqrt {\frac {\pi }{2}}})(erfi({\frac {x}{\sqrt {2}}}))+(2e^{\frac {x^{2}}{2}}+k_{2})\right]$

(Eq 3-6)

Part b

The equation 2 can be rewritten as:

$(x^{2}+1)y'+xy=2x$

(Eq 3-7)

(Eq 3-7) is divided by $(x^{2}+1)$ in order to have the coefficient of the leading derivative be unity.

$y'+{\frac {xy}{x^{2}+1}}={\frac {2x}{x^{2}+1}}$

(Eq 3-8)

$h(x)$ is then found by the following equation:

$h(x)=e^{\left[{\frac {x}{x^{2}+1}}+k_{1}\right]}={\sqrt {x^{2}+1}}$

(Eq 3-9)

(Eq 3-2) then becomes:

$y(x)={\frac {1}{\sqrt {x^{2}+1}}}\left[\int ^{x}{{\sqrt {s^{2}+1}}{\frac {2s}{s^{2}+1}}ds+k_{2}}\right]$

(Eq 3-10)

Which after evaluating the integral can be simplified to:

$2+{\frac {k_{2}}{\sqrt {x^{2}+1}}}$

(Eq 3-11)

R3.5: Commented problems from R2

Problem Statement

Redo commented problems from Report 2.

R2.4:Differential Element in Spherical Coordinates

Given Information

The following equation is given for ds²:

$ds^{2}=\sum _{i=1}^{3}(dx_{i})^{2}=h_{1}^{2}dr^{2}+h_{2}^{2}d{\theta }^{2}+h_{3}^{2}d{\varphi }^{2}$

(Eq 4-1)

x_i is written in spherical coordinates in (Eq 4-2) through (Eq 4-4).

$x_{1}=rcos{\theta }cos{\varphi }$

(Eq 4-2)

$x_{2}=rcos{\theta }sin{\varphi }$

(Eq 4-3)

$x_{3}=rsin{\theta }$

(Eq 4-4)

Problem Statement

Show that the infinitesimal length ds² in (Eq 4-1) can be written in spherical coordinates, (Eq 4-5), as given in Pb.R*7.3 on pp.39-[1,3] of sec.39 of the notes.

$ds^{2}=dr^{2}+r^{2}d{\theta }^{2}+r^{2}cos^{2}{\theta }d{\varphi }^{2}$

(Eq 4-5)

Solution

On our honor, we did this assignment on our own, without looking at the solutions in previous semesters or other online solutions.

The derivatives of x_i are given in (Eq 4-6) through (Eq 4-8).

$dx_{1}=cos{\theta }cos{\varphi }dr-rsin{\theta }cos{\varphi }d{\theta }-rcos{\theta }sin{\varphi }d{\varphi }$

(Eq 4-6)

$dx_{2}=cos{\theta }sin{\varphi }dr-rsin{\theta }sin{\varphi }d{\theta }+rcos{\theta }cos{\varphi }d{\varphi }$

(Eq 4-7)

$dx_{3}=sin{\theta }dr+rcos{\theta }d{\theta }$

(Eq 4-8)

(Eq 4-9) was developed by substituting the squared derivatives in the summation of (Eq 4-1) and combining like terms.

$ds^{2}=(sin^{2}{\theta }(cos^{2}{\varphi }+sin^{2}{\varphi })+cos^{2}{\theta })dr^{2}+(2rsin{\theta }cos{\theta }-2rsin{\theta }cos{\theta }(sin^{2}{\varphi }+cos^{2}{\varphi }))drd{\theta }+(r^{2}(sin^{2}{\theta }(sin^{2}{\varphi }+cos^{2}{\varphi })+cos^{2}{\theta })d{\theta }^{2}+(r^{2}cos^{2}{\theta }(sin^{2}{\varphi }+cos^{2}{\varphi }))d{\varphi }^{2}$

(Eq 4-9)

Using the Pythagorean Trigonometric Identity, (Eq 4-9) reduces to (Eq 4-10).

$ds^{2}=dr^{2}+r^{2}d{\theta }^{2}+r^{2}cos^{2}d{\varphi }^{2}$

(Eq 4-10)

The magnitude of the tangent vector h are then identified as follows:

$h_{1}=1$

(Eq 4-11)

$h_{2}=r$

(Eq 4-12)

$h_{3}=rcos{\theta }$

(Eq 4-13)

The definition of the Laplace operator, ${\triangle }u$, in cartesian coordinates:

${\triangle }u=\sum _{i=1}^{3}{\frac {{\partial }u}{{\partial }x_{i}}}$

(Eq 4-14)

The Laplace operator can be expressed in spherical coordinates from equation (2) from the class notes 39-2,

${\triangle }u={\frac {1}{h_{1}h_{2}h_{3}}}({\frac {\partial }{{\partial }r}}({\frac {h_{1}h_{2}h_{3}}{h_{1}^{2}}}{\frac {{\partial }u}{{\partial }r}})+{\frac {\partial }{{\partial }{\theta }}}({\frac {h_{1}h_{2}h_{3}}{h_{2}^{2}}}{\frac {{\partial }u}{{\partial }{\theta }}})+{\frac {\partial }{{\partial }{\varphi }}}({\frac {h_{1}h_{2}h_{3}}{h_{3}^{2}}}{\frac {{\partial }u}{{\partial }{\varphi }}}))$

(Eq 4-15)

After substituting the values of h from (Eq 4-11) through (Eq 4-13), the Laplacian takes the form:

${\triangle }u={\frac {1}{r^{2}cos{\theta }}}({\frac {\partial }{{\partial }r}}(r^{2}cos{\theta }{\frac {{\partial }u}{{\partial }r}})+{\frac {\partial }{{\partial }{\theta }}}(cos{\theta }{\frac {{\partial }u}{{\partial }{\theta }}})+{\frac {\partial }{{\partial }{\varphi }}}({\frac {1}{cos{\theta }}}{\frac {{\partial }u}{{\partial }{\varphi }}}))$

(Eq 4-16)

Factoring the $cos{\theta }$ from the first term, the partial derivative with respect to r, and the $sec{\theta }$ from the third term, the partial derivative with respect to ${\varphi }$, since both values are fixed in the partial derivatives, the Laplacian takes the form:

${\triangle }u={\frac {1}{r^{2}}}{\frac {\partial }{{\partial }r}}(r^{2}{\frac {{\partial }u}{{\partial }r}})+{\frac {1}{r^{2}cos{\theta }}}{\frac {\partial }{{\partial }{\theta }}}(cos{\theta }{\frac {{\partial }u}{{\partial }{\theta }}})+{\frac {1}{r^{2}cos^{2}{\theta }}}{\frac {{\partial }^{2}u}{{\partial }{{\varphi }^{2}}}}$

(Eq 4-17)

R2.5: Using Variation of Parameters To Find Second Solution To The Legendre Differential Equation

Given Information

The following equation is given for P₂(x), the Legendre Polynomial of degree n=2:

$P_{2}(x)={\dfrac {1}{2}}(3x^{2}-1)$

(Eq 5-1)

Problem Statement

Using variation of parameters, show that Equation 5-2 is a second solution to the Legendre differential equation in Equation 5-3 as given in Pb.R*6.11 on pp.37-4 in sec.37 of the notes.

$Q_{2}(x)={\dfrac {1}{4}}(3x^{2}-1)log({\dfrac {1+x}{1-x}})-{\dfrac {3}{2}}x$

(Eq 5-2)

Solution

On our honor, we did this assignment on our own, without looking at the solutions in previous semesters or other online solutions.

The Legendre differential equation of degree n=2 is given in equation 5-3.

$(1-x^{2})y''-2xy'+6y=0$

(Eq 5-3)

Using reduction of order, some function, v, times P₂(x) and the first two derivatives of the product are substituted into the Legendre differential equation in equation 5-7. Furthermore, since P₂(x) is a solution to the Legendre equation the coefficient of the first order term, v, reduces to zero. Since the lowest order term cancels out, the second order differential equation is reduced to a first order differential equation where w=v' and w'=v in Equation 5-8 which can be solved through separation of variables.

$R_{2}(x)=v({\dfrac {3}{2}}x^{2}-{\dfrac {1}{2}})$

(Eq 5-4)

$R'_{2}(x)=v'({\dfrac {3}{2}}x^{2}-{\dfrac {1}{2}})+v(3x)$

(Eq 5-5)

$R''_{2}(x)=v''({\dfrac {3}{2}}x^{2}-{\dfrac {1}{2}})+6v'x+3v$

(Eq 5-6)

$(1-x^{2})({\dfrac {3}{2}}x^{2}-{\dfrac {1}{2}})v''+(6x(1-x^{2})-2x({\dfrac {3}{2}}x^{2}-{\dfrac {1}{2}}))v'+(3(1-x^{2})-2x(3x)+6({\dfrac {3}{2}}x^{2}-{\dfrac {1}{2}}))v=0$

(Eq 5-7)

${\dfrac {dw}{w}}={\dfrac {(6x(1-x^{2})-2x({\dfrac {3}{2}}x^{2}-{\dfrac {1}{2}}))dx}{(1-x^{2})({\dfrac {3}{2}}x^{2}-{\dfrac {1}{2}})}}$

(Eq 5-7)

Through partial fraction decomposition, the integral is

$\int {\dfrac {dw}{w}}=\int {\dfrac {-6x}{{\dfrac {3}{2}}x^{2}dx-{\dfrac {1}{2}}}}+\int {\dfrac {2x}{1-x^{2}}}dx$

(Eq 5-8)

Through integration

$log(w)=-2log({\dfrac {3}{2}}x^{2}-{\dfrac {1}{2}})-log(1-x^{2})$

(Eq 5-9)

Using properties of the natural logarithm and then inverting it using the exponential function, equation 5-9 becomes

$w={\dfrac {1}{({\dfrac {3}{2}}x^{2}-{\dfrac {1}{2}})^{2}(1-x^{2})}}$

(Eq 5-10)

Substituting v' back into Equation 5-10 and again using partial fraction decomposition

$\int {dv}={\frac {1}{2}}\int {\dfrac {dx}{(1-x)}}+{\frac {1}{2}}\int {\dfrac {dx}{(1+x)}}+3\int {\dfrac {(3x^{2}+1)}{(3x^{2}-1)^{2}}}dx$

(Eq 5-11)

Using integration by parts

$\int {udv}=uv-\int {vdu}$

(Eq 5-12)

where u and v are defined as follows:

$dv=x(3x^{2}-1)^{-2}$

(Eq 5-13)

$u={\frac {(3x^{2}+1)}{x}}$

(Eq 5-14)

$v={\frac {-1}{2}}log(1-x)+{\frac {1}{2}}log(1+x)-3({\dfrac {(3x^{2}+1)}{6x(3x^{2}-1)}}-\int {\frac {dx}{6x^{2}}})$

(Eq 5-15)

Using the properties of the natural logarithm to combine the first two terms and performing the simple integration on the last term:

$v={\dfrac {1}{2}}log({\dfrac {1+x}{1-x}})-{\dfrac {3x}{3x^{2}-1}}$

(Eq 5-16)

Substituting v into Equation 5-4 then identifying the new solution as Q₂

$Q_{2}(x)={\dfrac {1}{4}}(3x^{2}-1)log({\dfrac {1+x}{1-x}})-{\dfrac {3x}{2}}$

(Eq 5-17)

R2.6: Testing Gauss Laguerre quadrature

The Bessel function integration test performed by Dr. Burkardt's code test the following integral:

$e^{-2}\int _{2}^{\infty }{\frac {\sin {x-1}}{\sqrt {x(x-2)}}}dx$

and has the following raw output:

        7                  0.162669
                1         0.19313       0.0304616
                2       0.0346675        0.128001
                4       0.0367188         0.12595
                8       0.0395037        0.123165
               16       0.0970831       0.0655858
               32        0.100708       0.0619605
               64        0.107105       0.0555637

The top row cow is the problem number, and the exact solution. The first column of numbers is the number of integration points, the second column is the approximated answer, and the third column is the absolute error.

This code did not test a Bessel integral of the form given in problem R2.7.

R3.6: Testing Gauss-Laguerre Quadrature Part 2

Problem Statement

Test the following integrals using Gauss-Laguerre quadrature.

$\int _{e^{\epsilon }}^{\infty }{\frac {1}{x^{2}\log {x}}}dx=\int _{\epsilon }^{\infty }e^{-u}{\frac {1}{u}}du=-Ei(-\epsilon )$

(Eq 6-1)

$\int _{e^{\epsilon }}^{\infty }{\frac {1}{x^{2}\log ^{2}{x}}}dx=\int _{\epsilon }^{\infty }e^{-u}{\frac {1}{u^{2}}}du=Ei(-\epsilon )+{\frac {e^{-\epsilon }}{\epsilon }}$

(Eq 6-2)

Where $\epsilon =0$

Find

Plot the convergence of the integrals in:

(a) (Eq 6-1)

(b) (Eq 6-2)

as the number of integration points are increased.

Solution

On our honor, we did this assignment on our own, without looking at the solutions in previous semesters or other online solutions.

Neither integral given converges, nor has a principle value for $\epsilon =0$. However, the exercise is carried out for a variety of decreasing values of $\epsilon$ ranging from $\epsilon =1$ to $\epsilon =10^{-6}$

The Gauss-Laguerre quadrature rule was given in Equation 7.7 Report2

$\int _{a}^{+\infty }e^{-y}g(y)\,dy=\int _{0}^{+\infty }e^{-(x+a)}g(x+a)\,dx\approx e^{-a}\sum _{i=1}^{n}w_{i}\,g(x_{i}+a)$

(Eq 7-7)

Plots were generated using the following numpy script.

Source:

################################################################################
# Name: gausslagtestEi.py
# Author: Cameron Stewart
# Date: 2/18/2014
#################################################################################
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import expi

# Test functions 1 and 2 are the functions tested by the G-L quad
def TestFunction1(x):
    y=1./x
    return y

def TestFunction2(x):
    y=1./(x**2)
    return y

#################################################################################
# Uses the numpy Laguerre Gauss quadrature routine:
# numpy.polynomial.laguerre.laggauss(deg)  to find roots and weight points.
# Output is tuple of 1d arrays for the weights and integration points.  
# Note: Not tested for deg>100
#################################################################################

weight=(np.polynomial.laguerre.laggauss(4))

# Array for the lower integral bound
a=np.array([1.0,0.1,0.01,0.001,0.0001,0.00001,0.000001])

#################################################################################
# Plot to show the convergence with increasing integration points
#################################################################################
# Function1
exact=-expi(-a)
for j in range(0,len(a)):
    t=np.linspace(2,100,99)
    error=[]
    for i in range(2, 101):
        weight=(np.polynomial.laguerre.laggauss(i))
        func=TestFunction1(weight[0]+a[j])
        integ=weight[1]*func
        integ=np.exp(-a[j])*np.sum(integ)
        error.append((np.abs(integ-exact[j]))/exact[j]*100.)
    plt.plot(t,error,label=r'$\epsilon$ = '+ str(a[j]))
    plt.title(r'Convergence test for Gauss-Laguerre quadrature: $f(x)=\frac{1}{x}$')
    plt.xlabel(r'Integration points: $n$')
    plt.ylabel(r'Error %')
    plt.legend(loc=1)
    plt.ylim(-5,100)
    plt.xlim(2,100)
#plt.show()
plt.clf()

# Function2
exact=expi(-a)+1./(np.exp(a)*a)
for j in range(0,len(a)):
    t=np.linspace(2,100,99)
    error=[]
    for i in range(2, 101):
        weight=(np.polynomial.laguerre.laggauss(i))
        func=TestFunction2(weight[0]+a[j])
        integ=weight[1]*func
        integ=np.exp(-a[j])*np.sum(integ)
        error.append((np.abs(integ-exact[j]))/exact[j]*100.)
    plt.plot(t,error,label=r'$\epsilon$ = '+ str(a[j]))
    plt.title(r'Convergence test for Gauss-Laguerre quadrature: $f(x)=\frac{1}{x^2}$')
    plt.xlabel(r'Integration points: $n$')
    plt.ylabel(r'Error %')
    plt.legend(loc=5)
    plt.ylim(-5,110)
    plt.xlim(2,100)
plt.show()

Part a

Figure 6.1: Percent error of ( 6-1) with increasing integration points for varying values of $\epsilon$

Part b

Figure 6.2: Percent error of ( 6-2) with increasing integration points for varying values of $\epsilon$

In both Figure 6.1 and Figure 6.2, the quadrature converge quickly to the correct value for larger values of $\epsilon$, but fail to converge to the correct answer for smaller $\epsilon$

Contributing Team Members

• Cameron Stewart Solved problem 6, part b of problem 4, and the part of problem 5 pertaining to R2.7. Reviewed problem 3.
• Elizabeth Bartlett
• Kevin Frost
• Christopher Neal

Problem Assignments
Problem #	Solved by	Reviewed by
1	Bartlett, Elizabeth
2	Bartlett, Elizabeth; Neal, Christopher
3	Bartlett, Elizabeth	Stewart, Cameron
4	Frost, Kevin; Stewart, Cameron
5	Bartlett, Elizabeth; Stewart, Cameron
6	Stewart, Cameron