Physics 230 -- Lab 4 (Differentiation)

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 computing limits, series expansions, and derivatives. We will use several moderately-difficult math and physics problems for context.

Clear["`*" ]

As a reminder, here are two useful shortcuts which have been mentioned before:
Shift-Enter(keypad) - evaluates cell and advances to the next evaluatable cell
Control-A followed by Enter(keypad) - evaluates all cells in the notebook.

Limits

Limit calculations

The Limit function can be used to determine limits for several types of situations, including: (a) Trivial limits, such as limit of (x^2) as x approaches 3. (b) Limits that involve L’Hospital’s rule such as when the numerator and denomimator of an expression both go to zero, example: sin(x)/x as x goes to zero. (c) Limits involving infinity, such as 1/x as x goes to infinity. The command format should make perfect sense to you by now:

Limit[x^2, x-> 3]

Limit[Sin[x]/x, x-> 0]

Limit[1/x, x-> Infinity]

Discontinuities. If a function is discontinuitous, the direction of the limit should be specified via a positive or negative number (positive = “from the right”, negative = “from the left”). Consider the Heaviside theta function, also known as the “step function” or Heaviside step function. It is zero when x is negative, and +1 when x is positive. It has no defined value at x=0.

f[x_] =HeavisideTheta[x];
Plot[f[x],{x,-1,1}, PlotStyle->{Blue, Thick}]
f[0] //N    (* no numerical result! *)

Limit[f[x],x->0, Direction-> 1]

Limit[f[x],x->0, Direction->-1]

Assignment 1. Fibonacci numbers

You are probably familiar with the Fibonacci numbers: 1, 1, 2, 3, 5, 8, 13, 21, ..., where each number in the sequence is the sum of the two previous numbers. In Mathematica you can access the numbers via the Fibonacci command. You may or may not know that the ratio of succesive numbers (such as 8/5, 13/8, and 21/13) approaches a constant as the sequence progresses, but you can see that from this next graph. (I’ve added the Floor function into the plot command so that Fibonacci won’t calculate fractional Fibonacci numbers for the graph... which Mathematica can do but which aren’t relevant to the discussion. You won’t need the Floor function in order to find the limit.) Find the exact value of this constant, then google the number to learn its name.

Plot[Fibonacci[Floor[n+1]]/Fibonacci[Floor[n]],{n,1,10},PlotRange->All,AxesOrigin->{0,0},AxesLabel->"Ratio of successive \n Fibonacci numbers"] (* the \n produces a line break *)

Series Expansions

Taylor series

As you should know, the Taylor series of a function about x=0  is a power series expansion in x, where the coefficient of x^n is (1/n!) times (d/dx)^n of f, evaluated at x=0. Visit the reference page for the Series function, which can automatically perform power-series expansions like that.  Browse through the Basic Examples and the Scope section to get a feel for the types of expansions that are possible. In particular notice how the Normal command can be used in conjunction with Series to strip off the final "+ order(x)" business at the end.

Assignment 2. The four Taylor series all physics majors need to know

(a) I tell students there are four Taylor series that all physics majors need to know. They are sin(x), cos(x), (e^x), and (((1+x))^n). Use Mathematica to give the series to enough terms that you can memorize the patterns, and memorize the first few terms such that the pattern is clear. Review some factorials if needed so that you can recognize them in the denominators. Call the T.A. over when you can demonstrate that you have generated the series and that you can write out the first few terms from memory (in terms of the factorials).

(b) Let’s now use three of those series to prove Euler’s identity: (e^ix)=cos(x) + i sin(x). Specifically, use Mathematica to plug ix in place of x, in your Taylor’s series for (e^x) (say, to 20 terms). Then take the Taylor series of sin, multiply by i, and add it to the Taylor series of cosine. Show that when you do that, the two resulting series are the same.

Assignment 3. Relativistic kinetic energy

