Log out Go to App
Go to App

Learning Materials

Features

Discover

Chapter 5: Problem 18

Rearranging Eq. (5.19) into a slightly more conventional form, we have: $$ \int_{a}^{b} f(x) \mathrm{d} x=h\left[\frac{1}{2} f(a)+\frac{1}{2} f(b)+\sum_{k=1}^{N-1} f(a+k h)\right]+\frac{1}{12} h^{2}\left[f^{\prime}(a)-f^{\prime}(b)\right]+\mathrm{O}\left(h^{4}\right) . $$ This result gives a value for the integral on the left which has an error of order \(h^{4}-a\) factor of \(h^{2}\) better than the error on the trapezoidal rule and as good as Simpson's rule. We can use this formula as a new rule for evaluating integrals, distinct from any of the others we have seen in this chapter. We might call it the "Euler-Maclaurin rule." a) Write a program to calculate the value of the integral \(\int_{0}^{2}\left(x^{4}-2 x+1\right) \mathrm{d} x\) using this formula. (This is the same integral that we studied in Example 5.1, whose true value is \(4.4\).) The order- \(h\) term in the formula is just the ordinary trapezoidal rule; the \(h^{2}\) term involves the derivatives \(f^{\prime}(a)\) and \(f^{\prime}(b)\), which you should evaluate using central differences, centered on \(a\) and \(b\) respectively. Note that the size of the interval you use for calculating the central differences does not have to equal the value of \(h\) used in the trapezoidal rule part of the calculation. An interval of about \(10^{-5}\) gives good values for the central differences. Use your program to evaluate the integral with \(N=10\) slices and compare the accuracy of the result with that obtained from the trapezoidal rule alone with the same number of slices. b) Good though it is, this integration method is not much used in practice. Suggest a reason why not.

Short Answer

Expert verified
Calculate the integral using the Euler-Maclaurin rule, compare with the trapezoidal rule, and discuss why this method is not commonly used.

Step by step solution

01

- Define the function and its derivative

