ABSTRACT
Boundary problems with two closed
boundaries
Inseparable problems
Many problems in physics and engineering that exhibit boundaries that can
not be described by coordinate lines of a coordinate system in which the
actual partial dierential equation is separable On the other hand there are
inseparable partial dierential equations and the boundary conditions could
be described by any coordinate system In such a situation one can try to nd
formal solutions to the partial dierential equation in cartesian coordinates
Furthermore some problems may be nonlinear in nature and might become
separable by linearization We give an example
Twodimensional plasma equilibria are described by the Shafranov equa
tion If the magnetic ux surfaces are given a priori the equi
librium equations and the transport equations can be decoupled But for a
priori given cross sections of the ux surfaces the constraints to be imposed
on the Shafranov equation make the problem awkward because the pres
sure function p and the current j must be speci ed as a function of
whose spatial dependence r z is not yet known In the case of the linear
Shafranov equation particular solutions can be superposed and the equilib
rium problem may be reformulated as an eigenvalue problem for an arbitrary
cross section of the magnetic surfaces
The Shafranov equation reads in cylindrical coordinates r z
z
r
r
r
r
p
j
Here r z is the magnetic ux p is the pressure and j the current
distribution Making the linearizing setup
p
a
j
b
we obtain with
r z
L
X
A
l
cos k
l
z B
l
sin k
l
zR
l
r
R
l
r
R
l
b
k
l
R
l
r
R
l
where b is the same eigenvalue for several R
l
k
l
is the separation constant
which must not be an integer
a and the R
l
r are the degenerate
eigenfunctions belonging to the same eigenvalue b The partial dierential
equation is separable but the boundary conditions cannot be satis ed on
curves of the cylindrical coordinate system r z
If we designate the two solutions of by R
l
and R
l
for given k
l
then the ux surfaces const are given by
r z
L
X
l
C
l
R
l
r D
l
R
L
r
A
l
cos k
l
z B
l
sin k
l
z
If now the cross section zr of the outer magnetic surface plasma
surface is given a priori by a set of N pairs r
i
z
i
i N then
constitutes a system of N linear homogeneous equations for the L N un
known coecients A
l
C
l
A
l
D
l
C
l
B
l
B
l
D
l
For a set of given k
l
the vanishing
of the determinant of the coecients yields the eigenvalue b Initially b is not
known An estimate has to be made in the beginning The determinant does
not vanish and becomes a transcendental function of b By trial and error and
by repeated integration of an exact value b can be found Appar
entlyMathematica is not able to solve The relevant partial dierential
equation has been separable but the boundary conditions on the magnetic
surface determined by the solution of could not be satis ed on
surfaces of the cylindrical coordinates in which has been separable
A similar problem arises when one has to calculate the eigenfrequencies
of a membrane of arbitrary form The Helmholtz equation itself is sepa
rable in cartesian coordinates but an arbitrary membrane boundary cannot
be described by the straight lines of a cartesian coordinate system Again
collocation methods may help The Faber theorem is also useful any
membrane of arbitrary form but same area should have a rst eigenvalue near
but greater than of a circular membrane see This oers a
criterion concerning the eigenvalue We now investigate a membrane with a
boundary described by a Cassini curve This algebraic curve of fourth order
is described by
F x yx
y
c
x
y
a
c
c a a c
or in polar coordinates
r
r
c
cos
q
a
c
sin
For c becomes a circle of radius a We rst need the location
of the four vertex points xmax xmin ymax ymin where
p
p
numerical integration
Now let Mathematica do the work We wrote a full program that is an
extension of the program Colmeigv discussed in section First we verify
the solution of the membrane equation
and giving
proves the solution in cartesian coordinates Then we de ne the boundary
given by and respectively and plot it using the arbitrary values
a c
we also use to calculate the vertex points The xcoordinates of the
n collocation points may be calculated from or but the
ycoordinates must be calculated from otherwise when using
multivalued or even complex y are resulting The combined plot pl1 de scribing the boundary in x y coordinates and plot pl2 showing the n collocation points is shown by pl3 in Figure
Figure
Cassini boundary with collocation points
Now we calculate the separation constants b and ll the matrix MM satis
fying the boundary condition
value remember the Faber theorem and as the limit of the largest separation
constant b[1] Then we start to calculate the eigenvalue as the root of the determinant of the matrix MM
f Second fk gg
The result k must now be communicated to Mathematica
k=2.99572; Having this eigenvalue which makes the linear system homogeneous for the
unknown amplitudes A we create an inhomogeneous rhs term bbf for the n linear equations rdutn=bbf x the last nth amplitude A
n
and solve
Then we rst check on the correctness of the solution of the linear equations
by inserting and using the symbol Amp for all n unknown amplitudes
The result is satisfying since
is numerically equivalent to zero for the
computer and also the last term due to the chosen A
n
is satisfying
f
g
Now we have to doublecheck the satisfaction of the Cassini boundary con
dition in cartesian coordinates
The rst numerical check is positive
Apparently our solution satis es with an accuracy
the Cassini bound
ary But we make a second graphic check We type
and obtain Figure
Figure
Cassini membrane
Boundary element methods may solve this problem but only nite
element or collocation methods are suitable if the domain is inhomogeneous
partial dierential equation with variable coecients An example is given
circular membrane with radially varying density r
R
r
has
been treated After separation of the membrane equation two ordinary dier
ential equations and result One may then use a collocation
method to satisfy the conditions and the density distribution The following
program does the job
Figure
Circle with collocation points
Next the density distribution will be de ned
Figure depicts this distribution
Figure
Density distribution
Since we do not use the analytic solution Airy functions but solve
and numerically so that the program may be used for quite arbitrary
density distributions we continue
This constitutes a fourfold loop that outputs the following results
fiteration omega Det
g
fiteration omega Det
g
fiteration omega Det
g
fiteration omega Det
g
fiteration omega Det
g
fiteration omega Det
g
To check the change of sign of the determinant a plot of p is made where
p is the iteration parameter p iterat p iterat This is eectuated
by the command
Clear[Detfct];Detfct=Table[{omeg[p],detf[p]},{p,1,3}]; ListPlot[Detfct,PlotJoined->True] resulting in Figure
Figure
Zero of the determinant as function of p
Figure
Improved eigenvalue om
The eigenvalue found om is conform to the Faber theorem and is
entered into the program om=2.3233 If one wants to see the next steps clearly it is recommended to give the
command MatrixForm[M] and to omit the concluding semicolon in the next commands
By the next two commands the solution for the amplitudes Amp is tested by a calculation of the boundary values which includes a test of the solution of
the linear equations for Amp
f g
g Finally the boundary values are tested
Nearly all toroidal problems are not separable Apparently only the Laplace
equation is separable in toroidal coordinates The general solution reads
U
h
AP
m
p
cosh BQ
m
p
cosh
i
C sin pD cos p E sinm F cosm
where the generalized spherical functions P Q satisfy the Legendre equation
d
Q
d
coth
dQ
d
p
m
sinh
Q
Mathematica uses the toroidal coordinate system u v phi a which is ob
tained by rotating bipolar coordinates u v z a about an axis perpendicular to
the axis connecting the two foci of bipolar coordinates The angle parame
terizes the rotation Bipolar coordinates are built around two foci separated
by a Holding u xed produces a family of circles that pass through both
foci Holding v xed produces a family of degenerate ellipses about one of the
foci The coordinate z describes the distance along the common foci default
value a The bipolar coordinates in the z plane are shown in
Figure
Using a trick even the vector Helmholtz equation can be solved
Expressing F
and F
r
by their cartesian coordinates F
F
x
sin
F
y
cos F
r
F
x
cosF
y
sin and deriving both components with respect
to yields
F
F
x
sin
F
y
cos F
r
F
r
F
x
cos
F
y
sin F
Applying the Laplace operator on F
and using F
x
k
F
x
F
y
k
F
y
sin
r
sin cos
r
cos one obtains
F
k
F
r
F
r
F
x
cos
F
y
sin
k
F
F
F
r
F
F
r
k
F
r
r
F
r
r
F
Due to the necessary uniqueness
F r z
F r z one has to assume
a dependence expim so that becomes
r
t
F
k
F
r
F
m
r
F
imF
r
where
t
m
r
is the transversal Laplacian operator Its transversal coordinates q
q
de
scribing the cross section of the torus laying in the cylindrical coordinate
system q
may be r q
may be z may describe an arbitrary torus cross
section Insertion of F
r
from into yields
t
m
r
k
t
m
r
k
F
m
r
F
t
r
m
k
t
r
m
k
F
This structure recalls We thus have two solutions
F
q
q
C
F
C
F
where
t
k
m
r
F
Assuming now q
r q
z and expiz one obtains
d
F
dr
r
dF
dr
k
F
m
r
F
and
F
r z
h
C
Z
m
p
k
r
C
Z
m
p
k
r
i
expiz
Here Z are cylinder functions either Bessel Neumann or Hankel For a
dependence expim one obtains Fock equations
F
r
r z t
A
mk
expi t im
r
k
m
im
rF
r z
r
rk
rF
r z
z
F
z
r z t
A
mk
expi t im
r
k
m
im
rF
r z
z
rk
rF
r z
r
the electric conductivity of the torus wall is high then the boundary condition
B
n
is valid for the magnetic eld normal component Then for an exciting
frequency kc a solution reads
B
r z
X
ml
a
ml
J
m
q
k
l
r
b
ml
N
m
q
k
l
r
cos
l
z cosm
and analogously for B
r
r z B
z
r z Then the unknown partial ampli
tudes a
ml
b
mk
have to be chosen to satisfy div
B Integrating n times
the dierential equations for the eld lines
dr
B
r
dz
B
z
rd
B
or
dz
d
rB
z
B
dr
d
rB
r
B
around the torus n one may nd the cross section of a torus
The totality of all points crossing the z r plane draws a picture of the cross
section in the z r plane Poincar
e map If after n n revolutions of
the eld line around the z axis an earlier point is hit then the toroidal surface
drawn is called resonant toroidal surface
Many toroidal problems have been solved either by collocation methods or
by numerical methods to In plasma physics arbitrary shapes of
the meridional cross section of the torus are of interest One possibility to
solve such problems is to construct the boundary curve after having received
a general solution the second possibility is oered by the collocation method
We rst discuss the possibility of constructing a given shape afterward One
may express the electromagnetic elds
E and
B by a complex quantity
F
E i
B
p
Then Maxwells equations give for cylindrical geometry where
is the angle around the z axis the solution
B
r z
b
J
r c
N
r
p
k
h
b
J
p
k
r
c
N
p
k
r
i
cos kz
B
r
r z k
p
k
h
b
J
p
k
r
c
N
p
k
r
i
sin kz
B
z
r z
b
J
r c
N
r
k
h
b
J
p
k
r
c
N
p
k
r
i
cos kz
where J
p
and N
p
are Bessel and Neumann functions respectively
Now the dierential equations for the eld lines in the r z plane are
dr
dz
r
z
bility We can express B
r
and B
z
by B
This gives
B
r
r
dr
B
r
z
dz dB
r
Thus the lines B
r const D are identical in form with the B
r
B
z
eld lines in the r z plane ie identical with the cross section of the toroidal
resonator To obtain a toroidal resonator of major radiusR and nearly circular
cross section of minor radius a of the form fr z r R
z
a
we have to determine the constants b
c
b
c
D in B
from
z
k
arccos
D b
r
J
r c
rN
r
p
k
b
rJ
p
k
r
c
rN
p
k
r
By choosing for a given k a set of collocation points z
i
zr
i
i
we may generate various arbitrary cross sections
Using similar methods the calculation of exact analytical forcefree three
dimensional toroidal equilibria of arbitrary cross section is possible
To demonstrate that the results are extremely insensitive to the arbitrary
choice of the separation constants
n
we solved the equations for an ax
isymmetric forcefree equilibrium
B
r
z
B
z
r
B
B
r
B
z
B
z
r
r
rB
in the cylindrical coordinate systems
B
N
X
n
cos
n
z
h
A
n
J
p
n
r
B
n
N
p
n
r
i
We assume a toroidal container of major radius R and some eective minor
radius a and de ned by a meridional cross section curve z z
r in the r z
plane This cross section curve may be a circle given by z
p
a
r R
an ellipse or a Cassini curve rR
z
b
rR
z
a
b
The boundary condition along the toroidal surface cross section curve z
r
is rB
const
One obtains for a circle
N
Ra Ra Ra
solved for the radial part of the vector Helmholtz equation for the toroidal
electric eld E
u
u
rr
r
u
r
u
zz
r
u k
u
To satisfy the homogeneous boundary condition ur z
R
on the surface
of an axisymmetric toroid with the meridional cross section curve z z
R
r
we use
and the four particular solutions namely
ur z A
J
kr A
N
kr
h
B
J
p
k
r
B
N
p
k
r
i
cos z
The boundary condition ur z
R
on the cross section curve z z
R
r
yields
z
i
arccos
A
J
kr
i
A
N
kr
i
B
J
p
k
r
i
B
N
p
k
r
i
A workstation delivers the result k A
A
B
B
if an oval de ned numerically by
i
r
i
z
i
is chosen as meridional cross section curve z
i
z
R
r
i
It seems that this
analytical method works quicker than any numerical nite dierences or nite
element method
The collocation method has been applied to many toroidal plasma problems
Thus the electromagnetic eld propagation and eigenfrequency in anisotropic
homogeneous and isotropic inhomogeneous toroidal plasmas of arbitrary cross
section as well as for anisotropic inhomogeneous axisymmetric plasmas has
been investigated In such a plasma the elements of the dielectric tensor are
functions of space
For a plasma magnetized by a toroidal magnetic eld
B
B
e
one has
rr
X
s
Ps
s
r
X
s
Ps
rz
X
i
s
Ps
z
s
Ps
space We have for a toroidal plasma s is the species index
s
s
Rr
s
e
s
B
m
s
where R is the major radius of the torus In toroidal experiments the plasma
density n
is mainly a parabolic function of the distance from the magnetic
axis Expressing n
in cylindrical coordinates r z we have
n
r z n
n
z
n
r R
since
z
r R
Here n
is the maximum density on the magnetic
axis r R z and n
n
a
is determined by the minor toroidal
radius a At z and r R a we have the wall a of the circular
toroidal vessel and the density vanishes In the axisymmetric case the system
of electromagnetic equations splits into three equations for the TE mode with
B
E
r
E
z
and for the TM mode with E
B
r
B
z
For axisymmetric modes we obtain for the TM mode
B
r
i
E
z
B
z
i
r
r
rE
E
z
E
r
r
E
r
r
E
E
Here E
r z is a toroidal electric mode depending on
r z
The solutions have to satisfy the electromagnetic boundary conditions on
the wall of the toroidal vessel On a perfectly conducting wall the tangential
electric component E
t
and the normal magnetic component B
n
must vanish
In our rotational coordinates the form of an arbitrary cross section of the
toroidal vessel is given by a function z zr Projecting into the meridional
plane a vector
A which signi es
E or
B on this curve zr we have
A
t
A
r
cosaA
z
sin a
A
n
A
r
sin aA
z
cosa
where tan a dzdr The electric boundary condition E
t
therefore yields
E
r
E
z
dz
dr
The magnetic boundary condition B
n
yields which represents
the equation for the magnetic eld lines due to B
n
the wall is a special
magnetic eld line and which is a condition for E
Inserting we obtain after
integration
rE
const
This describes the projection of the magnetic eld lines onto the meridional
Ps
we obtain
E
z
E
r
r
E
r
r
E
E
a br cr
cz
E
where
a
m
I
m
E
m
I
n
m
I
n
R
m
E
n
m
E
n
R
b
n
e
R
m
I
m
E
m
I
m
E
R m a m
c
e
n
m
I
m
E
m
I
m
E
n
m
B
kG
Here m
I
and m
E
are the proton ion and electron mass respectively To
solve we make the ansatz E
r z RrUz and obtain
R
rR
r
R
k
a br cr
R
which is a B
ocher equation and
U
cz
U k
U
of the same kind as the Weber equation Here k is the separation constant
which must not be an integer To satisfy the boundary conditions we integrate
and numerically To do this we need a rst arbitrary ap
proximation of the eigenvalue This can be taken from the empty resonator
for a b c Since both equations are homogeneous the eigenvalue
does not depend on the initial conditions We choose therefore two dierent
arbitrary initial conditions R
r R a R
r R a U
U
and
R
r Ra R
r R a U
U
We thus obtain a solution in the
form
E
r z AR
BR
CU
DU
Choosing dierent values k
i
of the separation constant i N we
i
P
N
i
ily given cross section zr by assuming the coordinates r
p
z
p
of p P
points lying on zr we obtain P N linear homogeneous equations
N
X
i
h
A
i
R
i
r
p
B
i
R
i
r
p
i
h
C
i
U
i
z
p
D
i
U
i
z
p
i
The R
i
l
U
i
l
l have been obtained by numerical integration assum
ing an approximate value of the still unknown eigenvalue The system
has then and only then a nontrivial solution for the N unknown
coecients A
i
B
i
C
i
D
i
if the determinant of the coecients of the
system vanishes Since had been assumed arbitrarily the determinant D
will not vanish Calculating D for varying by integrating and
several times and using the regula falsi method a better value of
can be obtained Renewed integration with the better value and repetition
of the procedure nally yields the exact value of As soon as the correct
value of is found the system can be solved for the N unknowns
A
i
B
i
C
i
whereas for example one D
i
will be chosen to be unity out
of i N
We assume a b c or B
T n
m
PE
PI
E
I
Various aspect ratios of circular toroidal vessels and other noncircular cross
sections can be described For i we have r
and z respectively and k we obtain for instance
an ovaloid cross section Now or
cycles and A B C D For k
we obtain and another oval Many other forms have been
produced by assuming various r
i
z
i
The collocation method can also be used for the lowfrequency MHD waves
Starting from the equations of continuity of motion Ohms law and Maxwells
equations we used the usual linearization Elimination of B
r
B
z
B
v
and
a time dependence expi t yields two ordinary dierential equations of
rst order with variable coecients
i
dv
r
dr
c
A
nc
s
kv
z
c
A
nc
s
k
m
r
k
c
A
iv
r
r
nc
s
r
d
dr
c
A
dv
z
dr
r
m
c
A
i
dv
r
dr
krc
A
v
z
r
d
dr
rk
c
A
m
c
A
r
s
mode number in the direction and the Alfv
en speed is given by c
A
B
r
r so that n and c
A
depend on r Due to the gyrotropy of a
magnetized plasma the phase factor expi i appears so that we make
a snapshot at time t We then have two possibilities to make the
dierential equations real modes of type replace v
r
iv
r
and modes of
type replace v
z
iv
z
B
v
are the equilibrium quantities c
s
is
the adiabatic sonic speed Taking circular cylinder coordinates we assumed
r
B
r B
re
B
r RB
Rr R is the major and a the
minor torus radius Furthermore a dependence expim ikz has been
assumed for all wave quantities
If the meridional cross section in the r z plane of the containing toroidal
surface is described by a curve z
zr and if tan dz
dr then v
n
v
r
sin v
z
cos Inserting the real solutions v
z
and v
r
in the form for
mode type we get
v
z
r z
KM
X
kms
A
s
km
v
km
zs
r cos kz cosm
v
r
r z
KM
X
kms
A
s
km
v
km
zs
r sin kz cosm
where s indicates two dierent but arbitrary initial conditions used in
the numerical integration of One obtains for P collocation
points r
i
z
i
i P a system of P equations for the unknown partial am
plitudes A
s
km
Since the sum over s gives two values and if summation
over m gives M values and over k one might have K values the number N of
unknown partial amplitudes A
s
km
is given by P M K The number
P of collocation points determines the accuracy of the solution For xed
P it has been shown earlier that the accuracy of the results is practically
independent of the arbitrary choice of the separation constants k Accord
ing to our experience with axisymmetric MHD modes we make the choice
k a a a Ka where a is the eective minor torus radius for ar
bitrary cross section When the kvalues are given then the condition v
n
with v
r
v
z
inserted constitutes a system of P linear homogeneous equations
for the M K unknowns A
s
km
To solve this system the determinant
D of the coecients known at r
i
z
i
must vanish This condition determines
the global eigenfrequency In order to be able to integrate
the density distribution of the plasma and the frequency have to be known
We thus integrate the dierential equations which deliver the elements of
the determinant with an initial guess for and calculate D Looking for
a root of the function D we obtain an improved value for and a new
global eigenfrequency is exhibited best by the distribution
x A expx
x r Ra
As soon as the solutions are known it is possible to cal
culate the electric eld E
r
r z and E
z
r z When they are known the
dierential equations for the electric eld lines dzE
z
r z drE
r
r z have
to be integrated numerically The results obtained are summarized Here the
following parameters have been chosen major radius R m B
Tesla
T keV c
s
eTm a m
Curve b m Frequency type Frequency type
Circle a E Hz E Hz
Circle a E Hz E Hz
Circle a E Hz E Hz
Circle a E Hz E Hz
Ellipse E Hz E Hz
Ellipse E Hz
We see that for mode type the frequency increases but for mode type
decreases with increasing inhomogeneity With increasing mode number m
the frequency decreases for type but increases for type Probably for ro
tating patterns the frequencies converge Also the dependence of the global
eigenfrequency on the aspect ratio has been investigated As it is expected
for large aspect ratio A one obtains the solution for a cylindrical plasma
Problems
Modify the program Cassmem to describe a circular membrane of radius
R and calculate the rst four eigenvalues Compare the results with
the values of the roots of the Bessel functions given earlier
Demonstrate that the program VARMEM yields om for
constant surface density of the membrane
Modify VARMEM for an elliptic membrane of semiaxes a b
and constant density Is it possible to reproduce the eigenvalues of the
function of x y and z Solve nx y z B
x y zux y z for
nboundary for a cube and a sphere
Use the following collocation program Homcircplate to calculate the
lowest eigenvalue of a circular plate of radius R in cartesian coor
dinates Read in section equations and remember
the result Is the collocation program able to reproduce
this value Execute the plot commands Try to vary the separation
constants within the interval b
n
fst
Does the Faber theorem hold for plates too Yes
Compare the k above with
Holes in the domain Two boundaries belonging to
dierent coordinate systems
The theorem that the solution of an elliptic partial dierential equation is
uniquely determined by ONE closed boundary is valid for analytic solutions
only A partial dierential equation of second order has however two arbi
trary functions in its general solution The Laplace equation has two so
lutions and the Helmholtz equation has two solutions J
p
and N
p
but we excluded the latter due to its singularity at r Now let us
investigate the combination of the two particular solutions
We rst consider the Laplace equation in polar coordinates r
Assuming the inhomogeneous boundary conditions on a circular ring of radii
R
and R
the conditions read
Ur
rR
f
Ur R
f
The solution is then given by
Ur a
a
R
ln
r
R
X
k
a
k
R
k
kR
k
a
k
r
k
kR
k
a
k
R
k
a
k
r
k
k
R
k
R
k
R
k
R
k
cos k
X
b
k
R
k
kR
k
b
k
r
k
kR
k
b
k
R
k
b
k
r
k
k
k
k
k
sin k
a
Z
f
d a
k
Z
f
cos kd
a
k
Z
f
cos kd a
Z
f
d
b
k
Z
f
sin kd b
k
Z
f
sin kd
One may remember the Laplacian singularity ln r which we know from
equation
A sector of a circular ring with the radii R
R
R
and the central angle
and subjected to the inhomogeneous boundary conditions
UR
UR
U
Ur Ur
has for the Laplace equation the solution
Ur
U
X
n
rR
n
R
r
n
R
R
n
R
R
n
sin
n
n
n
n
The solutions of the Helmholtz equation are of greater interest
We consider the homogeneous boundary problem of a circular ring membrane
with the radii r R
R
R
There are now two boundary conditions
uR
uR
In plane polar coordinates r Mathematica has given the solution
which we shall use now The Neumann function Y is the second solution of
the Bessel equation It may be de ned by
Y
x
J
x cos J
x
sin
where is not an integer For p integer the Hospital rule yields
Y
p
x
ln
x
J
p
For u p The result may also be derived from
The two boundary conditions demand
AJ
p
c
R
BY
p
R
c
AJ
p
R
BY
p
R
be solved if the determinant vanishes
J
p
R
c
Y
p
R
c
Y
p
c
R
J
p
c
R
For given R
R
and c Mathematica gives the solution for p by
This gives x
For a sector of a circular ring membrane the solution analogous to
is found by replacing the powers of r by Bessel functions Thus an annulus
a r b with a sector a has solutions of the form
J
m
j
mn
r
a
sin
m
Bessel functions of fractional order appear and the order is determined by
the angle
Now we understand the role of singularities they cut out a hole within the
domain circumscribed by a closed boundary The method worked quite well
in plane polar coordinate systems r which has a natural singularity at
r Although we had some success with the solution function cosh of the
Laplace equation in Figure and equation we are sceptic since
the function cosh tends to assume large values Therefore we still use polar
coordinates Let Mathematica do the work
<<Calculus!VectorAnalysis! SetCoordinates[Cylindrical] which gives CylindricalRr Ttheta Zz
Now we make a setup and solve the Laplace equation
u[Rr_,Ttheta_,0]:=R[Rr]*Cos[m*Ttheta] LL=Laplacian[u[Rr,Ttheta,0]] Collect[LL,Cos[m*Ttheta]] OK the result looks nice We obtain an ordinary dierential equation for
the rdependence
Cosm Ttheta
m
RRr
Rr
R
Rr Rr R
RrRr
and ask Mathematica to solve it
DSolve[R[r]+R[r]/r-mˆ2*R[r]/rˆ2==0, R[r],r] which yields a not very nice expression Rr C Coshm Logr i C Sinhm Logr
So we would like to test our own solution
a c Logr
nn
X
n
cn r
n
bn r
n
Cosn
To verify that this is really a solution of the Laplace equation in polar co
ordinates r we type Simplify[D[U[r,],{r,2}]+D[U[r,],r]/r+ D[U[r,],{,2}]/rˆ2] which is apparently too dicult yielding
P
nn
n
r
r
cn r
n
bn r
n
Cosn
r
nn
X
n
r
r
cn r
n
bn r
n
Cosn
P
nn
n
r
r
cn r
n
bn r
n
Cosn
r
OK then we simplify and set
Clear[U,a,b]; U[r,]=ao+c0*Log[r]+(c[n]*rˆn+b[n]*rˆ(-n))*Cos[n*] and the next command veri es our solution
Simplify[D[U[r,],{r,2}]+D[U[r,],r]/r+ D[U[r,],{,2}]/rˆ2] yielding Now we are able to solve the Laplace equation with two boundaries
belonging to two dierent coordinate systems As the outer boundary we
choose a circle of radius R with the homogeneous boundary condi
tion Ur R For the inner boundary we choose
a square of lateral length with the inhomogeneous boundary condition
Usquare U
const x y The homogeneous bound
ary condition on the circle demands
a
c
lnR and c
n
R
n
b
n
R
n
so that the solution becomes
Ur c
ln
r
R
X
n
c
n
r
n
R
n
r
n
cosn
This solves the problem of the electromagnetic potential between an in nite
square prism and a circular cylinder
Since the c
c
n
are arbitrary they can be used to satisfy the inner in
homogeneous boundary condition on the square We do this using a col
location method Avoiding the corner points we choose the collocation
points x
i
i y
i
and r
i
p
x
i
y
i
i
arctany
i
x
i
Then the inhomogeneous inner boundary con
dition reads
Ur
i
i
c
ln
r
i
R
X
n
c
n
r
n
i
R
n
r
n
i
cosn
i
U
i
the unknowns c
c
c
c
c
c
Now let Mathematica do the work
Figure
Capacitor consisting of square and circle
the coordinates of the collocation points
c Log ri c Cosi
ri
ri
c Cos i
ri
ri
c Cos i
ri
ri
c Cos i
ri
ri
c Cos i
ri
ri
NowMathematica calculates the matrix and solves the linear system for the
unknowns c c
R
Cosi
R
ri
ri
Cos i
R
ri
ri
Cos i
R
ri
ri
Cos i
R
ri
ri
Cos i
R
ri
ri
g fi g
R=1.;U0=10.; MatrixForm[M]; b={U0,U0,U0,U0,U0,U0}; B=LinearSolve[M,b]; We now help Mathematica By copy and paste of the results we inform the
program on the results obtained for c c etc
Sum[c[n]*(rˆn-Rˆ(2*n)*rˆ(-n))*Cos[n*],{n,5}]; Table[V[r[i],[i]],{i,1,6}] Table[V[R,[i]],{i,1,6}] The result of checking the satisfaction of both boundary conditions is delightful
f g
f g
Thus we have demonstrated that collocation methods are able to solve an
elliptic partial dierential equation even for two closed boundaries even in
the case that one boundary value problem is homogeneous and the other one
inhomogeneous
Now we exchange the boundaries the outer inhomogeneous boundary is
described by a rectangle and the inner homogeneous boundary is given
by a circle of radius Then the boundary conditions read
Ur Ux y cosx Ux y
We use very detailled commands and choose collocation points and identify
them using a symbol
point
symbol q l t u v w p
r
i
and r
i
p
x
i
y
i
i
arctany
i
x
i
l
n
n
p
n
p
cosn arctan
For ln r
i
the abbreviation a
i
ln r
i
will be used Mathematica does the work
At rst one has to dene a matrix
Then the boundary conditions create the rhs term of the linear system
b=N[{0, Cos[Pi/8],Cos[Pi/4],Cos[3*Pi/8],0,0, 0}] Then we solve the system by the command
n
Informing Mathematica about these values and the denitions
Figure
Top view of
Figure
Threedimensional view of
dierent coordinate systems can be solved by other methods than by colloca
tion methods A method of solving these nonuniform boundary problems has
been known since but it is cumbersome Tinhofer has applied
this method on our problem of a rectangle a b with a circular hole of radius
c For the Laplace equation the following boundary conditions have been
assumed
Ux a y Ux y b
Ur c
P
n
D
n
cos n
Now the cartesian solution will be expressed by the polar solutions
Ur
a
X
c
r
a
cos b
sin
To do this the formulae x r cos y r sin
cosx
e
ix
e
ix
e
ir cos
e
ir cos
coshy
e
y
e
y
e
r sin
e
r sin
cos
n x
a
cosh
n y
a
Re
exp
i
n
r
a
e
i
exp
i
n
r
a
e
i
are helpful
We see that the function cos as well as cosh are both expressed by their com
mon mother function exp recall the grandmother hypergeometric equation
Furthermore one has to use
e
x
X
s
x
s
s
e
ire
i
X
s
i
s
s
r
s
e
is
s
and s s i
s
s
so that from results
Re
P
s
s
i
s
s
n
r
a
s
cos s i sin s
P
s
s
i
s
s
n
r
a
s
cos s i sin s
P
s
n
r
s
cos s
cos
n y
b
cosh
n x
b
P
s
n
r
b
s
s
s"
cos s
Now it is possible to rewrite the solution in the form
Ux y Ur ln
s
x
y
a
b
X
n
A
n
X
s
n r
a
s
s
s"
cos s
X
n
B
n
X
s
n r
b
s
s
s"
cos s
where A
n
and B
n
are known from Introducing the abbreviations
C
s
P
n
A
n
n
a
s
s
s"
P
n
B
n
n
b
s
s
s"
the solution can be rewritten using
s
for s but otherwise
as
Ux y Ur
X
s
ln
r
p
a
b
s
C
s
r
s
cos s
To satisfy the boundary condition on the circle one uses the fact that the
derivative of a solution of a dierential equation is again a solution of the
dierential equation We build new solutions
U
x y
x
Ux y
ln r
x
x y
a b
X
s
C
s
r
s
cos s
where
C
s
P
n
A
n
n
a
s
s
s"
B
n
n
s
s
A
n
a cosh
n b
a
a
Z
a
x
ln r
x y
a b
cos
n
x
a
dx
B
n
b cosh
n a
b
b
Z
b
x
ln r
x y
a b
cos
n
y
b
dy
Using
ln r
x
"r
cos
the new partial solutions assume the form
U
x y U
r
"
p
a
b
cos
"
r
cos
X
s
C
s
r
cos
Using
#
r
X
s
ln
r
p
a
b
s
C
s
r
s
cos s
#
r
"
p
a
b
cos
C
X
s
"
r
s
C
s
r
s
cos s
one can write
Ur
X
S
#
r
Here the expansion coecients S
are still unknown they depend on the
boundary condition which reads now
P
s
D
s
cos s The coe
cients S
have to be calculated from
X
s
S
D
s
ln
c
p
a
b
C
"
p
a
b
cos
C
s
C
s
c
s
s
h
"c
s
C
s
c
s
i
Some of these nonuniform problems may be very important In the nuclear
power station KAHL it has been crucial that the neutronabsorbing control
rods could control the chain reaction In this nuclear reactor the control rods
were arranged on the surface of a cone The boundary value problem of the
neutron diusion equation could not be solved numerically with the necessary
accuracy since two boundaries on cylinder and cone had to be satis ed and
the vertex of the cone presented a singularity
Membrane domains could exhibit holes too We rst investigate if collo
cation methods can be avoided A clamped square membrane
a
x
a
b
a
y
a
can again be described by a modi ed Helmholtz equation
of type
uE R ru u
a
y
u
x
a
where is a given parameter This cuts out a circular hole of radius R
Perturbation theory
u u
u
u
E E
E
E
with
u
a
cos
x
a
cos
y
a
u
X
c
nm
u
nm
u
nm
a
cos
n
a
x cos
m
a
y
E
E
a
Z Z
R r cos
x
a
cos
y
a
dxdy
E
R
a
J
R
a
J
p
R
a
might be the answer E Rietsch For R a one obtains
E
a
a
One might also think to expand the point forces represented by R r with
respect to
cos
x cos
y
variable density Many problems concerning two boundaries need
to know how to handle corners in a boundary Problems will thus be oered
in the next section Anyway try to solve this problem The program DUM
solves an inner homogeneous boundary on a square and an outer homogeneous
bondary on a circle Modify the program to calculate an inner homogeneous
boundary on the square and an outer inhomogeneous boundary on the circle
Solution the modi cations have to be made in steps and In step the
resulting two output lines should then be exchanged numerically
Problems
Solve the inhomogeneous problem of a rectangular membrane with
an inner homogeneous circular boundary of radius The inhomoge
neous values on the rectangle should be
p1=ListPlot[Table[{x[l],y[l]},{l,1,18}], PlotStyle->PointSize[1/40],AspectRatio->0.5, PlotRange->{{-2.,2.},{0.,2.}}]; Here N guarantees that one gets numerical values from the function
arctan The command ListPlot generates Figure showing the lo cation of all collocation points PointSize determines the diameter of the points by of the dimension of the plot AspectRatio->0.5 xes the ratio of the y to the x dimension and PlotRange describes the limits of the plot
Figure
List plot of collocation points
Now one has to de ne the matrix m and the two boundary values
It might be that you get a message concerning a badly conditioned
matrix In this case remember problem of section
BC=LinearSolve[RM,RP] A=Table[BC[[lk]],[lk,1,n}] The command Rationalize[M,p] transforms all elements of the ma trix M into rational numbers fractions using a tolerance of p in the approximation
g Observe the condition
/;xˆ2+yˆ2>=1 after the de nition of the function g[x_,y_] which should de ne the function gx y for values only outside and on the
circular inner boundary
The following commands should generate Figures
Figure
Contourplot of rectangle with circular hole
Figure
Rectangular membrane with circular hole
Figure
Rectangular membrane from default viewpoint
Figure
Membrane sideways
Try to write a program to calculate the eigenvalue of a circular mem
brane with a rectangular hole Use polar or cartesian coordinates So
lution very dicult approximation very rough Collocation methods
often fail to solve eigenvalue problems with two dierent boundaries
Corners in the boundary
Corners in the boundary of a domain present new mathematical diculties
In recent years quite a great number of papers have been published on this
subject Due to the large amount of material we only may quote a small
selection A collocation method using boundary points has been discussed for
a triangular plate as early as and the same author treated vibrations
of a conical bar Triangular and rhombic plates were also discussed
Vekua drew the attention to the fact that smoother boundaries yield
better approximations Fox stated that it is necessary to take the origin
of the polar coordinate system as the vertex of a corner with the angle
and that a function
ur
N
X
n
C
n
J
n
p
r sinn
where n
n
n N satis es the boundary condition on both sides of
the angle and has the correct singularity at the corner For an Lshaped
membrane one has an angle and so that Bessel functions of
fractional order appear Such reentrant corners also have been investigated
tion methods and others but state that often engineering applications solved
by the nite element method are of limited accuracy A collocation
method applied on a nonconvex domain with one boundary generally yields
poor accuracy too It is therefore interesting to know the computable
bounds for eigenvalues of elliptic operators Another review on prob
lems on corner domains has been given by Dauge Finally Ukrainian
authors draw the attention to the obstacles hindering the use of variational
methods as well as of nite dierence nite elements and boundary elements
for such boundary problems and present the theory of R functions In
section Figure showed a rectangle with an incision The angle of
that incision is apparently given by reentrant corner We rst
consider the laplace equation in cartesian coordinates with homogeneous
boundary conditions
u y y
ux x
We assume diagonal symmetry around the diagonal straight line y x For
the remaining boundaries we assume inhomogeneous conditions
u y sin
y y
ux x
These conditions guarantee continuity in the points and Avoiding
cartesian solutions with the disagreable cosh we use polar coordinates To
satisfy the boundary conditions along the lines x y thus
and along x y thus r we set up
ux y ur
N
X
n
a
n
r
n
sin
n
This is an exact solution of the Laplace equation in polar coordinates com
pare similar solutions like
We now use the partial amplitudes a
n
to satisfy the boundary conditions
For this purpose several methods can be used We prefer the collocation
method choosing collocation points
Using the abbreviations
r
i
q
x
i
y
i
i
arctan y
i
x
i
fx
i
y
i
s
in
sin n
i
u x
i
y
i
X
a
n
r
n
i
s
in
partial amplitudes a
n
and plot ux y This can be done by a Mathematica
program You should Input x[i],y[i],r,f,m b={N[Sin[Pi*y[1]/2]],N[Sin[Pi*y[2]/2]], N[Sin[Pi*y[3]/2]],N[Sin[Pi*y[4]/2]],1,1,1,1} f g
g
Figure
Reentrant corner
To solve a similar boundary value problem for the Helmholtz equation
we use the corresponding solutions
ur
N
X
A
n
J
n
kr sin
n
b={Sin[Pi*y[1]/2],Sin[Pi*y[2]/2], Sin[Pi*y[3]/2],Sin[Pi*y[4]/2],Sin[Pi*y[5]/2],1,1,1,1} Far more complicated are problems with two rectangular boundaries For
the Laplace equation Still private communication has solved such a prob
lem For u he assumes a homogeneous boundary problem for the inner
square and an inhomogeneous problem u on the outer square Due to the
symmetry only half of the rst quadrant will be considered A trial function
is taken in the form
ur a
a
ln r
n
X
k
a
k
r
k
a
k
r
k
cosk
A collocation method with N n points showed heavy oscillations between
the collocation points The method obviously suered from the corners with
an angle where is not dierentiable In fact it allows a singular
expansion around this point r
by
ur
r
sin
Therefore the trial function had been modi ed by adding the term
n
X
k
a
k
r
a
k
r
k
sink
a
n
r
sin
This trial function gave slightly better results Furthermore Still also con
sidered conformal mapping
An even more general problem has also been considered by Still The
boundary of a triangle may consist of two straight lines
and a curve described by
gr The Laplace equation solution should
satisfy the boundary values g
r on
for One can assume
g
r
N
X
n
a
n
r
n
for g
r
N
X
n
a
n
r
n
for
or simpler
g
r a
r R
for g
r a
R
for
g
r a
R
sin
Br
sin
for
R
r R
Then the singularity free setup may be suggest
ing itself
ur u
r u
r a
N
X
n
a
n
cosn b
n
sinn r
n
N
X
n
n
r
n
sinn
where If the boundary conditions would contain a discontinuity
the addition of singular terms like c cr coslog r sin which solve the
Laplace equation could help The
n
in have to be determined by
the boundary condition on $
g
r u
r
N
X
n
n
r
n
sinn
whilst for $
$
one has
u
u
r
If instead of the Laplace the Helmholtz equation has to be solved one has
to make the replacement
r
n
cosn J
n
kr cosn
To execute numerical calculations we choose R
R
a
a
B Satisfying the boundary conditions results in
g
r r ur a
N
X
n
a
n
r
n
g
r ur
a
N
X
n
a
n
cosn b
n
sinn r
n
g
r sin
Br
sin
a
N
X
n
a
n
cosn b
n
sinn r
n
N
X
n
n
r
n
sinn
The equations determine the coecients a
a
n
and b
n
in
$
and $
Equation may yield the
n
This may be done using a
collocation method First gives a
a
a
n
n
From one obtains for
b
cos
b
n
sin
i
Br
i
sin
i
cos
i
N
X
n
n
r
n
i
sin
n
i
i
i N
Now choosing and four collocation points on $
determining
one obtains Table
Table The expansion coecients
i x
i
y
i
b
i
I
By similar methods and ngers crossed corner problems of this type can
be solved These might include rectangular domains with rectangular or trian
gular or hexagonal holes or domains without a hole but exhibiting polygonal
boundaries
Problems
Find the lowest eigenvalue of an hexagonal membrane with no holes and
homogeneous boundary conditions Assume an hexagon with the same
area as a circular membrane of radius so that the Faber theorem
may help to control your results You may choose collocation points
and a radius of the circumcircle equal to
q
p
guaranteeing the
same area of the circular membrane Then you should obtain the
eigenvalue of the circular membrane Try to nd
the eigenvalue by plotting Detm as function of eigenvalue by narrowing
the search interval
Discuss a rectangular membrane with a rectangular hole Boundary
conditions might be value on the outer inner rectangle and ho
mogeneous clamped membrane on the inner outer rectangle Use
polar coordinates Make a PlotD of the solution