In Physics 123 and/or 222, you learned that the true kinetic energy formula is  KE=(Subscript[E, tot])-(Subscript[E, rest]), where (Subscript[E, tot])=((γmc)^2) and (Subscript[E, rest])=(mc^2).  In case you’ve forgotten, γ is the Lorentz factor, γ=1/sqrt(1-(β^2)) , and β is the velocity expressed as a fraction of the speed of light, β=v/c.  Use the Series command to calculate the kinetic energy for β<<1.  Show that when you replace beta with v/c, the first nonzero term in the expansion gives the familiar KE=(1/2)(mv^2) formula; if β is really small, that is all that’s important. What’s the second term in the expansion?

Derivatives

As mentioned in earlier labs, derivatives can be calculated either via tick marks like  '  or  '', or via the D function. We’ll now investigate them in more detail. 

At this point, let me warn you again about the difference between between functions that have been “trained” to take arguments and those that have not. By functions that have been trained to take arguments, I mean normal function definitions such as y[x_]=x^2 or y[x_]:=x^2). These are sometimes called “pure functions” or “argumentative functions”. By those that have not, I mean commands such as y=x^2. Recall that entering in a command like y=x^2 just sets the variable y equal to the algebraic expression x^2, and does not really define y as a function of x.  Those two different ways of dealing with functions can often be used interchangeably, but not always. Always think carefully about the results of your derivative commands, and if a command doesn’t give you a result that makes sense, think about using the other type of function definition. Argumentative functions are nearly always more reliable in this context, but longer to type in and therefore more susceptible to typos.

Basic derivatives

Look at the help file on D. In particular, notice how to do the following:
1. Take a first derivative with respect to x.
2. Take a higher deriviative with respect to x.
3. Take a first derivative with respect to x, followed by a first derivative with respect to y.
4. Take a higher derivative with respect to x, followed by a higher derivative with respect to y.

All of those can be done with a single D command.

Assignment 4. Calculating derivatives

Do the following exercises, using a single D command in each case:
(a) Find the numerical value of the second derivative of f(x) = (Sqrt[1+(x^2)]), evaluated at x = 4.

(b) Plot (Subscript[J, 0])(x) along with its first derivative, from 0 to 20. Plot the (Subscript[J, 0]) function in blue and the derivative in red. Note that (Subscript[J, 0])(x) is a special function called a “zeroth-order Bessel function of the first kind” (see BesselJ). The subscript 0 indicates the order; the J indicates that it’s a Bessel function of the “first kind”. These kind of Bessel functions often arise in wave situations in optics and acoustics where there is cylindrical symmetry, for example the vibrations of a drum head. The derivatives of Bessel functions are other Bessel functions just like the derivatives of sines and cosines are other sines and cosines. 

Hint: a quick way to plot a function and its derivative is to use this syntax: 
f = function_of_x
df = D[f,x]
Plot[{f,df}, various_plot_options]

Here you don’t have to use pure argumentative functions (f[x_]=...), because simple algebraic expressions will work with the Plot command.

Something that won’t work when trying to plot a derivative, however, is sticking the derivative command inside the Plot command, like this command shown with sin(x) as the function for demonstration purposes: 
Plot[D[Sin[x],x]}, {x,0,10}]

