(*^
::[	Information =

	"This is a Mathematica Notebook file.  It contains ASCII text, and can be
	transferred by email, ftp, or other text-file transfer utility.  It should
	be read or edited using a copy of Mathematica or MathReader.  If you 
	received this as email, use your mail application or copy/paste to save 
	everything from the line containing (*^ down to the line containing ^*)
	into a plain text file.  On some systems you may have to give the file a 
	name ending with ".ma" to allow Mathematica to recognize it as a Notebook.
	The line below identifies what version of Mathematica created this file,
	but it can be opened using any other version as well.";

	FrontEndVersion = "Macintosh Mathematica Notebook Front End Version 2.2";

	MacintoshStandardFontEncoding; 
	
	fontset = title, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeTitle, center, M7, bold, e8,  24, "Times"; 
	fontset = subtitle, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeTitle, center, M7, bold, e6,  18, "Times"; 
	fontset = subsubtitle, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeTitle, center, M7, italic, e6,  14, "Times"; 
	fontset = section, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, grayBox, M22, bold, a20,  18, "Times"; 
	fontset = subsection, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, blackBox, M19, bold, a15,  14, "Times"; 
	fontset = subsubsection, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, whiteBox, M18, bold, a12,  12, "Times"; 
	fontset = text, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  14, "Times"; 
	fontset = smalltext, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  10, "Times"; 
	fontset = input, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeInput, M42, N23, bold, L-5,  12, "Courier"; 
	fontset = output, output, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, L-5,  12, "Courier"; 
	fontset = message, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, R65535, L-5,  12, "Courier"; 
	fontset = print, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, L-5,  12, "Courier"; 
	fontset = info, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, B65535, L-5,  12, "Courier"; 
	fontset = postscript, PostScript, formatAsPostScript, output, inactive, noPageBreakInGroup, nowordwrap, groupLikeGraphics, M7, l34, w405, h197,  12, "Courier"; 
	fontset = name, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, italic,  10, "Geneva"; 
	fontset = header, inactive, noKeepOnOnePage, preserveAspect, M7,  12, "Times"; 
	fontset = leftheader, inactive, L2,  12, "Times"; 
	fontset = footer, inactive, noKeepOnOnePage, preserveAspect, center, M7,  12, "Times"; 
	fontset = leftfooter, inactive, L2,  12, "Times"; 
	fontset = help, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  10, "Times"; 
	fontset = clipboard, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  12, "Times"; 
	fontset = completions, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  12, "Times"; 
	fontset = special1, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  12, "Times"; 
	fontset = special2, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  12, "Times"; 
	fontset = special3, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  12, "Times"; 
	fontset = special4, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  12, "Times"; 
	fontset = special5, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  12, "Times"; 
	paletteColors = 128; currentKernel; 
]
:[font = section; inactive; Cclosed; noPageBreak; preserveAspect; cellOutline; backColorRed = 58981; backColorGreen = 58981; backColorBlue = 58981; startGroup]
Example  1: Potential and a Rectangular Groove.  V[x,0]=x
:[font = text; inactive; noPageBreak; preserveAspect]
Consider a rectangular groove that runs from z=-Infinity to +Infinity and  is open in the  positive y direction. The groove is bounded by three walls: y=0, x=0 and x=a.   The walls at x=0 and x=a are at zero potential and the end at y=0 is  at a specified potential V[x,0] (independent of the z-coordinate).
      
a) Find a series expansion for the potential inside the groove for an arbitrary V[x,0].

b)  Write the first few terms of the potential if V[x,0]=x.
 Plot the potentials, and use VEPlot to plot the equipotential and electric field lines. 
 
