Physics 230 -- Lab 5 (Integration)

Fall 2015 - Transtrum

Labs developed by John S. Colton, BYU Physics & Astronomy, based in part on past labs by Branton Cambell and Steve Turley

This lab is designed to familiarize you with Mathematica's tools for symbolic and numerical integration in the context of moderately-challenging math and physics problems.

Clear["`*" ]

Remember these two often useful shortcuts:
Shift-Enter(keypad) - evaluates cell and advances to the next evaluatable cell
Control-A followed by Enter(keypad) - evaluates all cells in the notebook.

1D Integrals

Basic integration

Most of your basic integration needs will be handled by the Integrate command. The help file provides many examples and options which are worth glancing over, but since you’ve already learned some basic integration in Lab 1, let’s just jump in.

Assignment 1. Two basic integrals

(a) Integrate (Subscript[J, 1])(x) as an indefinite integral. Plot the integral and integrand together, over the range {x, 0, 20}.  As mentioned in the Differentiation lab, this is a Bessel function (of the “first kind”), used by physicists in problems involving waves with circular boundaries. It’s obtained in Mathematica via BesselJ. The integral of (Subscript[J, 1]) is related to (Subscript[J, 0]). 
(b) Compute (∫_(-∞)^( ∞)) (e^(-(x^2)))dx, a definite integral.  This represents the total area under the (e^(-(x^2))) function, so I’d like you to also plot that function using the Filling -> Axis option to visualize it. Have your plot go from -5 to 5.

The length of an arc

When you plot a curve as y(x), we often use integrals to calculate the “area under the curve”. However, you can also use integration to calculate the length of the curve. In other words, if you were an ant following the curve on a piece of paper, how long is the path that you walked? 

To do that calculation, note that from the Pythagorean Theorem the length of an infinitesimal section of a curve is (Sqrt[(dx^2)+(dy^2)]) ,  which equals  (Sqrt[1+((((dy/dx)))^2)]) dx, so the total length of the curve must the integral of that: path length=∫(Sqrt[1+((((dy/dx)))^2)])dx, integrated over the appropiate range in x.

Assignment 2. How long is a parabola?

(a) Consider the section of parabola defined by the equation y = 1-(x^2) over the range -1 < x < 1. Make a plot of y(x), and give a rough estimate (but as accurate as you can) of the total length of that curve.
(b) Find the derivative, (dy/dx), and use it in the path length formula to determine the precise value of the path length (in decimal format).

Alternate Input Format

By this time you may have realized that the “nice formatting” that you can obtain with e.g. the Basic Math Assistant palette can also be used for entering in commands. I really dislike doing things that way myself, but some students like it. You can use it if you wish. See how these two commands are equivalent.

Integrate[x^2 Sin[x], x]

∫(x^2)Sin[x]\[DifferentialD]x

Assumptions

Assumptions in integrals

Periodically you run into integrals whose results depend on details of the parameters involved in the expression you are integrating. For example, typically when you are integrating a power of x, say (x^n), you raise the power by 1 and divide by the new power: ∫(x^n)dx =(1/(n+1))(x^(n+1)).  However, if n = -1, that rule doesn’t hold (otherwise you’d be dividing by 0). For that case, the integral is given by ∫(x^(-1)) dx =ln(x). Therefore, if you tell Mathematica to integrate (x^n), without specifying n, it doesn’t quite know what to do. Try it:

Integrate[x^n,{x,0,2}]

The ConditionalExpression result is a sign that Mathematica hasn’t really given you the answer to the integral. It means that the answer, ((2^(1+n))/(1+n)), is only going to be valid when the condition is true (in this case, that the real part of n is > -1).

The way to handle situations like that is through specifying what Assumptions Mathematica should make while calculating the integral. There are at least two nearly equivalent ways to do that: one using the Assuming command, and one using the Assumptions specifier within the Integrate command itself. I personally prefer the first method.

Assuming[n==-1, Integrate[x^n,{x,0,2}]]  (* Note the double equals. Also, integral doesn't converge because you'd have to take ln(0) *)
Assuming[n>-1, Integrate[x^n,{x,0,2}]]  (* works *)

Integrate[x^n,{x,0,2},  Assumptions->n == -1] 
Integrate[x^n,{x,0,2},  Assumptions->n > -1]

