Solution of Laplace's equation using interpolation functions rather than FEM / FEA

Ken Hill
kohill@rogers.com (613) 592-2809

Ref. Silvester and Ferrari "Finite Elements for Electrical Engineers" 3rd Edition

Needs["Miscellaneous`Units`"]
Needs["Miscellaneous`PhysicalConstants`"]
Needs["Miscellaneous`SIUnits`"]
<< Calculus`VectorAnalysis`
<< Graphics`MultipleListPlot`
<< Graphics`Legend`

The potential of two long parallel wires a distance 2a apart with charge densities ±q,
respectively, is given by
U(x, y) =(q / 2π) ln( ((x + a)^2 + y^2) / ((x - a)^2 + y^2 ) )

Turn off warning messages triggered by similar variable names.

Off[General::spell];
Off[General::spell1];

SetCoordinates[Cartesian[x, y, z]]

dim1 = 50;
dim2 = 50;

The potential of two long parallel wires a distance 2a apart with charge densities ±q

ThePot[x_, y_, a_, q_] := (q/(2  π)) 
        Log[((x + a)^2 + y^2)/((x - a)^2 + y^2)]

The symbolic "Potential" (note that some of the a's are changed to b's; in this example two a's are changed, each to represent a wire)

fn = 
        Table[{i, j, a[i, j]}, {i, 1, dim1}, {j, 1, 
              dim2}] /. {a[23, 25] -> b[23, 25], 
            a[27, 25] -> b[27, 25]};

Dimensions[fn]

The Flattened "Potential"

ne = Flatten[fn, 1];

Selecting all the a's data (at some point some of the a[i,j]s will be inputed as b[i,j]s, i.e. constants -- see above)

avalues = 
        Cases[ne, ({_, _, a[x_, y_]} -> a[x, y])];

The interpolation function on the potential

re = Interpolation[ne, InterpolationOrder -> 1]

Gradient of re using Mathematica

Grad[re[x, y]]

Dot product of the gradient of re with itself

ss[x_, y_] = DotProduct[Grad[re[x, y]], Grad[re[x, y]]]

We work out the gradient and evaluate the integral of the dot product of the gradient with itself (see Silvester)

Using the well-known principle of minimum potential energy we minimize the stored energy over the whole region. We need to integrate the dot product of the gradient with itself to compute the stored energy.

We have been able to work out the stored energy over the whole region by defining a new interpolation function over the dot product of the gradient of the potential with itself and integrating the function.

Generate a new interpolation function for the dot product of the gradient with itself (proportional to the energy density) -- we should then be able to integrate it symbolically

energyFN = 
        Table[{i, j, ss[i, j]}, {i, 1, dim1}, {j, 1, dim2}];

flatenergyFN = Flatten[energyFN, 1];

energyfunction = 
      Interpolation[flatenergyFN, InterpolationOrder -> 1]

StoredEnergy = 
        Integrate[
          energyfunction[x, y], {x, 1, dim1}, {y, 1, dim2}];

The equations for determining the potential values

TheEquationsForExtremum = 
        Table[D[StoredEnergy, avalues[[i]]] == 0, {i, 
            1, Dimensions[avalues][[1]]}];

Lets give values to the b's -- in this case we are looking at a dipole

b[23, 25] = 1
b[27, 25] = -1

TheSol = 
        Flatten[Solve[TheEquationsForExtremum, avalues], 1];

Let's Plot the interpolation function on the potential

DensityPlot[re[x, y] /. TheSol, {x, 1, dim1}, {y, 1, dim2}, 
      PlotPoints -> 36]

ContourPlot[re[x, y] /. TheSol, {x, 1, dim1}, {y, 1, dim2}, 
      PlotPoints -> 36]

The potential looks OK

Let's plot the explicit form of the potential function

DensityPlot[
      ThePot[x - dim1/2, y - dim2/2, a, q] /. {a -> 2, 
          q -> -1}, {x, 1, dim1}, {y, 1, dim2}, 
      PlotPoints -> 36]

ContourPlot[
      ThePot[x - dim1/2, y - dim2/2, a, q] /. {a -> 2, 
          q -> -1}, {x, 1, dim1}, {y, 1, dim2}, 
      PlotPoints -> 36]

Plot along the x-axis of the numerical result (re) and the explicit result -- note: the explicit result has singularities at (dim1 / 2 +/-  a )

Plot[{re[x, dim2/2] /. TheSol, 
        ThePot[x - dim1/2, 25 - dim2/2, a, q] /. {a -> 2, 
            q -> -1}}, {x, 1, dim1}, 
      PlotStyle -> {RGBColor[0, 0, 1], RGBColor[1, 0, 0]}]

Plot along the x-axis of the numerical result (re) and the explicit result (with a little offset of y to fit the peaks)

Plot[{re[x, dim2/2] /. TheSol,ThePot[x - dim1/2,25.16-dim2/2,a,q]/.{a -> 2, q -> -1}},{x, 1, dim1},PlotStyle -> {RGBColor[0, 0, 1], RGBColor[1, 0, 0]}]