c) Sum the series for V[x,0]=x and   get a closed form for the solution.
:[font = subsubsection; inactive; noPageBreak; preserveAspect]
Remarks and Outline:
:[font = text; inactive; noPageBreak; preserveAspect]
1)   Use Cartesian coordinates and express the answer as a product X[x] Y[y].  
2)   The boundary conditions  suggest that X[x] is proportional to  Sin, and  Y[y]  is proportional to a decreasing exponential. The orthogonality properties of Sin  allow the expansion coefficients to be expressed as integrals over the boundary at y=0.
3) For V[x,0]=x the command  Sum will analytically sum the infinite potential series.
:[font = text; inactive; noPageBreak; preserveAspect]
Using Mathematica's ability to work with symbolic integrals, we can actually solve an entire class of problems with only two lines of code.  Having solved the two 1-dimensional problems, we can reconstruct these back into a single 2-dimensional problem, and then graph the solution to verify its accuracy.  We show intermediate results using series with varying numbers of terms. 
:[font = text; inactive; noPageBreak; preserveAspect]
Compare the accuracy of the solutions with the specified boundary conditions as a function of the number of terms used in the series. Note that it is particularly difficult to reproduce the "sharp corners" using sine and exponential expansions.  Experiment with other choices for the specified potential such as  Sin[n x] and Cos[n x]. 
:[font = input; noPageBreak; preserveAspect]
Clear["Global`*"];
:[font = subsubsection; inactive; noPageBreak; preserveAspect; startGroup]
Required Packages:
:[font = input; noPageBreak; preserveAspect; endGroup]
Needs["Algebra`SymbolicSum`"] 
Needs["Algebra`Trigonometry`"] 
Needs["Calculus`VectorAnalysis`"]
Needs["Graphics`PlotField`"] 
:[font = subsubsection; inactive; Cclosed; noPageBreak; preserveAspect; startGroup]
Solution
:[font = text; inactive; noPageBreak; preserveAspect; cellOutline]
a) Find a series expansion for the potential inside the groove for an arbitrary V[x,0].
:[font = text; inactive; noPageBreak; preserveAspect]
a) The boundary conditions suggest that we use  Cartesian coordinates, and that the solution is independent of z. We can express the potential as an expansion of orthogonal functions of variables x and y 
:[font = input; noPageBreak; preserveAspect]
pot[f_,n_]:= Sum[A[f,m ] Sin[Pi m x/a] Exp[-Pi m y/b],{m,1,n}]
:[font = text; inactive; noPageBreak; preserveAspect]
with coefficients A[f,m] which are determined using the orthogonality property of the Sin function:
:[font = input; noPageBreak; preserveAspect]
A[f_,m_]:= 2/a Integrate[f Sin[Pi m x/a],{x,0,a}]
:[font = text; inactive; noPageBreak; preserveAspect]
The form of the potential is chosen to ensure that it vanishes at x=0 and x=a.

 The function pot[f,n] gives the first n terms of the potential if the end wall at y=0 has a potential V[x,0]. (The operator pot only works for the special case that the integral in A can be done analytically; otherwise, Integrate must be replaced by NIntegrate.)
:[font = text; inactive; noPageBreak; preserveAspect; cellOutline]
b)  Write the first few terms of the potential if V[x,0]= x.
 Plot the equipotential and electric field lines. 
:[font = text; inactive; noPageBreak; preserveAspect; startGroup]
Consider the case when the potential on the end of the grove is  V[x,0]=x. The first 4   terms of the potential are  (this takes a while)
:[font = input; noPageBreak; preserveAspect]
pot1= pot[x,4]
:[font = text; inactive; noPageBreak; preserveAspect]
To see this visually, let us choose some values, and plot the potential in the {x,y} plane:
:[font = input; noPageBreak; preserveAspect]
values={a->1,b->1,V0->1};