Doing it the first way you just “wrap” the Assuming command around the Integrate command, specifying the assumptions as the first argument of the Assuming command and the Integrate command as the second argument. The second method requires you to provide the assumptions directly as part of the Integrate command. In both cases you can specify multiple assumptions at once via a list.

Warning:  Do not trust Mathematica blindly! This is good advice for dealing with any type of computational computer program--including any programs you may write yourself, for example later this semester or in Physics 330 and 430. The program is only as smart as the people who have programmed it, and there’s no guarantee that they have always considered things completely. In this particular case, I discovered that Mathematica does not always tell you when it's making assumptions on its own. Consider this:

Integrate[x^n,x]

That is clearly not the right answer for n = -1! So that command should really return a conditional expression--but it doesn't.
The moral: Do not blindly trust everything Mathematica tells you! Mathematica is not infallible!

Assumptions in other commands

Assumptions can be made via other commands as well. As another example, consider this integral:  (∫_0^∞) (e^sx)dx. The function (e^sx) will be a rising exponential if s is positive, in which case there’s an infinite area under the curve so the integral doesn’t converge; the function will be a decaying exponential if s is negative, in which case there’s a finite area and the integral does converge.

myIntegral=Integrate[E^(s x), {x,0,Infinity}]

Some commands, such as Simplify, allow you to simplify expressions--including conditional expressions such as that--by including assumption-like parameters inside the command. In the case of Simplify, you don’t need to type out the word Assumptions or Assuming, you just need to add the information itself as an additional argument.

Simplify[myIntegral,s<0]

Types of assumptions

There are a wide range of assumptions that can be made--basically anything that Mathematica can use as a test for equality.

Simplify[myIntegral, s==-2]
Simplify[myIntegral,s∈Reals]  (* you can get the "is an element of" symbol, ∈, via <esc>elem<esc> *)
Simplify[myIntegral, {s∈Reals, s<0}]  (* You can provide multiple assumptions via a list *)

Tip: My experience is that typing only “s<0” by itself is enough to force both assumptions in that last command (i.e. that s is a real number and also <0). I generally try the simplest possible assumption that might work, then add additional assumptions if that isn’t enough.

Modifying the global assumptions

Two ways to specify assumptions were given above (three ways, if you count the Simplify command). There's another way, which is to edit Mathematica's main assumptions setting, called the “global” assumptions. This could be useful if you  want to make assumptions about a variable that will apply to all subsequent calculations. But it’s dangerous, because then the assumptions don’t automatically reset when you’re done with the integral. The global assumptions are set via $Assumptions. 

(There are many other similar global parameters that that you may run across, which are also indicated with the character “$”  at the front. They are sometimes called “system variables”. Many of those system variables can be modified this same way.)

$Assumptions = {s<0};
 Integrate[Exp[s x],{x,0,Infinity}]

Do that only with great care, though. If you now try to specify a different assumption for s, you will get a “Contradictory assumptions” error message:

Simplify[myIntegral, s==0]

To view the global assumptions variable, you can enter this :

$Assumptions

To reset the global assumptions variable, you can set it equal to an empty list like this:

$Assumptions = {}

Quitting the kernel also resets all global variables.

Application: Orthogonal functions

Two functions are considered to be “orthogonal” if their product integrates to zero over a particular range. That’s by analogy to vectors, where two vectors are called orthogonal if their dot product equals zero. Orthogonal functions arise in a wide variety of physics contexts, including electromagnetic waves, heat flow, quantum mechanics, and so forth.

As an example, sin(x) and sin(2x) are orthogonal over the range of 0 to 2π.

f1=Sin[x]
f2=Sin[2x]
Integrate[f1*f2,{x,0,2 Pi}]

You can see graphically why this is the case: when you multiply them together there is as much area above the x-axis as there is below it.

Plot[f1 f2,{x,0,2 Pi},Filling->Axis]

By contrast, if you integrate sin(2x)*sin(2x), all the area in the product is above the x-axis and so the integral is non-zero. So sin(2x) is not orthogonal to itself.

Plot[f2 f2,{x,0,2 Pi},Filling->Axis]
Integrate[f2 f2,{x,0,2Pi}]

Assignment 3. Orthogonal functions: sines and cosines

