<jupyter_start><jupyter_text>Short Exercises Using the trapezoid rule and Simpson's rule estimate the following integrals with the following number of intervals: $2,4,8,16,\dots, 512$. - $\int_0^{\pi/2} e^{\sin x} \,dx \approx ~3.104379017855555098181$- $\int_0^{2.405} J_0(x) dx  \approx 1.470300035485$, where $J_0(x)$ is a Bessel function of the first kind given by $$ J_\alpha(x) = \sum_{m=0}^\infty \frac{(-1)^m}{m! \, \Gamma(m+\alpha+1)} {\left(\frac{x}{2}\right)}^{2m+\alpha}.  $$   Solution The file $\texttt{ch15.py}$ contains all of the functions contained in the Chapter 15 notes. It will be imported in order to use the $\texttt{trapezoid}$ and $\texttt{simpsons}$ functions. It is also acceptable to paste the functions individually.- We will then define the endpoints $\texttt{a}$ and $\texttt{b}$ for each part, and the function $\texttt{f}$ to be integrated for each part. In addition, the variable $\texttt{intervals}$ is defined to store the intervals we want to solve for.- The functions are then iterated and printed.<jupyter_code>import numpy as np
import math
from ch15 import *

# Define array of intervals
intervals = np.array([2,4,8,16,32,64,128,256,512])

# Define bounds and function for part a
f = lambda x: np.exp(np.sin(x))
a = 0
b = np.pi/2

# Calculate for part a and print
print('Estimating int_0^pi/2 of e^sin(x) dx\n')
print("\t Trapezoid\t\t Simpson's")
i = 0
for pieces in intervals:
    ansTrap = trapezoid(f,a,b,pieces)
    ansSimp = simpsons(f,a,b,pieces)
    print(intervals[i],"\t %.15f" % ansTrap,"\t %.15f" % ansSimp)
    i += 1
    
# Define bounds and function for part b
def f(x, M = 100):
    """Order zero Bessel function of the first-kind
    evaluated at x
    Inputs:
        alpha:  value of alpha
        x:      point to evaluate Bessel function at
        M:      number of terms to include in sum
    Returns:
        J_0(x)
    """
    total = 0.0
    for m in range(M):
        total += (-1)**m/(math.factorial(m)*math.gamma(m+1))*(0.5*x)**(2*m)
    return total 
a = 0
b = 2.405

# Calculate for part b and print
print('\nEstimating int_0^2.405 of J_0(x) dx\n')
print("\t Trapezoid\t\t Simpson's")
i = 0
for pieces in intervals:
    ansTrap = trapezoid(f,a,b,pieces)
    ansSimp = simpsons(f,a,b,pieces)
    print(intervals[i],"\t %.12f" % ansTrap,"\t %.12f" % ansSimp)
    i += 1<jupyter_output><empty_output><jupyter_text>As expected, the answers converge to the exact integral.  Inverse Fourier Transform Consider the neutron diffusion equation in slab geometry an infinite, homogeneous medium given by $$-D \frac{d^2}{dx^2}\phi(x) + \Sigma_\mathrm{a} \phi(x) = \delta(x),$$where $\delta(x)$ is the Dirac delta function. This source is equivalent to a planar source inside the slab at $x=0$.  One way to solve this problem is to use a Fourier transform.  The Fourier transform of a function can be defined by$$\mathcal{F}\{ f(x)\}= \hat{f}(k) = \frac{1}{\sqrt{2\pi}} \int\limits_{-\infty}^\infty dx\, f(x) (\cos kx - i \sin kx).$$The Fourier transform of the  diffusion equation above is $$(D k^2+ \Sigma_\mathrm{a}) \hat{\phi}  = \frac{1}{\sqrt{2\pi}}.$$We can solve this equation for $\hat{\phi}(k)$, and then apply the inverse Fourier transform:$$\mathcal{F}^{-1}\{ \hat{f}(k)\}= {f}(x) = \frac{1}{\sqrt{2\pi}} \int\limits_{-\infty}^\infty dk\, \hat{f}(k) (\cos kx + i \sin kx).$$This leads to the solution being defined by $$\phi(x) =  \int\limits_{-\infty}^\infty  \frac{\cos kx\,dk}{{2\pi} (D k^2 + \Sigma_\mathrm{a})}  + i  \int\limits_{-\infty}^\infty  \frac{\sin kx\,dk}{{2\pi} (D k^2 + \Sigma_\mathrm{a})}  .$$The imaginary part of $\phi$ is  zero because $\phi$ is real. You can see that this is so because the integrand of the imaginary part is odd and the integral is symmetric about 0.Your task is to compute the value of $\phi(x)$ at various points using $D = \Sigma_\mathrm{a} = 1$. Because you cannot integrate to infinity you will be computing integrals of the form$$\int\limits_{-L}^L f(x)\,dx,$$for large values of $L$. - Compute value of $\phi(x)$ at 256 points in $x \in [-3,3]$ using Simpson's and the trapezoidal rule with several different numbers of intervals (pieces) in the integration {\em and} using different endpoints in the integration, $L$. Plot these estimates of $\phi(x)$.- Plot the error between your estimate of $\phi(1)$ and the true solution of $\frac{1}{2} e^{-1}$. Make one graph each for trapezoid and Simpson's rule where the $x$-axis is $h$ and the $y$-axis is the absolute error.  On each plot show a curve for the error decay for $L=10, 1000, 10^5, 10^8.$- Give your best estimate, using numerical integration, for the absorption rate density of neutrons at $x=2.$ The file $\texttt{ch15.py}$ contains all of the functions contained in the Chapter 15 notes. It will be imported in order to use the $\texttt{trapezoid}$ and $\texttt{simpsons}$ functions. It is also acceptable to paste the functions individually.- First, all necessary constants, the domain of $x$ values, the exact function, intervals, and $L$s to integrate at are defined. Three dimensional arrays, $\texttt{phiSimp}$ and $\texttt{phiTrap}$ are defined to store the integrals for each given $L$, interval, and $x$ value.- A set of nested $\texttt{for}$ loops are then defined to solve for each value as desired. Note that the function for $f(x)$ must be re-defined for each $x$.<jupyter_code>import numpy as np