Plot3D[ pot1 //.values //Evaluate ,{x,0,1},{y,0,0.5}];
:[font = text; inactive; noPageBreak; preserveAspect]
Clearly, the boundary conditions are satisfied on all surfaces except y=0.
If we increase the number of terms in the potential series, the solution at the y=0 surface is closer to the prescribed function.
:[font = text; inactive; noPageBreak; preserveAspect]
If we increase the number of terms, the potential at the y=0 surface is closer to the prescribed value:
:[font = input; noPageBreak; preserveAspect]
pot2=pot[x,15] //.values;
Plot3D[ pot2 ,{x,0,1},{y,0,0.5}];
:[font = input; noPageBreak; preserveAspect]
We now plot the equipotential lines and electric field lines.
;[s]
1:0,1;62,-1;
2:0,12,10,Courier,1,12,0,0,0;1,15,12,Times,0,14,0,0,0;
:[font = input; noPageBreak; preserveAspect]
plot1= ContourPlot[ pot2 //.values ,{x,0,1},{y,0,0.5},ContourShading->False];
:[font = input; noPageBreak; preserveAspect]
plot2= PlotGradientField[ -pot2 //.values ,{x,0,1},{y,0,0.5}
,Graphics`PlotField`ScaleFunction->(1&)];
:[font = input; noPageBreak; preserveAspect; endGroup]
Show[plot1,plot2];
:[font = text; inactive; noPageBreak; preserveAspect; cellOutline]
c) Sum the series for V[x,0]=x and   get a closed form for the solution.
:[font = text; inactive; noPageBreak; preserveAspect; startGroup]
The case V[x,0]=x can be summed exactly with the help of the package  Algebra`SymbolicSum`. 
(Note, Algebra`SymbolicSum` is obsolete in version 3.0.)
;[s]
3:0,0;93,1;149,0;150,-1;
2:2,16,12,Times,0,14,0,0,0;1,16,12,Times,2,14,0,0,0;
:[font = input; noPageBreak; preserveAspect]
pot3=((( pot[x,Infinity]
         //ComplexToTrig 
       ) //Simplify
         //TrigToComplex
       ) //Simplify
      )
:[font = text; inactive; noPageBreak; preserveAspect]
(Note, we have to help Sum  with TrigToComplex.)
We can see that the result agrees exactly with V[x,0]=1 on the y=0 boundary. 
:[font = input; noPageBreak; preserveAspect; endGroup; endGroup; endGroup]
Plot3D[ pot3 //.values //Evaluate ,{x,0.001,0.999},{y,0,0.5}];
:[font = section; inactive; noPageBreak; preserveAspect; cellOutline; backColorRed = 58981; backColorGreen = 58981; backColorBlue = 58981; startGroup]
Problem  2: Potential and a Rectangular Groove.  V[x,0]=1
:[font = text; inactive; preserveAspect; fontColorRed = 65535; endGroup]
Copy and modify from above. 
:[font = section; inactive; noPageBreak; preserveAspect; cellOutline; backColorRed = 58981; backColorGreen = 58981; backColorBlue = 58981; startGroup]
Problem  3: Potential and a Rectangular Groove.  V[x,0]=Sin[2 Pi x/a]^2
:[font = text; inactive; preserveAspect; fontColorRed = 65535; endGroup]
Copy and modify from above. 
:[font = section; inactive; Cclosed; noPageBreak; preserveAspect; cellOutline; backColorRed = 58981; backColorGreen = 58981; backColorBlue = 58981; startGroup]
Example  4:  Conducting Cylinder with a Potential on the Surface 
:[font = text; inactive; noPageBreak; preserveAspect]
Consider a long cylinder of radius r0 that is separated by a small lengthwise gap on each side.  Use cylindrical coordinates {r,phi,z} and let the length of the cylinder be along the z-axis. One-half of  the cylinder (0<phi<pi)  is kept at a potential V0 and the other half has a potential  -V0.  
a) Express the potential inside and outside the sphere as a power series in r and 1/r, respectively.
b) Plot the value of the potential as a function of theta for fixed r. Compare the series expansion at r=r0. 
c) Make a 3-dimensional plot of the potential as a function of r and theta.
:[font = subsubsection; inactive; noPageBreak; preserveAspect]
Remarks and Outline:
:[font = text; inactive; noPageBreak; preserveAspect]
      The solution of Laplace's equation can be represented by a product of two functions R[r] P[phi].  The symmetry of the problem requires P[phi] to be of the form Sin[m phi]. The term R[r] is of the form  {r^m, 1/r^m} depending on whether the field point is inside or outside the cylinder.  