Note: Using Assumptions inside an Integrate command is apparently not quite identical to using the Assuming command. When using the former, I sometimes have to additionally use the Refine command to get it fully simplified, as in “Refine[%, {n∈Integers, m ∈Integers}]”. Unfortunately I haven’t yet been able to figure out a pattern for when that will/will not be neccesary. It was necessary a couple of times for me when I did the following assignment.

(a1) Show that all sin(nx) and sin(mx) functions are orthogonal to each other over the integral 0 to 2π, when n does not equal m. Use Assuming or Assumptions to specify that n ≠ m, and that both n and m are integers. 
(a2) Find the value of the integral of sin(nx) sin(mx) when n and m are equal. (Again, use an assumption to specify that n==m. Note the double equals sign there.) 
(b) Repeat part (a), for cos(nx) and cos(mx).
(c) Repeat, for sin(nx) and cos(mx).

Summary: You should have found that sine and cosine functions are always orthogonal, that sine functions are orthogonal to each other when n and m are not equal, and that cosine functions are similarly orthogonal when n ≠ m.

Final words: Orthogonal functions are incredibly important in physics applications and will be a major part of your study in Phys 318, 441, and 452, just to name a few courses. Sines/cosines are just one set of orthogonal functions. Two other important sets include the Bessel functions and the Legendre polynomials (see below for an extra credit problem involving the Legendre polynomials).

Multivariable Integrals

In physics classes we frequently perform integrals of several variables, such as:  ∫∫((x^2)+(y^2))dx dy. If you are unfamiliar with that notation, it means “Integrate (x^2)+(y^2) with respect to x, then integrate the result with respect to y.”  Actually, the order generally does not matter except for certain strange cases you’re not likely to run into in physics classes. You can therefore do the integral with respect to y first, if you want.

Just as the integral of a function of x gives the area under the curve defined by the function, the integral of a function of x and y gives the volume under the surface defined by the function. Mathematica can do that in one step, as such:

Clear[x,y]   (* in case you've set x or y equal to anything *)
Integrate[x^2+y^2,x,y]

If the integrand has a constant value of 1, then the volume under the surface is the same as the area enclosed by the region of integration, ignoring units. Hopefully you can visualize why that is the case--ask a TA if you are not sure.

Similarly, a function of three variables can be defined and integrated: ∫∫∫f(x,y,z)dx dy dz. In this case f(x,y,z) often has the physical meaning of a density such as mass per unit volume or charge per unit volume. The integral is then the total mass or charge enclosed by the volume. Or, if the function f  just has a constant value of 1, then the integral is the volume of the region of space.

Assignment 4. Perform three multivariate integrals

Compute the following integrals.
(a) Compute the indefinite integral ∫∫sin(x y)dx dy.  Plot3D the integrand and the integral on the same graph over the range {x, -π, π} and  {y, -π, π}.  Use the PlotStyle option to give the two surfaces different colors.  (The CosIntegral function that should appear in your result is another of Mathematica’s many built-in functions. You can see what shape it takes by viewing the help file.)
(b) Compute the definite integral ∫∫((x^2)+(y^2))dx dy, with x and y both ranging from -2 to 2. Plot3D the integrand over the same range and use the PlotRange->All and Filling->Axis options to visualize the region of space whose volume you have just calculated.
(c) In spherical coordinates, dx dy dz translates to (r^2)sinθ dϕ dθ dr.  Here r is the distance from the origin to an arbitrary point (x,y,z).  θ is the angle the (x,y,z) vector makes with respect to the z-axis, and ϕ is the angle the projection of the (x,y,z) vector onto the x-y plane makes with respect to the x-axis. (Careful: this is the usual definition for physics classes but math classes often define θ and ϕ the opposite way.) Thus to cover all points inside a sphere of radius R: r ranges from 0 to R, θ ranges from 0 to π, and ϕ ranges from 0 to 2π. Let the function you are integrating be equal to 1 everywhere, and compute the integral ∫∫∫dx dy dz expressed in spherical coordinates to obtain the well-known expression for the volume of a sphere.

Regional integrals

As mentioned above, areas (volumes) can be calculated via 2D (3D) integrals where the function being integrated is just equal to a constant 1. Here’s a trick to be able to define interesting areas. It involves the RegionPlot command to specify the area over which an inequality or set of inequalities are true. For example, the inequality  4<2(x^2)+(y^2)<9   (which is really two inequalities, 4<2(x^2)+(y^2) and 2(x^2)+(y^2)<9) specifies the regions between two ellipses because the first half is only true for the region of space outside a small ellipse and the second half is only true for the region of space outside a larger ellipse. RegionPlot illustrates this, as such:

inequality =   4<2(x^2)+(y^2)<9
RegionPlot[inequality,{x,-4,4},{y,-4,4}]
RegionPlot[(4 < 2 x^2 + y^2) &&  ( 2 x^2 + y^2 < 9),{x,-4,4},{y,-4,4}]  (* An alternate, equivalent way of specifying things. The symbol && means both the first inequality AND the second one must be true. *)

Now here’s the trick: the Boole function takes an inequality as input, and returns 1 at points where inequality is satisfied or 0 otherwise.  Such a function is best visualized with a plot.

Plot3D[Boole[inequality],{x,-4,4},{y,-4,4}, Filling->Axis]

Therefore if we do a 2D integral of the Boole function defined that way, we will obtain the volume of the non-zero region which (ignoring units) is exactly the same as the area of the region illustrated in the RegionPlot command previously. This method can be extended to volumes of 3D space or lengths of 1D space. In this case, here’s the answer:

Integrate[Boole[inequality],{x,-4,4},{y,-4,4}] //N

(You can check this against the area of two ellipses, one subtracted from the other, if you wish.)

Going one step farther: It is also possible to integrate a separate function f  over the specific region where inequality is satisfied.  To do this, we simply multiply f by Boole[inequality] and integrate.  In the cell below, we use our inequality to define a region of space, combined with a function, 2(x^2) + (y^2)-4, that increases away from the origin, to define a “football stadium” function of x and y. Then we integrate f over that region to find the volume under the bleachers. (Warning: this does take a bit of time to execute.)

f = 2x^2+y^2-4;
Plot3D[f,{x,-4,4},{y,-4,4}]
stadium = f*Boole[inequality];
Plot3D[stadium,{x,-4,4},{y,-4,4},Filling->Bottom,PlotRange->{0,10}]
Integrate[stadium,{x,-4,4},{y,-4,4}] //N

Assignment 5. The area of a circle, two ways

(a) Determine the area of a circle (radius 1) by using a conventional 1D integral in x. (Hint: use the equation of a circle, (x^2)+ (y^2)=1 to solve for y.) Think carefully about your limits of integration and what the integral represents. A Plot command with Filling->Axis may help.
(b) Determine the area of the same circle by using RegionPlot (and Plot3D if you wish) to visualize the inquality (x^2)+(y^2)≤ 1, then Boole and Integrate to determine the area inside the circle.

Numerical Integration With NIntegrate

We discussed NIntegrate a bit in a previous lab. As was mentioned previously, it’s different than using the Integrate command followed by the N command. The basic idea is that NIntegrate evaluates the integral using strictly numerical techniques, calculating the area under a curve by adding up a bunch of smaller areas using numerical techniques such as the trapezoidal rule. By contrast, Integrate followed by N first does the integral analytically (if possible), and then converts the result to a number. Some more information is given in the Numerical Mathematics in the Wolfram Language tutorial.

Since it may take Mathematica a while to do the integral symbolically, or it may “spin its wheels” before deciding it cannot do the integral symbolically, NIntegrate may well be faster if all you care about is a numerical result. Of course, if the integral is not super hard for Mathematica to do symbolically, Integrate may be the way to go because it will generally still be fast and will give you an exact answer instead of a numerical estimate. 

Here’s an example of an easy integral where there’s essentially no difference between the two commands. I’ll use the Timing command (in postfix form) to display how long the calculations take. The first number that Timing spits out is the time of the calculation, in milliseconds; the second number is the result of the calculation. Both commands are very fast.

Integrate[x^2,{x,0,10}]//N  // Timing
NIntegrate[x^2,{x,0,10}]// Timing

In this next example, however, NIntegrate takes significantly less time.

Integrate[1/Sqrt[x^2 + y^2+z^2],{x,-1,1},{y,-1,1},{z,-1,1}]   //Timing  (* took 25 seconds on my office computer *)
NIntegrate[1/Sqrt[x^2 + y^2+z^2],{x,-1,1},{y,-1,1},{z,-1,1}]  //Timing  (* took 0.1 seconds on my office computer *)

Because Integrate solves the integral exactly, you can always take the result out to more decimal places.

Integrate[1/Sqrt[x^2 + y^2+z^2],{x,-1,1},{y,-1,1},{z,-1,1}]   (* takes a while, but gives an exact answer *)
N[%,50]

You can’t necessarily do that with NIntegrate.

NIntegrate[1/Sqrt[x^2 + y^2+z^2],{x,-1,1},{y,-1,1},{z,-1,1}]  
N[%,50]

The next section talks about how to use NIntegrate’s settings to get accurate results to as many decimal places as you desire.

Basic NIntegrate Options: WorkingPrecision, PrecisionGoal, AccuracyGoal

The tutorials on Numerical Precision and Controlling the Precision of Results discuss how Mathematica treats precision and accuracy in numerical calculations, including the WorkingPrecision, PrecisionGoal, and AccuracyGoal options for NIntegrate. They are short. Read them and look at the examples.

For those who didn’t follow the last instruction, here’s a summary. But you should really read through the tutorials.

Precision	the number of significant digits stored internally
Accuracy	the number of significant digits that are to the right of the decimal place (although it can be fractional, somehow)

WorkingPrecision	the number of internal digits Mathematica should use for computations
PrecisionGoal		the number of digits of precision Mathematica should try to get
AccuracyGoal	the number of digits of accuracy Mathematica should try to get

Therefore the number of decimal places you want to obtain in a numerical integration can be set via the WorkingPrecision, PrecisionGoal, and AccuracyGoal options.

Actually, most of the time you can get away with just setting WorkingPrecision. When you use that inside a command like NIntegrate, Mathematica will also automatically adjust PrecisionGoal and AccuracyGoal to values that seem reasonable to it. Be advised, however, that things don’t always work out the way you want. So while this command, for example, should get a numerical answer to the integral above that’s accurate to about 15 digits, it doesn’t really--as you can see by comparing against the exact answer from the Integrate command above.

NIntegrate[1/Sqrt[x^2 + y^2+z^2],{x,-1,1},{y,-1,1},{z,-1,1},WorkingPrecision->15]  (* numerical integration, in theory accurate to ~15 digits *)
N[-2 (π+Log[64]-12 Log[1+(Sqrt[3])]),15] (* exact answer *)

Intermediate Options: MaxPoints, recursion, and specifying singularities

The main reason why the previous example didn’t give the correct answer to 15 decimal places is because there is a singularity at the origin--when x, y, and z are zero, the function we were integrating blows up to infinity. Functions that have singularities or other sharp features in them are the bane of numerical integration. Quickly scan through the basic tutorial on Numerical Integration to learn how Mathematica handles such situations. Look for examples that show you how to use MaxPoints, MaxRecursion, and MinRecursion achieve greater control over the integration process. MinRecursion forces the algorithm to implement a minimum number of recursions; hopefully you remember about recursion from the plotting lab. MaxRecursion limits the total number of recursions. And MaxPoints limits the number of points where the function is sampled.

Assignment 6. Dealing with a singularity

(a) Make a plot of (e^(-(((x-1))^2))), from -1000 to 1000. It’s not really a singularity (as you can see by making the plot range much smaller), but on that large of a scale it acts like one.

(b) If you try to NIntegrate this function over the range fom -1000 to 1000, it will give you an error message and a wrong result.

NIntegrate[Exp[-(x-1)^2],{x,-1000,1000}]

(I know the result is wrong, because I happen to know the answer should be extremely close to sqrt(π). That’s because when the limits are from -∞ to ∞, it’s essentially the integral you did for assignment 1b.)

The main way to automatically deal with singularities in integrals is to specify a point or points that Mathematica should pay special attention to as it’s going through its numerical algorithm. The Numerical Integration tutorial linked to above shows how to do that. After consulting the tutorial, redo the numerical integral, this time specifying x=1 as one such special point.

(c) Unfortunately, even after you specify the point x=1 as a possible singularity, for this particular problem Mathematica still gives you the wrong answer. And it probably won’t even give you any error messages as a warning this time. Remember, Mathematica is not infallible! For critical calculations you should always think about, and possibly double-check, the answers it gives you. That being said, Mathematica is still the top-of-the-line computational mathematical program. It’s not perfect, but it’s very, very good at what it does. 

To make an even better calculation, try using the MaxPoints, MinRecursion, and/or MaxRecursion options. One of them, or possibly a combination, will give you the correct result.

Advanced options using Mathematica’s “integration strategies”.

The purpose of this section is to make you aware of the extensive tutorial on Advanced Numerical Integration -- a great resource when facing a difficult numerical-integration problem. When confronted with an integral, Mathematica attempts to intelligently choose an algorithm. However, you can force it to use a different algorithm via the Method option. There are a dizzying number of algorithms available to choose from. In the advanced numerical integration tutorial, go to the section on NIntegrate Integration Strategies (here’s a direct link as well) to get a feel for how extensive Mathematica’s capabilities are in this regard.

As an example, consider the function: f(x,y)=(1/(Sqrt[(x^2)+(y^2)])).  I wanted to do a 2-dimensional integral with NIntegrate, with x and y both going from 0 to 1, and with a precision of 30 decimal places. Here’s the first command I used to attempt that. (Execute the command yourself.)

NIntegrate[1/Sqrt[x^2+y^2],{x,0,1},{y,0,1}, WorkingPrecision->30]//Timing

After waiting more than 10 seconds, I eventually gave up and aborted the command with Alt-period. (You can do the same.)

It’s clear that the function has a singularity at (0,0), one of the corners of the integration domain. Since singularities often create problems for numerical integration techniques, I decided to hunt around in the NIntegrate Integration Strategies tutorial to see if I could get better results if I specified a particular method to use. After searching for the word “corner”, I discovered the “Duffy Coordinates” method. Here’s what I found: “Duffy's coordinates is a technique that transforms an integrand over a square, cube, or hypercube with a singular point in one of the corners into an integrand with a singularity over a line, which might be easier to integrate.”  Hey! a singular point at one of the corners is exactly what the issue is in this situation! Maybe it will help! And in fact, it did help. 

That’s an example of what you may need to do at times--if you detect issues with your numerical integration, do a bit of searching to see if any of the advanced strategies deal with the type of problem you are having, then try them out to see if they improve matters.

Assignment 7. Solving a problem with an advanced Method

Find the section on Duffy Coordinates in the integration strategies tutorial, and figure out how to use that method (specifying the problematic corner) to get an answer to the 1/Sqrt[x^2+y^2] integral in about a second or less, still using a WorkingPrecision of 30.

Additional Assignments

Application: Electric potential from a spherical shell charge distribution

Consider a uniformly-charged hollow spherical shell of radius R and total charge Q. The charge is uniformly distributed over the whole surface. Suppose you want to find the electric potential from this shell at an arbitrary point (x,y,z). Due to the spherical symmetry this problem can actually be solved quite quickly through Gauss’s Law, but with a less symmetric charge distribution on the shell you’d need to solve the problem through direct integration of the infinitesimal contribution to V from each tiny amount of charge. So let’s do it that way. Because the electric potential from a point charge is V=((k Q)/distance), where k is 1/(4(Subscript[πϵ, 0])), each tiny bit of charge on the shell dQ makes a small contribution to the electric potential at (x,y,z): dV=((k dQ)/distance), where the denominator is the distance from the tiny charge to the point (x,y,z). If you add up all of the small contributions, i.e. integrate, you should obtain the total V(x,y,z). (Note that the since the charge per area is constant, this ratio is true,  dQ/dA = Q/A,  which means that  dQ=(Q/(4((πR)^2)))  dA .) 

Let’s call the center of the sphere (0,0,0). The position where we want to know the field at is (x,y,z). The hard part about setting up the integral is that we’re integrating over the surface of the sphere, which is NOT (x,y,z)  (because (x,y,z) is the point where we want to know the field). Let’s call the integration coordinates (x’,y’,z’). A picture may help. As we integrate, we’re changing the location of (x’,y’,z’) (but it’s always on the sphere), which means that we’re changing the distance from our dQ to (x,y,z).

Note: (x,y,z) could be inside the sphere. The equations don’t care--I put it outside because it was easier to draw that way.

We can simplify things a bit for this case by noticing that due to symmetry the potential at (x,y,z) will be the same for all points that are that same distance away. Let’s then consider a specific point, which is on the z-axis the same distance r away. Equivalently, you can think of this as just rotating the ball and the point (x,y,z) until the point lies on the z-axis.

The critical distance is marked on the figure: for a given point at (0,0,r), how far is it from an arbitrary bit of charge at (x’,y’,z’)? That can be obtained by the Law of Cosines, which because the sphere has radius R, says: (distance^2) =(R^2)+(r^2)-2Rrcosθ. 

Thus the total potential at (x,y,z) is V(r) = ∫((k dQ)/distance)=∫((k  Q/(4((πR)^2))  dA)/(Sqrt[(R^2)+(r^2)-2R r cosθ])). In spherical coordinates,  dA = (R^2)sinθ dθ dϕ, and the integral goes from 0 to 2π in ϕ and from 0 to π in θ.  (Caution: as mentioned in one of the sections above, math classes often reverse the θ and ϕ. In physics classes, though,  θ is nearly always the angle from the z-axis and ϕ the angle as you go around the z-axis.)  I used a capital R in the dA expression because you are integrating only over the surface of the sphere, where the radius is a constant R.

Assignment 8. Electric potential from a spherical shell charge distribution

(a) Perform the integral to get the potential as a function of r.  Hint: Remember to let Mathematica know that r > 0 and R > 0. If you find the integral taking more than a few seconds, or if it gives you bizarre output, try reversing the order you tell Mathematica to perform the θ and ϕ integrals. While that shouldn’t matter--and didn’t to me--it seemed to, to some students in a previous semester.
(b) Now use variable replacement to evaluate your potential for the simple case where {k -> 1, Q -> 1, R -> 1}.  The resulting expression should depend only on r. Plot it over the range 0 < r < 5, using PlotRange->{0, 1.1} as an option.  Interpret this (perhaps surprising) result to your TA. Incidentally, since gravity is also an inverse square law and since field is essentially the derivative of the potential, your result should prove that if the earth were a hollow sphere, there would be no gravitational field on the inside.

Monte Carlo Integration

The most important feature of a numerical integration algorithm is its choice of sample points. The more sample points you can evaluate the function at, the better the numerical integral will approximate the correct answer. However, more sample points means more time spent doing the calculation, and in some cases (as with singularities discussed above) it can take a LOT of time before the numerical integration result gets close to the correct answer. Therefore, the key to any numerical integration strategy is to cleverly select the sample points. 

One of the built-in Methods is the Monte Carlo method, as you can see in the Details and Options section of the NIntegrate reference page. This type of integration is so important in many research areas that it’s worth exploring in more depth. The basic idea is simple: you randomly select the sample points throughout the integration region; the title “Monte Carlo” is a reference to the randomness of gambling, a popular pastime in Monte Carlo. 

Using the Monte Carlo approach, a one-dimensional integral is computed as:
 (∫_a^b)f(x)dx = (((b-a))/n)∑f((Subscript[x, i])), 
where the sum runs over n evaluations of f at random values of (Subscript[x, i]).  
  
The formula basically applies the rectangle method (a simpler version of the trapezoidal rule) to a set of rectangles that have random widths. The larger the number of random samples, the more accurate the integration will be.

Evaluate the cells below (modified from the Neat Examples section of the NIntegrate reference page) to reveal the points sampled while integrating the volume of a particular 3D shape. You don’t have to understand the details of the commands, just obseve the outputs. Specifically, notice how regular the points are spaced when using  Method->”TrapezoidalRule”  compared to how randomly they are spaced when using  Method->”MonteCarlo”. (Monte Carlo doesn’t actually do a great job here, but it works well in many other applications, in particular large dimensional integrals.)

(* setting things up *)
func[{x_,y_}] = Exp[-2(x^2+y^2)];

revealpoints:= Point[ReplacePart[#,3->func[#[[1;;2]]]+0.01]]&
fplot = Plot3D[func[{x,y}],{x,-2,2},{y,-2,2},BoxRatios->{1,1,1},Mesh->False, AxesLabel->{"x","y","z"},PlotRange->All]
Integrate[func[{x,y}],{x,-2,2},{y,-2,2}]//N  (* the answer *)

(* Trapezoidal method *)
options = {Method->"TrapezoidalRule",PrecisionGoal->2,MaxPoints->1000};
samples = Reap[NIntegrate[func[{x,y}],{x,-2,2},{y,-2,2},Evaluate[options],EvaluationMonitor:>Sow[{x,y,0}]]][[2,1]];
Show[fplot,{PointSize[0.015],revealpoints/@samples}//Graphics3D]
NIntegrate[func[{x,y}],{x,-2,2},{y,-2,2},Evaluate[options]]  (* the answer given by the trapezoidal rule *)

(* Monte Carlo method *)
options = {Method->"MonteCarlo",MaxPoints->1000};
samples = Reap[NIntegrate[func[{x,y}],{x,-2,2},{y,-2,2},Evaluate[options],EvaluationMonitor:>Sow[{x,y,0}]]][[2,1]]//Quiet;
Show[fplot,{PointSize[0.015],revealpoints/@samples}//Graphics3D]
NIntegrate[func[{x,y}],{x,-2,2},{y,-2,2},Evaluate[options]]//Quiet  (* the answer given by the Monte Carlo method *)

Assignment 9. Monte Carlo integration of a circle

The purpose of this assignment is to create your own Mathematica function which performs a Monte Carlo integral, without using NIntegrate.

(a) Create the function circ[x_] := Sqrt[1-x^2] and plot it over the range {x, -1, 1} with AspectRatio -> Automatic and Filling -> Axis.  NIntegrate the function over the same range. It’s half of a circle with radius 1, so the area should obviously be π/2.
(b) Now, create a function called montecirc[n_]. The purpose of this function is to implement the Monte Carlo method using the summation formula in the description above (i.e. not by using NIntegrate). The function argument n indicates the number of random points to use in the summation. Specifically, your function should create a list of n random numbers inside the domain, evaluate circ(x) at those points, sum the function evaluations together with Total, and finally multiply everything by (b-a)/n. Calculate montecirc[n] for n = (10^1), (10^2), (10^3), (10^4), (10^5) and (10^6). For each result calculate the fractional error between your result and the true value of the integral: (result - true)/true. Your fractional error should decrease as you use more and more random points in your function.

Orthogonal functions, revisited

Assignment 10, Extra Credit. Legendre polynomials

Another set of orthogonal functions are the “Legendre polynomials” (see LegendreP). The Legendre polynomials (Subscript[L, m])(x) and the “associated Legendre polynomials” are commonly encountered in quantum mechanics when dealing with atomic wave functions. In this case, the orthogonality integral must be done from -1 to 1 instead of from 0 to 2π as was the case for the trig functions. 

(a) Generate a list of the first six Legendre polynomials (i.e. m runs from 0 to 5) and plot them all on the same graph over the range {x, -1, 1}. They are polynomials which always go to the values 1 or -1 at the endpoints of that range.
(b) Pick two of the polynomials (different ones), and plot their product with the Filling->Axis option, to convince yourself it’s plausible that the total area above the graph equals the total area below the graph. Then do the integral with those two specific polynomials to show that it does actually equal zero.
(c) Repeat, using two of the same polynomial instead of two different ones. In this case, the integral won’t equal zero. 
(d) Now construct a table of values of the integral for all values of m and n. Do this with the Table command, using {n,0,5} and {m,0,5} as options. (You don’t have to use any assumptions.) Your table should tell you at a glance that the integral is zero when n and m are not equal. By looking at the pattern figure out what the integral equals when n=m, in general.

Wrap-up

The following functions were discussed in today's lab, or are very closely related:

Symbolic and some numerical integration: Integrate
How to specify assumptions: Assuming, $Assumptions (global setting), Assumptions (option to various commands such as Integrate and Simplify), Refine
Trick to integrate over strange regions: RegionPlot, Boole
How to analyze your speed: Timing
Numerical integration: NIntegrate;  with options: WorkingPrecision, PrecisionGoal, AccuracyGoal, MaxRecursion, MinRecursion, MaxPoints (no help entries for those two)
Special functions: BesselJ, LegendreP

Potentially useful summary/tutorial pages:
Tutorial: OrthogonalPolynomials
Tutorial: Numerical Mathematics in the Wolfram Language 
Tutorials: Numerical Precision and Controlling the Precision of Results 
Tutorial: Advanced Numerical Integration and NIntegrate Integration Strategies