If you try that command, you’ll get error messages such as “0.0002042857142857143` is not a valid variable.”  That’s because it’s choosing x=0.0002042857142857143 as a point for plotting the function before it tries to calculate the derivative of sin(x) with respect to x. So it’s attempting to calculate D[Sin[0.0002042857142857143],0.0002042857142857143], which clearly makes no sense.

(c) Plot (Sqrt[1+(x^2)]) along with ((d^3)/(dx^3))(Sqrt[1+(x^2)]) , from -4 to 4. Plot the function in blue and the derivative in red. 
(d) Plot3D cos(x+y) along with (∂/(∂y))((∂/(∂x))cos(x+y)), from -π to π on both the x- and the y-axes. Plot the function in blue and the derivative in red. Note that if you're unfamiliar with partial derivatives, the curvy d in the d/dx sign means "take the derivative with respect to x, while treating y as a constant". Same sort of thing for the curvy d/dy.  Shortcut notation for that is ((∂^2)/(∂x∂y))cos(x+y) or ((∂^2)/(∂y∂x))cos(x+y)  (the order that you take the x and y derivatives generally does not matter).

Applications

Electric field from point charges

In Physics 220 you’ve learned that the electric potential from a point charge at the origin is V(x, y, z) =(1/4(Subscript[πϵ, 0])) q/r, where r is the distance from the origin to the point at which you want to know the field, r=((((x^2)+(y^2)+(z^2)))^0.5). Similarly, the electric potential at the point (x,y,z) due to a charge located at (x0,y0,z0) is the same equation, with x replaced by x-x0, y replaced by y-y0, and z replaced by z-z0. 

You may not yet have learned that the components of the electric field can always be calculated from the electric potential via: (Subscript[E, x])=-((δV)/(δx)), (Subscript[E, y])=-((δV)/(δy)), and (Subscript[E, z])=-((δV)/(δz)). In other words, the electric field is the negative gradient of the electric potential.

Assignment 5. Calculating the field from the potential: electric dipole

Consider the same electric dipole we looked at in the Plotting lab. For simplicity, again set 1/(4(Subscript[πϵ, 0])) equal to 1. And we’ll again set z = 0 to just focus on the x-y plane.
(a) Define a function (or expression) that gives you the electric potential for two charges of opposite signs, a positive charge q = +1 at x = 0.5 and a negative charge q = -1 at x = -0.5. Remember to shift the potential from each charge away from the origin using the guidance from the introductory paragraph above (x replaced by x-x0, etc.).
(b) Define the components of the electric field as the negative derivatives of the potential. Verify that StreamPlot gives you the same plot of the field lines it did in the last lab. Review the Plotting lab and/or consult the help file for the StreamPlot command to remind yourself what the input arguments need to be.
(c) Make a ContourPlot of the electric potential in the xy plane with ContourShading->None, and plot it on the same graph as the StreamPlot (the Show command will be helpful). You should be able to verify that the electric field lines are always perpendicular to the equipotential contours, as taught in Physics 220.

Assignment 6. Ten random charges

Now that you know what you're doing, you should be able to fairly easily modify your code to use the potential function given below. It will give you the potential due to 10 random charges (positive and negative) at random positions. Plot the equipotential contours along with the field lines. Make two plots: one so you can see things close to the charges (say, x and y both go from -1 to 1), and one so you can see things far from the charges (say, x and y both go from -10 to 10). Why do the potential and field lines behave the way they do on the “zoomed out” graph?

Vbase = q/Sqrt[x^2+y^2];
V = Sum[Vbase /. {q-> RandomReal[{-1,1}],x-> x-RandomReal[{-1,1}], y-> y-RandomReal[{-1,1}]}, {n,1,10}];

Air resistance

The projectile equations,  x(t)=(Subscript[x, 0])+(Subscript[v, 0])cos(θ)t and y(t)=(Subscript[y, 0])+(Subscript[v, 0])sin(θ)t- (1/2)g (t^2) arise from an analysis of an object in free-fall when gravity in the y-direction is the only force acting on an object. Thus (Subscript[a, x])=0 and (Subscript[a, y])=-g. In Physics 121 you likely completely neglected air resistance, aka “viscous drag”, because it complicates the equations tremendously; however, with Mathematica we can fairly easily add the effects from drag into the equations. Let’s suppose there is some drag coefficient b which we know, a number which would depend on the nature of the fluid (e.g. air vs. water) as well as the nature of the object moving through the fluid (such as the cross-section area). In many cases the drag force is proportional to the velocity: (Subscript[F, drag])=-b v. (The units of b are kg/s.)

Let’s consider only 1D motion, such as throwing a ball straight up in the air. Newton’s second law, (Subscript[F, net]) = ma, is:   
	-m g-b (( dy)/dt)=m (((d^2)y)/(d (t^2))) 

Mathematica can solve such an equation (called a “differential equation”) using the DSolve command. You just need to tell DSolve the equation to be solved, along with the initial position y[0] and velocity y'[0] (all three in a list). Since many/most of you may not yet be familiar with differential equations, don’t worry about learning how to use DSolve; just evaluate this box to get the properly defined function for y(t) for a projecticle with air resistance.

Clear["`*"]   (* just in case you've defined y as something that it won't like *)
DSolve[  {-m g -b y'[t] == m y''[t] , y[0]==y0, y'[0]==v0}  , y[t], t]
y[t_] = y [t]/. %[[1]]  //Expand   (* this is needed to define y as a function of t because the output of DSolve is in a weird format *)