:[font = text; inactive; noPageBreak; preserveAspect]
We now take advantage of Mathematica's ability to work in different coordinate systems. Again, we plot the result to verify its accuracy as a function of the terms in the series.  
:[font = text; inactive; noPageBreak; preserveAspect]
Compare the accuracy of the solutions with the specified boundary conditions as a function of the number of terms used in the series. Note that it is particularly difficult to reproduce the "sharp corners" using a finite number of terms in the expansion.

Experiment with other boundary potential functions such as Sin[n x] and Cos[n x].  For more complicated functions, replace Integrate by NIntegrate to generate purely numerical solutions. 

Try dividing the cylinder into 4 parts, and re-compute the expansion coefficients. 
:[font = input; noPageBreak; preserveAspect]
Clear["Global`*"];
:[font = subsubsection; inactive; noPageBreak; preserveAspect]
Required Packages:
:[font = input; noPageBreak; preserveAspect]
Needs["Calculus`VectorAnalysis`"]
Needs["Graphics`ParametricPlot3D`"]
Needs["Graphics`PlotField`"] 

:[font = text; inactive; noPageBreak; preserveAspect; cellOutline]
a) Express the potential inside and outside the sphere as a power series in r and 1/r, respectively.
:[font = text; inactive; noPageBreak; preserveAspect; startGroup]
The symmetry of the problem suggests we use cylindrical coordinates, and that the solution is independent of z. The solution of Laplace's equation can be represented by a product of two functions R[r] P[phi] where R[r] is of the form   r^m or r^-m,   and P[phi] is of the  form Sin[m phi].  (The dependence on the r variable is determined by requiring the terms in the potential to be finite.)
       Let us first construct the potential inside the cylinder. Inside the cylinder the potential must be finite at r=0 so  the  function R[r] is of the form r^m. The potential is:
