Using Taylor Series to check the behavior of a function in different limits

Using a Coulomb integral to determine the potential or electric field for a distribution of charge can be tricky. How do we know we got the right answer?

One of the first things we should do with a new result is make sure it is consistent with things we already know.

If you think you understand the potential at a distance r from a blob of charge, does your result look like the potential for a point charge when r is very large? What if you approach very close to the surface? Is the behavior what you’d expect in that case?

We can use a Taylor series to see what our result looks like when some quantity in our result is very big or very small compared to other quantities. For example, if the charge distribution’s size is characterized by a length L, then the potential or electric field should behave like a point charge when r>>L, which is the same thing as L/r << 1. How do we do this in Mathematica?

Let’s use the uniformly charged disk example from class. The constant surface charge density is σ=Q/(π (R^2)), and the potential at a point above or below the center of the disk is

V[0,0,z] = (σ/(2(Subscript[ϵ, 0])))((Sqrt[(R^2)+(z^2)])-(Sqrt[(z^2)]))= (Q/(2π (Subscript[ϵ, 0])(R^2)))((Sqrt[(R^2)+(z^2)])-(Sqrt[(z^2)]))

Let's define a function V[z] that describes this potential. First, let me tell Mathematica that R (the radius of the disk) is a positive quantity:

$Assumptions = R>0 ;

We don't have to do this. But Mathematica doesn’t assume that quantities are positive or negative or even real, so stating this now will let Mathematica simplify things like (Sqrt[(R^2)])=R without us always having to remind it. Now let’s define our function:

V[z_]:= (Q/(2π (Subscript[ϵ, 0])(R^2)))((Sqrt[(R^2)+(z^2)])-(Sqrt[(z^2)]))

The notation “:=” assigns a relationship between V[z] and the expression on the right-hand side. The underscore in the argument – writing “V[z_]” instead of “V[z]” – tells Mathematica that z could be anything. That way, when we invoke our function, Mathematica will replace z in the definition with whatever we put in square brackets. Defining the function this way lets us evaluate it for different values of z.

V[0]//Simplify

V[R]//Simplify

I told Mathematica to Simplify the result, so it takes advantage of knowing R>0 and replaces (Sqrt[(R^2)]) with R.

Next, let’s check that this potential reproduces the electric field we found at points along the z axis. In general, the electric field is minus the gradient of the potential. We could take the gradient here, but since we only know the potential along a line above and below the disk, we don’t actually know dV/dx or dV/dy. It’s better to just evaluate d/dz of our potential, which we can trust to tell us the electric field at points along the z axis.

-D[V[z],z]//Simplify

This is the result we got in class for (Subscript[E, z])(0,0,z). Notice that Mathematica doesn’t know anything about z – it doesn’t even know that z is a real number. Since z is ambiguous it doesn’t try to replace (Sqrt[ (z^2)]) with z. We can tell it that z is real, in which case it replaces(Sqrt[(z^2)]) with z and z/(Sqrt[(z^2)]) with Sign[z]=±1.

Assuming[z ∈ Reals,-((Q z (-(1/(Sqrt[(z^2)]))+(1/(Sqrt[(R^2)+(z^2)]))))/(2 π (R^2) (Subscript[ϵ, 0])))//Simplify]

Or we can tell Mathematica that z is positive or negative:

Assuming[z >0,-((Q z (-(1/(Sqrt[(z^2)]))+(1/(Sqrt[(R^2)+(z^2)]))))/(2 π (R^2) (Subscript[ϵ, 0])))//Simplify]

Assuming[z <0,-((Q z (-(1/(Sqrt[(z^2)]))+(1/(Sqrt[(R^2)+(z^2)]))))/(2 π (R^2) (Subscript[ϵ, 0])))//Simplify]

We could have declared any of these at the beginning when we specified that R>0, but being too specific can clutter Mathematica’s output with little warning and caveats. We’ll just ask it to treat z as positive or negative as needed.

V[z]//FullSimplify

If we explicitly tell Mathematica that z is positive or negative, it replaces (Sqrt[(z^2)]) with ±z, respectively.

Assuming[z>0,V[z]//FullSimplify]

Assuming[z<0,V[z]//FullSimplify]

So how do we check the behavior of V[z] for z>>R? I can use the built-in “Series” command to perform a Taylor series for any variable, around any point, to any order. For example, I can ask Mathematica for just the leading term when z approaches +∞:

Series[V[z],{z,∞,1}]

This is the correct result for a point charge.  Since z was approaching positive infinity, I got a result appropriate for positive z. If instead I expand around z-> -∞, I get an expression appropriate for negative z.

Series[V[z],{z,-∞,1}]

The “((O[1/z])^2)” term is an object that Mathematica tracks and understands. I can get rid of it with the “Normal” command. The Series command can go inside “Normal[...]”

Normal[Series[V[z],{z,∞,1}]]

Or I can use “postfix” notation to apply Normal to the output of Series. Either way, I get the same thing.

Series[V[z],{z,∞,1}]//Normal

We can also ask Mathematica to give us more terms in the expansion. When we expand z around infinity, Mathematica understands that we’re asking for terms out to a certain order in 1/z.

Series[V[z],{z,∞,5}]

We can also expand around a finite value of z, like z=0. Here are the leading term and the first sub-leading term (the next most important term).

Series[V[z],{z,0,1}]//Normal

Since we haven’t specified anything about z, Mathematica hasn’t simplified any further.

What if we want to Taylor expand the electric field around z=0?

-Grad[V[z],{x,y,z}]//FullSimplify

Series[-((Q z (-(1/(Sqrt[(z^2)]))+(1/(Sqrt[(R^2)+(z^2)]))))/(2 π (R^2) (Subscript[ϵ, 0]))),{z,0,0} ]//Normal

Tell Mathematica that z is a real number and simplify the result, so we explicitly get plus or minus for the z component of the electric field on either side of the disk.

Assuming[z∈ Reals,((Q z)/(2 π (R^2) (Sqrt[(z^2)]) (Subscript[ϵ, 0])))//FullSimplify]

Is the leading term correct? Let’s replace Q with σ*π (R^2) (the surface charge density times the area of the disk) to compare with the result we got in class:

(((Q Sign[z])/(2 π (R^2) (Subscript[ϵ, 0])))/.Q-> σ*π*(R^2))//Simplify

This is correct! The field is in opposite directions on either side of z=0. Mathematica even knows how to take the limit from either side so we can check the discontinuity in the field.

Limit[{0,0,((Q z)/(2 π (R^2) (Sqrt[(z^2)]) (Subscript[ϵ, 0])))},z->0,Direction->"FromAbove"]-Limit[{0,0,((Q z)/(2 π (R^2) (Sqrt[(z^2)]) (Subscript[ϵ, 0])))},z->0,Direction->"FromBelow"]/.Q-> σ*π*(R^2)

This is correct result. The discontinuity in the electric field is in the (OverHat[z]) direction, and the value on the two sides differs by σ/(Subscript[ϵ, 0]).