import matplotlib.pyplot as plt
from ch15 import *
%matplotlib inline
# Define given constants
D = 1.0
Sig_a = 1.0

# Define domain
points = 256
xPoints = np.linspace(-3.0,3.0,points)

# Define exact function
exact = lambda x: 0.5*np.exp(-np.absolute(x))
xExact = exact(xPoints)

# Define multiple intervals and endpoints
powers = np.array([2,3,4])
Ls = 10**powers
intervals = 10**powers

# Create empty arrays to store approximations of phi
phiSimp = np.zeros([Ls.size,intervals.size,xPoints.size])
phiTrap = np.zeros([Ls.size,intervals.size,xPoints.size])

# Loop through endpoints
nL = 0
for L in Ls:
    # Loop through intervals
    nInterval = 0
    for interval in intervals:
        nx = 0
        # Loop through values of x
        for x in xPoints:
            # Define function for given value of x
            f = lambda k: np.cos(k*x)/(2*np.pi*(D*k**2 + Sig_a))
            # Solve Simpson's
            phiSimp[nL,nInterval,nx] = simpsons(f,-L,L,interval)
            # Solve trapezodial
            phiTrap[nL,nInterval,nx] = trapezoid(f,-L,L,interval)
            # Increment x index
            nx += 1
        # Increment interval index
        nInterval += 1
    # Increment L index
    nL += 1<jupyter_output><empty_output><jupyter_text>Now that we have all of the numerical solutions, we will plot them.-Seperate plots are made for each value of $L$, and for each method. A $\texttt{for}$ loop is used to loop through each value of $L$, and then two plots are made for each loop (one for Simpson's and one for the trapezodial rule). In addition, the exact solution is plotted on each plot as a comparison.<jupyter_code># Loop through each value of L
nL = 0
for L in Ls:
    # Plot for Simpson's for Ls[i]
    plt.title("Simpson's approximation of $\phi(x)$ for $L = 10^" + str(powers[nL]) + "$")
    plt.xlabel('$x$')
    plt.ylabel('$\phi$')
    nInterval = 0
    # Plot the line for each set of intervals
    for interval in intervals:
        plt.plot(xPoints,phiSimp[nL,nInterval,:],label="$10^" + str(powers[nInterval]) + "$ intervals")
        nInterval += 1
    plt.plot(xPoints,xExact,'--',color="black",label="Exact solution")
    plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
    plt.show()
    
    # Plot for trapezoidal for Ls[i]
    plt.title("Trapezoidal approximation of $\phi(x)$ for $L = 10^" + str(powers[nL]) + "$")
    plt.xlabel('$x$')
    plt.ylabel('$\phi$')
    nInterval = 0
    # Plot the line for each set of intervals
    for interval in intervals:
        plt.plot(xPoints,phiTrap[nL,nInterval,:],label="$10^" + str(powers[nInterval]) + "$ intervals")
        nInterval += 1
    plt.plot(xPoints,xExact,'--',color="black",label="Exact solution")
    plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
    plt.show()
    nL += 1<jupyter_output><empty_output><jupyter_text>The oscillatory motion can be clearly seen when the number of intervals is low. $\br$Of interest is that for the approximation with the smallest value of $L$, the solution with $10^3$ intervals approximations closely matches the exact solution. With larger values of $L$, the solution with $10^3$ intervals does not follow the exact solution. With $10^4$ intervals, the solution for each set of $L$ follows the exact solution well. It would seem that there is a ratio between $L$ and the number of intervals that allows for a better approximation. This will be seen later.$\br$Next, we will evaluate the numerical solution at $x = 1$ with varying values of $L$ and the number of intervals in order to evaluate the errors at that specific point.$\br$The exact solution and function $f(x)$ are then redefined for this evaluation. In addition, new intervals and endpoints are defined. $\br$As done previously, variables are defined to store the approximations. The variables $\texttt{errorSimp}$ and $\texttt{errorTrap}$ are used to store the errors at each point.$\br$A set of two $\texttt{for}$ loops is used to run through each set of $L$ and each set of intervals. The resulting set of errors is then plotted for each $h$.<jupyter_code># Define constant x and exact solution
x = 1
exact = np.exp(-1)/2

# Re-define f with new constant x
f = lambda k: np.cos(k*x)/(2*np.pi*(D*k**2 + Sig_a))

# Define multiple intervals and endpoints
Lpowers = np.array([1,3,5,8])
Ls = 10**Lpowers
intervals = 10**np.array([0,1,2,3,4,5])

# Create empty arrays to store approximations of phi
phiSimp = np.zeros([Ls.size,intervals.size])
phiTrap = np.zeros([Ls.size,intervals.size])
errorSimp = np.zeros([Ls.size,intervals.size])
errorTrap = np.zeros([Ls.size,intervals.size])
h = np.zeros([Ls.size,intervals.size])

# Loop through endpoints
nL = 0
for L in Ls:
    # Loop through intervals
    nInterval = 0
    for interval in intervals:
        # Define h
        h[nL,nInterval] = 2*L/interval
        # Solve Simpson's and determine error
        phiSimp[nL,nInterval] = simpsons(f,-L,L,interval)
        errorSimp[nL,nInterval] = np.fabs(exact-phiSimp[nL,nInterval])
        # Solve trapezodial and determine error
        phiTrap[nL,nInterval] = trapezoid(f,-L,L,interval)
        errorTrap[nL,nInterval] = np.fabs(exact-phiTrap[nL,nInterval])
        # Increment interval index
        nInterval += 1
    # Increment L index
    nL += 1

# Plot for Simpson's
plt.title("Absolute error for Simpson's")
plt.xlabel("$h$")
plt.ylabel("Error")
nL = 0
# Loop through values of L
for L in Ls: 
    plt.loglog(h[nL,:],errorSimp[nL,:],'o-',label="$L = 10^ " + str(Lpowers[nL]) + "$")
    nL += 1
plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
plt.show()

# Plot for trapezodial
plt.title("Absolute error for trapezodial rule")
plt.xlabel("$h$")
plt.ylabel("Error")
nL = 0
# Loop through value sof L
for L in Ls: 
    plt.loglog(h[nL,:],errorTrap[nL,:],'o-',label="$L = 10^ " + str(Lpowers[nL]) + "$")
    nL += 1
plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
plt.show()<jupyter_output><empty_output><jupyter_text>It appears that there is a range where the error is roughly first order for each method, and then for any larger number of intervals it flatlines at a constant error. $\br$Lastly, we will see how close we can determine the exact solution at $x = 2$. Any reasonable attempt at this is acceptable. By looking at the plot of errors and the general plots in part 1, you should notice that a ratio of $L$ to intervals of 0.1 returns a reasonable error with minimum calculation.<jupyter_code># Define constant x and exact solution
x = 2
exact = np.exp(-2)/2

# Re-define f with new constant x
f = lambda k: np.cos(k*x)/(2*np.pi*(D*k**2 + Sig_a))

# Define L and intervals
L = 10**7
intervals = 10**8

# Determine estimation
ans = trapezoid(f,-L,L,intervals)

# Print to user
print('For phi(2):')
print("Numerical solution: %.15f" % ans)
print("Exact solution: %.15f" % exact)<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup