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