Assignment 7. Projectile with air resistance

(a) Test to make sure that y(t) function makes sense by evaluating it in the limit as b-> 0. It should match the Phys 121 equation. 
(b) Make plots of y(t), v(t), and a(t) for a 1 kg object, starting at y0 = 0 and with initial velocity v0 = 100 m/s. (Remember (Subscript[a, y])=-g, so g =+9.8 (m/(s^2)).)  Let t range from 0 to 20 s. Make one set of plots for b = 0.00001 kg/s, and another set for b = 0.3 kg/s. 
(c) For the condition where b = 0.3, find the time it takes the object to reach the maximum height, the time it takes to return to the ground, and the velocity it has when it reaches the ground (note how kinetic energy is not conserved!).

Hint: Here's a slick way to generate three plots at once, assuming f1(t), f2(t), and f3(t) are previously-defined pure functions:  Plot[# ,{t,0,20}] & /@{f1[t], f2[t], f3[t]}.   Or if f1, f2, and f3 are expressions, you can do it this way: Plot[# ,{t,0,20}] & /@{f1, f2, f3}.  It creates a three element array, each element of which is a plot. Remember the /@ notation refers to the Map function, and is extremely helpful in getting another function (e.g. Plot) to be evaluated for each of the elements in a list.

Maxwell-Boltzmann velocity distribution

In Phys 123 you learned about the Maxwell-Boltzmann velocity distribution which describes the distribution of speeds in an ideal gas:  f(v)=sqrt((2/π)((((m/kT)))^3))(v^2)exp(-((mv^2)/(2kT))).  In the equation T is the temperature, k is Boltzmann’s constant, m is the mass of an individual molecule, and the function f describes the probability a given molecule will be at speed v. You should also have learned about these three special speeds: the most probable speed: (Subscript[v, mp])=sqrt(((2kT)/m));  the average speed: (Subscript[v, ave])=sqrt(((8kT)/(π m))); and the root mean square speed, i.e. the square root of the average (v^2) value: (Subscript[v, rms])=sqrt(((3kT)/m)).

Assignment 8. Using the Maxwell-Boltzmann velocity distribution

(a) To get a better feel for the Maxwell-Boltzmann distribution, make a plot of it for nitrogen molecules at temperatures of 77K (the condensing point), 300K (room temperature), and 3000K (roughly the temperature of the filament inside an incandescent light bulb). Hint: the molecular mass of nitrogen molecules is 28.01 AMU; you’ll need to convert that to kg.
(b) The phrase “most probable speed” refers to where the f(v) distribution is peaked. Without using any numbers at all now, take the derivative of f(v), set it equal to zero, and solve for v. You should obtain the (Subscript[v, mp]) formula given above.

Wrap-up

The following functions were discussed in today's lab, or are very closely related:
Mathematica calculus-related commands:  Limit, D (shortcut  ' ),  Series (often with Normal)
Interesting mathematical functions: HeavisideTheta, Fibonacci, BesselJ, Erf
Statistical functions: RandomReal, MaxwellDistribution
Plot functions: StreamPlot, ContourPlot
Misc:  /@ is the Map function
Used in the lab, but you don’t need to learn for this course: DSolve

Potentially useful summary/tutorial pages:
(I don’t believe there were any today.)