Define the function to be integrated as \( f(x) = x^4 - 2x + 1 \) and its derivative \( f'(x) \). The derivative can be found using the central difference method: \[ f'(a) = \frac{f(a + \epsilon) - f(a - \epsilon)}{2\epsilon} \] \[ f'(b) = \frac{f(b + \epsilon) - f(b - \epsilon)}{2\epsilon} \] Use \( \epsilon = 10^{-5} \) as suggested.
02

- Discretize the interval

Divide the interval \([0,2]\) into \(N=10\) slices. The step size \( h \) is given by: \[ h = \frac{b - a}{N} = \frac{2 - 0}{10} = 0.2 \].
03

- Evaluate the main sum

Calculate the sum of function values at the slice points other than the endpoints, weighted by \( h \): \[ \sum_{k=1}^{N-1} f(a+kh) = f(0.2) + f(0.4) + ... + f(1.8) \].
04

- Compute the Euler-Maclaurin formula

Combine the components to evaluate the integral using the Euler-Maclaurin rule: \[ I = h \left[ \frac{1}{2} f(a) + \frac{1}{2} f(b) + \sum_{k=1}^{N-1} f(a+kh) \right] + \frac{1}{12} h^2 \left[f'(a) - f'(b)\right] + O(h^4) \].
05

- Compute the trapezoidal rule

Evaluate the integral using the trapezoidal rule for comparison: \[ I_{trap} = h \left[ \frac{1}{2} f(a) + \frac{1}{2} f(b) + \sum_{k=1}^{N-1} f(a+kh) \right] \].
06

- Compare results

Compare the results obtained from both the Euler-Maclaurin rule and the trapezoidal rule to the true value of the integral.
07

- Discuss reason for limited use

Discuss why this method is rarely used in practice, despite its accuracy. The reason is primarily its complexity and the availability of simpler and equally effective methods like Simpson's rule.

Unlock Step-by-Step Solutions & Ace Your Exams!

  • Full Textbook Solutions

    Get detailed explanations and key concepts

  • Unlimited Al creation

    Al flashcards, explanations, exams and more...

  • Ads-free access

    To over 500 millions flashcards

  • Money-back guarantee

    We refund you if you fail your exam.

Start your free trial

Over 30 million students worldwide already upgrade their learning with Vaia!

Key Concepts

These are the key concepts you need to understand to accurately answer the question.

Euler-Maclaurin formula
The Euler-Maclaurin formula bridges the gap between discrete sums and continuous integrals. This formula is invaluable in numerical integration because it considerably improves the accuracy of simple methods like the trapezoidal rule. The Euler-Maclaurin formula takes a sum and adds correction terms involving derivatives to approximate an integral more precisely.
For example, in the exercise, the Euler-Maclaurin formula is given by: \[ \int_{a}^{b} f(x) \, \text{d} x = h \left[ \frac{1}{2} f(a) + \frac{1}{2} f(b) + \sum_{k=1}^{N-1} f(a+kh) \right] + \frac{1}{12} h^2 \[ f^{\prime}(a) - f^{\prime}(b) \] + \mathcal{O}(h^4) \]
Here, the first part is simply the trapezoidal rule, but the formula also includes a term with the second derivative \(f^{\prime}(a)-f^{\prime}(b)\). This term significantly reduces the error, making the integration much more accurate. However, in practice, it's not widely used due to its complexity and the need for computing the derivatives accurately.
Trapezoidal rule
The trapezoidal rule is a simple yet powerful method for approximating the definite integral of a function. The idea behind it is straightforward: you approximate the region under the curve as a series of trapezoids rather than rectangles.
Mathematically, for a function \(f\) over an interval \[a, b\], the trapezoidal rule is expressed as: \[ \int_{a}^{b} f(x) \, dx \approx \frac{h}{2} [ f(a) + 2f(a+h) + 2f(a+2h) + ... + 2f(b-h) + f(b) ] \]
Here, \(h \) is the width of each slice. In the exercise, it breaks down into: \[ h = \frac{b - a}{N} \] \[ I_{trap} = h \left[ \frac{1}{2} f(a) + \frac{1}{2} f(b) + \sum_{k=1}^{N-1} f(a+kh) \right] \]
Although the trapezoidal rule is straightforward and easy to apply, its accuracy is limited. This is why methods like the Euler-Maclaurin formula and Simpson's rule, which provide higher accuracy at similar computational costs, are often preferred.
Central difference method
The central difference method is essential when calculating derivatives, especially in numerical settings. It gives an approximation for the derivative of a function, and it's particularly useful in methods like the Euler-Maclaurin formula.
The central difference formula for the first derivative of a function \(f\) is given by: \[ f^{\prime}(x) \approx \frac{f(x + \epsilon) - f(x - \epsilon)}{2\epsilon} \]
In this formula, \(\epsilon\) is a small number used to calculate the difference. In the exercise, \(\epsilon = 10^{-5}\) is used to ensure a good approximation. This method is simple yet effective as it reduces truncation errors, making it suitable for high-accuracy numerical differentiation. By using \`centered differences\` around points \(a\) and \(b\), we can accurately compute the necessary derivatives for the Euler-Maclaurin formula.

One App. One Place for Learning.

All the tools & learning materials you need for study success - in one app.

Get started for free

Most popular questions from this chapter

5.17 The gamma function: A commonly occurring function in physics calculations is the gamma function \(\Gamma(a)\), which is defined by the integral $$ \Gamma(a)=\int_{0}^{\infty} x^{n-1} \mathrm{e}^{-x} \mathrm{~d} x . $$ There is no closed-form expression for the gamma function, but one can calculate its value for given \(a\) by performing the integral above numerically. You have to be careful how you do it, however, if you wish to get an accurate answer. a) Write a program to make a graph of the value of the integrand \(x^{n-1} \mathrm{e}^{-x}\) as a function of \(x\) from \(x=0\) to \(x=5\), with three separate curves for \(a=2,3\), and 4 , all on the same axes. You should find that the integrand starts at zero, rises to a maximum, and then decays again for each curve. b) Show analytically that the maximum falls at \(x=a-1\). c) Most of the area under the integrand falls near the maximum, so to get an accurate value of the gamma function we need to do a good job of this part of the integral. We can change the integral from 0 to \(\infty\) to one over a finite range from 0 to 1 using the change of variables in Eq. (5.67), but this tends to squash the peak towards the edge of the \([0,1]\) range and does a poor job of evaluating the integral accurately. We can do a better job by making a different change of variables that puts the peak in the middle of the integration range, around \(\frac{1}{2}\). We will use the change of variables given in Eq. (5.69), which we repeat here for convenience: $$ z=\frac{x}{c+x} . $$ For what value of \(x\) does this change of variables give \(z=\frac{1}{2}\) ? Hence what is the appropriate choice of the parameter \(c\) that puts the peak of the integrand for the gamma function at \(z=\frac{1}{2}\) ? d) Before we can calculate the gamma function, there is another detail we need to attend to. The integrand \(x^{n-1} \mathrm{e}^{-x}\) can be difficult to evaluate because the factor \(x^{2-1}\) can become very large and the factor \(\mathrm{e}^{-x}\) very small, causing numerical overflow or underflow, or both, for some values of \(x\). Write \(x^{n-1}=\mathrm{e}^{(a-1) \ln x}\) to derive an alternative expression for the integrand that does not suffer from these problems (or at least not so much). Explain why your new expression is better than the old one. e) Now, using the change of variables above and the value of \(c\) you have chosen, write a user-defined function gamma (a) to calculate the gamma function for arbitrary argument \(a\). Use whatever integration method you feel is appropriate. Test your function by using it to calculate and print the value of \(\Gamma\left(\frac{3}{2}\right)\), which is known to be equal to \(\frac{1}{2} \sqrt{\pi} \simeq 0.886\). f) For integer values of \(a\) it can be shown that \(\Gamma(a)\) is equal to the factorial of \(a-\) 1. Use your Python function to calculate \(\Gamma(3), \Gamma(6)\), and \(\Gamma(10)\). You should get answers closely equal to \(2 !=2,5 !=120\), and \(9 !=362880\).

Quantum uncertainty in the harmonic oscillator In units where all the constants are 1 , the wavefunction of the \(n\)th energy level of the one-dimensional quantum harmonic oscillator-i.e., a spinless point particle in a quadratic potential well-is given by $$ \psi_{n}(x)=\frac{1}{\sqrt{2^{11} n ! \sqrt{\pi}}} \mathrm{e}^{-x^{2} / 2} H_{n}(x), $$ for \(n=0 \ldots \infty\), where \(H_{n}(x)\) is the \(n\)th Hermite polynomial. Hermite polynomials satisfy a relation somewhat similar to that for the Fibonacci numbers, although more complex: $$ H_{n+1}(x)=2 x H_{n}(x)-2 n H_{n-1}(x) . $$ The first two Hermite polynomials are \(H_{0}(x)=1\) and \(H_{1}(x)=2 x\). a) Write a user-defined function \(\mathrm{H}(\mathrm{n}, \mathrm{x})\) that calculates \(H_{n}(x)\) for given \(x\) and any integer \(n \geq 0\). Use your function to make a plot that shows the harmonic oscillator wavefunctions for \(n=0,1,2\), and 3, all on the same graph, in the range \(x=-4\) to \(x=4\). Hint: There is a function factorial in the ath package that calculates the factorial of an integer. b) Make a separate plot of the wavefunction for \(n=30\) from \(x=-10\) to \(x=10\). Hint If your program takes too long to run in this case, then you're doing the calculation wrong-the program should take only a second or so to run. c) The quantum uncertainty in the position of a particle in the \(n\)th level of a harmonic oscillator can be quantified by its root-mean-square position \(\sqrt{\left\langle x^{2}\right\rangle}\), where $$ \left\langle x^{2}\right\rangle=\int_{-\infty}^{\infty} x^{2}\left|\psi_{n}(x)\right|^{2} \mathrm{~d} x . $$ Write a program that evaluates this integral using Gaussian quadrature on 100 points, then calculates the uncertainty (i.e., the root-mean-square position of the particle) for a given value of \(n\). Use your program to calculate the uncertainty for \(n=5\). You should get an answer in the vicinity of \(\sqrt{\left\langle x^{2}\right\rangle}=2.3 .\)

Period of an anharmonic oscillator The simple harmonic oscillator crops up in many places. Its behavior can be studied readily using analytic methods and it has the important property that its period of oscillation is a constant, independent of its amplitude, making it useful, for instance, for keeping time in watches and clocks. Frequently in physics, however, we also come across anharmonic oscillators, whose period varies with amplitude and whose behavior cannot usually be calculated analytically. A general classical oscillator can be thought of as a particle in a concave potential well. When disturbed, the particle will rock back and forth in the well: The harmonic oscillator corresponds to a quadratic potential \(V(x) \propto x^{2}\). Any other form gives an anharmonic oscillator. (Thus there are many different kinds of anharmonic oscillator, depending on the exact form of the potential.) One way to calculate the motion of an oscillator is to write down the equation for the conservation of energy in the system. If the particle has mass \(m\) and position \(x\), then the total energy is equal to the sum of the kinetic and potential energies thus: $$ E=\frac{1}{2} m\left(\frac{\mathrm{d} x}{\mathrm{~d} t}\right)^{2}+V(x) $$ Since the energy must be constant over time, this equation is effectively a (nonlinear) differential equation linking \(x\) and \(t\). Let us assume that the potential \(V(x)\) is symmetric about \(x=0\) and let us set our anharmonic oscillator going with amplitude \(a\). That is, at \(t=0\) we release it from rest at position \(x=a\) and it swings back towards the origin. Then at \(t=0\) we have \(\mathrm{d} x / \mathrm{d} t=0\) and the equation above reads \(E=V(a)\), which gives us the total energy of the particle in terms of the amplitude. a) When the particle reaches the origin for the first time, it has gone through one quarter of a period of the oscillator. By rearranging the equation above for \(\mathrm{d} x / \mathrm{d} t\) and then integrating with respect to \(t\) from 0 to \(\frac{1}{4} T\), show that the period \(T\) is given by $$ T=\sqrt{8 m} \int_{0}^{a} \frac{\mathrm{d} x}{\sqrt{V(a)-V(x)}} $$ b) Suppose the potential is \(V(x)=x^{4}\) and the mass of the particle is \(m=1 .\) Write a Python function that calculates the period of the oscillator for given amplitude a using Gaussian quadrature with \(N=20\) points, then use your function to make a graph of the period for amplitudes ranging from \(a=0\) to \(a=2\). c) You should find that the oscillator gets faster as the amplitude increases, even though the particle has further to travel for larger amplitude. And you should find that the period diverges as the amplitude goes to zero. How do you explain these results?

Diffraction gratings: Light with wavelength \(\lambda\) is incident on a diffraction grating of total width \(w\), gets diffracted, is focused with a lens of focal length \(f\), and falls on a screen: Theory tells us that the intensity of the diffraction pattern on the screen, a distance \(x\) from the central axis of the system, is given by $$ I(x)=\left|\int_{-w / 2}^{w / 2} \sqrt{q(u)} \mathrm{e}^{i 2 \pi x u / \lambda f} \mathrm{~d} u\right|^{2} $$ where \(q(u)\) is the intensity transmission function of the diffraction grating at a distance \(u\) from the central axis, i.e., the fraction of the incident light that the grating lets through. a) Consider a grating with transmission function \(q(u)=\sin ^{2} \alpha u\). What is the separation of the "slits" in this grating, expressed in terms of \(\alpha\) ? b) Write a Python function \(q(u)\) that returns the transmission function \(q(u)=\sin ^{2} \alpha u\) as above at position \(u\) for a grating whose slits have separation \(20 \mu \mathrm{m}\). c) Use your function in a program to calculate and graph the intensity of the diffraction pattern produced by such a grating having ten slits in total, if the incident light has wavelength \(\lambda=500 \mathrm{~nm}\). Assume the lens has a focal length of 1 meter and the screen is \(10 \mathrm{~cm}\) wide. You can use whatever method you think appropriate for doing the integral. Once you've made your choice you'll also need to decide the number of sample points you'll use. What criteria play into this decision? Notice that the integrand in the equation for \(I(x)\) is complex, so you will have to use complex variables in your program. As mentioned in Section 2.2.5, there is a version of the math package for use with complex variables called cmath. In particular you may find the exp function from cmath useful because it can calculate the exponentials of complex arguments. d) Create a visualization of how the diffraction pattern would look on the screen using a density plot (see Section 3.3). Your plot should look something like this: e) Modify your program further to make pictures of the diffraction patterns produced by gratings with the following profiles: i) A transmission profile that obeys \(q(u)=\sin ^{2} \alpha d \sin ^{2} \beta u\), with \(\alpha\) as before and the same total grating width \(w\), and \(\beta=\frac{1}{2} \alpha\). ii) Two "square" slits, meaning slits with \(100 \%\) transmission through the slit and \(0 \%\) transmission everywhere else. Calculate the diffraction pattern for non-identical slits, one \(10 \mu \mathrm{m}\) wide and the other \(20 \mu \mathrm{m}\) wide, with a \(60 \mu \mathrm{m}\) gap between the two.

Heat capacity of a solid Debye's theory of solids gives the heat capacity of a solid at temperature \(T\) to be $$ C_{V}=9 V \rho k_{B}\left(\frac{T}{\theta_{D}}\right)^{3} \int_{0}^{\otimes_{D} / T} \frac{x^{4} \mathrm{e}^{x}}{\left(e^{x}-1\right)^{2}} \mathrm{~d} x $$ where \(V\) is the volume of the solid, \(\rho\) is the number density of atoms, \(k_{B}\) is Boltzmann's constant, and \(\theta_{D}\) is the so-called Debye temperature, a property of solids that depends on their density and speed of sound. a) Write a Python function \(\mathrm{cv}(\mathrm{T})\) that calculates \(C_{V}\) for a given value of the temperature, for a sample consisting of 1000 cubic centimeters of solid aluminum, which has a number density of \(\rho=6.022 \times 10^{28} \mathrm{~m}^{-3}\) and a Debye temperature of \(\theta_{D}=428 \mathrm{~K}\). Use Gaussian quadrature to evaluate the integral, with \(N=50\) sample points. b) Use your function to make a graph of the heat capacity as a function of temperature from \(T=5 \mathrm{~K}\) to \(T=500 \mathrm{~K}\).
See all solutions

Recommended explanations on Physics Textbooks

Translational Dynamics

Read Explanation

Dynamics

Read Explanation

Fields in Physics

Read Explanation

Atoms and Radioactivity

Read Explanation

Physical Quantities And Units

Read Explanation

Space Physics

Read Explanation
View all explanations

What do you think about this solution?

We value your feedback to improve our textbook solutions.

Study anywhere. Anytime. Across all devices.

Sign-up for free