:[font = input; noPageBreak; preserveAspect]
termIn[m_]= A[m] (r/r0)^m Sin[m phi] 
:[font = input; noPageBreak; preserveAspect]
potIn[n_]:= Sum[ termIn[m] ,{m,1,n}]
:[font = text; inactive; noPageBreak; preserveAspect]
We can determine the coefficients A[m] using the orthogonality properties of the sine function.  At the r=r0 boundary, the potential will be given by a sum of terms of the form:
:[font = input; noPageBreak; preserveAspect]
eq1= V[phi] == termIn[m] /.{r->r0}
:[font = text; inactive; noPageBreak; preserveAspect]
(There is an implied sum on the right-hand-side (rhs), {it i.e. V_0= sum_m A_m Sin(m phi).) We can determine A[m] by multiplying both sides of eq1 by Sin[n phi], and then integrating.  We do this first to the rhs. Since the sine functions are orthogonal, the only surviving term will be when n=m, so we spare Mathematica the effort of integrating all terms except this one. (That is, int , Sin(n phi) , Sin(m phi)=0 except for n=m.) The rhs is then:
:[font = input; noPageBreak; preserveAspect]
rhs= ( Integrate[ eq1[[2]] Sin[m phi] ,{phi,0,2 Pi}]
        //.{Sin[4 m Pi]->0}
      )
:[font = text; inactive; noPageBreak; preserveAspect]
Likewise, for the left-hand-side (lhs), we have:
:[font = input; noPageBreak; preserveAspect]
lhs=(  ( Integrate[ +V0 Sin[m phi] ,{phi, 0,1 Pi}] 
       + Integrate[ -V0 Sin[m phi] ,{phi,Pi,2 Pi}] 
     ) //.{Cos[m Pi]-> (-1)^m, Cos[2 m Pi]-> 1} 
       //Simplify
    ) 
       
:[font = text; inactive; noPageBreak; preserveAspect]
We have made use of the fact that the potential for  0<phi< pi is V0  and for pi< phi< 2 pi the surface potential is  -V0.
:[font = text; inactive; noPageBreak; preserveAspect]
Setting the lhs equal to the rhs, we find A[m].
:[font = input; noPageBreak; preserveAspect]
aSol= Solve[rhs==lhs,A[m]][[1]] //Simplify;
A[m_]= A[m] /. aSol 
:[font = text; inactive; noPageBreak; preserveAspect]
and thus, the potential:
:[font = input; noPageBreak; preserveAspect]
potIn[n]
:[font = text; inactive; noPageBreak; preserveAspect]
The potential outside follows in a similar fashion. The A[m] coefficients will match, but the expansion will be in powers of 1/r such that V[r->Infinity]=0
:[font = input; noPageBreak; preserveAspect]
termOut[m_]= A[m] (r0/r)^m Sin[m phi] 
:[font = input; noPageBreak; preserveAspect; endGroup]
potOut[n_]:= Sum[ termOut[m] ,{m,1,n}]
potOut[n]
:[font = text; inactive; noPageBreak; preserveAspect; cellOutline]
b) Plot the value of the potential as a function of theta for fixed r. Compare the series expansion at r=r0. 
:[font = text; inactive; noPageBreak; preserveAspect; startGroup]
We plot the value of the potential as a function of r and phi for r0=1 and V0=1 
:[font = input; noPageBreak; preserveAspect]
values={r0->1, V0->1};
:[font = text; inactive; noPageBreak; preserveAspect]
by combining the graphics for the inside and outside region. To do this, we creat a new function pot with the appropriate conditional properties:
:[font = input; noPageBreak; preserveAspect]
Clear[pot];
pot[n_,x_ /; x<1 ]:=potIn[ n];
pot[n_,x_ /; x>=1]:=potOut[n];
:[font = text; inactive; noPageBreak; preserveAspect]
To see the potential as a function of r, we plot a table of r values as a function of phi. First, we examine values for r<r0=1.
:[font = input; noPageBreak; preserveAspect]
Plot[Table[ pot[20,r/r0] //.values,{r,0,1,0.2}] //Evaluate
 ,{phi,0,2Pi}];
:[font = text; inactive; noPageBreak; preserveAspect]
Note, as r approaches r0, the potential approaches that on the surface. The effect of only summing a finite number of term in the series is obvious as we approach r=r0. Next, we examine values for r>r0=1.
:[font = input; noPageBreak; preserveAspect]
Plot[Table[ pot[20,r/r0] //.values,{r,1,2,0.2}] //Evaluate
 ,{phi,0,2Pi}];
:[font = text; inactive; noPageBreak; preserveAspect; endGroup]
Note, as r approaches Infinity, the potential approaches zero. The effect of only summing a finite number of term is again obvious.
:[font = text; inactive; noPageBreak; preserveAspect; cellOutline]
c) Make a 3-dimensional plot of the potential as a function of r and theta.
:[font = text; inactive; noPageBreak; preserveAspect; startGroup]
We can view the potential more directly using CylindricalPlot3D. 
:[font = input; noPageBreak; preserveAspect]
CylindricalPlot3D[ pot[20,r/r0] //.values
  ,{r,0,3}
  ,{phi,0,2Pi}
  ,BoxRatios->{1,1,1}
  ,ViewPoint->{4, -2, 0}
  ];
:[font = input; noPageBreak; preserveAspect; endGroup; endGroup]
CylindricalPlot3D[ pot[20,r/r0] //.values
  ,{r,0,1}
  ,{phi,0,2Pi}
  ,BoxRatios->{1,1,1}
  ,ViewPoint->{4, -2, 0}
  ,PlotPoints->30
  ];
^*)