Supplement to Lab 6: Legendre PolynomialsRecall from recent lectures that Legendre's diferential equation is:restart; with(plots):assume(n, integer);PkkpbGVnZW5kcmVHNiIvLUklZGlmZkclKnByb3RlY3RlZEc2JComLCYiIiJGLCokSSJ4R0YkIiIjISIiRiwtRic2JC1JInlHRiQ2I0YuRi5GLEYuLCQqKEkibkdGJEYsLCZGOEYsRixGLEYsRjNGLEYwor, in Maple formatlegendre := diff((1-x^2)*diff(y(x),x),x) = -n*(n+1)*y(x);Let's have Maple solve it instead of going once again through all the work we did in class. We do that by typing:dsolve(legendre,y(x));This is telling us that we can form a solution to Legendre's equation by taking arbitrary linear combination of the two functions LegendreP and LegendreQ. As long as n is an integer, the first one of these is a relatively simple function: it's just a polynomial. If n is not an integer, things aren't so nice. In fact, the function called Q is pretty messy even if n is an integer. Let's take a look at some of them.plot({LegendreP(0,x),LegendreP(1,x),LegendreP(2,x),LegendreP(3,x),LegendreP(4,x)},x=-1..1);If n is not an integer, the Legendre polynomials go to infinity at x = \342\200\2231. plot({LegendreP(1/2,x),LegendreP(3/2,x)},x=-1..1);It means that the Legendre functions for non-integer n are not useful in the physical applications we'll be considering. Why? Well, typically the argument x in the Legendre polynomail is x = cos(\316\270) where \316\270 is a "latitude"-like variable. We don't expect physical properties to necessarily blow up to infinity at the "poles" (when cos(\316\270) = \342\200\2231).Let's look at the Q functions:plot({Re(LegendreQ(1,x)),Re(LegendreQ(2,x))},x=-1..1);They all go to infinity at x = \302\2611, even when n is an integer! So, although they're mathematically consistent solutions to Legendre's equation, they're not generally useful in physical applications. Henceforth, we will not discuss them.The first couple of Legendre functions are pretty simple:LegendreP(0,x);LegendreP(1,x);The next one is just a quadratic. Maple's default mode is to keep the next one hidden when you type the following:LegendreP(2,x);but we can force it to display the function:taylor(LegendreP(2,x),x);
If we want to, we can force Maple to write out explicitly all of the terms in any Legendre polynomial by declaring a one-line procedure as we have done several times before:LegendrePexplicit:= n-> convert(taylor(LegendreP(n,x),x,n+1),polynom);
An now we can explicitly show all the Legendre polynomials:LegendrePexplicit(2);LegendrePexplicit(3);LegendrePexplicit(4);LegendrePexplicit(25);
Actually, though, we don't need to go to all this trouble. Maple knows a more efficient method of calculating the Legendre polynomials. It's contained in a package called "orthopoly":with(orthopoly):P(2,x);P(3,x);P(4,x);P(25,x);plot({P(0,x),P(1,x),P(2,x),P(3,x),P(4,x),P(5,x),P(6,x),P(7,x),P(8,x),P(9,x)},x=-1..1);
For this week's lab exercise you are going to be using Legendre series to fit the "Heaviside" function, which is basically a function which takes on a value of zero for x<0 and values of one for x>0. So, for example:heavy := Heaviside(x); plot(heavy,x=-1..1,thickness=2);
And now some physics:
As discussed in class, Legendre polynomials come up when dealing with spherically symmetric potentials (such as those for gravity, electromagnetism, etc.). The electric potential for a charge distribution can be written as a series of terms containing Legendre Polynomials. The nth term in such a series (starting with the n=0 term) can be written:V:= n-> r^(-(n+1))*P(n,cos(theta));and get a solution to the first few terms in this infinite series. Here are the first few. The n=0 (monopole) term is proportional to the electric potential due to a point charge at the origin. V(0); The n=1 term is proportional to the potential due to a dipole.V(1); The n=2 term is proportional to the potential due to a quadropole.V(2); Let's see what that these potentials looks like. UNfortunately they are in spherical coordinates and the plotting routines prefer cartesian (x,y,z) coordinates.
Therefore we must first write a multi-line procedure, called convert_to_xyz which will take something written in spherical coordnates and it will make the substitutions:
r = LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYlLUkjbWlHRiQ2I1EhRictSSZtc3FydEdGJDYjLUYjNigtSSVtc3VwR0YkNiUtRiw2JVEieEYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYkLUkjbW5HRiQ2JFEiMkYnL0Y+USdub3JtYWxGJ0ZGLyUxc3VwZXJzY3JpcHRzaGlmdEdRIjBGJy1JI21vR0YkNi1RIitGJ0ZGLyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR0ZRLyUpc3RyZXRjaHlHRlEvJSpzeW1tZXRyaWNHRlEvJShsYXJnZW9wR0ZRLyUubW92YWJsZWxpbWl0c0dGUS8lJ2FjY2VudEdGUS8lJ2xzcGFjZUdRLDAuMjIyMjIyMmVtRicvJSdyc3BhY2VHRmpuLUY1NiUtRiw2JVEieUYnRjpGPUZARkhGSy1GNTYlLUYsNiVRInpGJ0Y6Rj1GQEZIRkZGRg==\316\270 = LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYlLUkjbWlHRiQ2JVEnYXJjY29zRicvJSdpdGFsaWNHUSZmYWxzZUYnLyUsbWF0aHZhcmlhbnRHUSdub3JtYWxGJy1JKG1mZW5jZWRHRiQ2JC1GIzYkLUkmbWZyYWNHRiQ2KC1GLDYlUSJ6RicvRjBRJXRydWVGJy9GM1EnaXRhbGljRictRiM2JS1GLDYjUSFGJy1JJm1zcXJ0R0YkNiMtRiM2KC1JJW1zdXBHRiQ2JS1GLDYmUSJ4RidGQC8lMGZvbnRfc3R5bGVfbmFtZUdRKDJEfk1hdGhGJ0ZCLUYjNiQtSSNtbkdGJDYlUSIyRidGVEYyRjIvJTFzdXBlcnNjcmlwdHNoaWZ0R1EiMEYnLUkjbW9HRiQ2LlEiK0YnRlRGMi8lJmZlbmNlR0YxLyUqc2VwYXJhdG9yR0YxLyUpc3RyZXRjaHlHRjEvJSpzeW1tZXRyaWNHRjEvJShsYXJnZW9wR0YxLyUubW92YWJsZWxpbWl0c0dGMS8lJ2FjY2VudEdGMS8lJ2xzcGFjZUdRLDAuMjIyMjIyMmVtRicvJSdyc3BhY2VHRl5wLUZPNiUtRiw2JlEieUYnRkBGVEZCRldGZ25Gam4tRk82JS1GLDYmRj9GQEZURkJGV0ZnbkYyRjIvJS5saW5ldGhpY2tuZXNzR1EiMUYnLyUrZGVub21hbGlnbkdRJ2NlbnRlckYnLyUpbnVtYWxpZ25HRl9xLyUpYmV2ZWxsZWRHRjFGMkYyRjI=\317\206 = arctan (y/x)
and then tell Maple to simplify the results to get an expression in for a potential in spherical coordinates into Cartesian (x,y,z) coordinates. convert_to_xyz:= proc (f,r,theta,phi,x,y,z)local g;g:= eval(eval(f,{r=sqrt(x^2+y^2+z^2),theta=arccos(z/sqrt(x^2+y^2+z^2)),phi=arctan(y,x)}));g:=simplify(g);return(g);end;Let's convert the dipole potential V(1) into cartesian coordinates:V1xyz:=convert_to_xyz(V(1),r,theta,phi,x,y,z);
To get the electric field corresponding to this potential, we need to take the negative of the gradient (don't worry if you haven't seen this, its in E&M, in a nutshell to turn an electric potential into a electric field, you have to take the derivative of it).with(linalg):E1xyz:=grad(-V1xyz,[x,y,z]);fieldplot3d(E1xyz,x=-10..10,y=-10..10,z=-10..10,axes=boxed,arrows=THICK,scaling=constrained);Hmmm. It's hard to tell what's going on there. The field is so strong near the origin that all of the arrows except the few closest ones are too small to see.
Here's one way around that: let me make all of the points within a sphere of radius five equal to zero and then replot this:Ecut:= evalm(Heaviside(x^2+y^2+z^2-5^2)*E1xyz);fieldplot3d(Ecut,x=-10..10,y=-10..10,z=-10..10,axes=boxed,arrows=THICK,scaling=constrained);