MEMOIRS OF
NUMERICAL
MATHEMATICSNUMBER 6

1979

ED rrOH I AL BOARD

HIROSHI FUJITA (Chairman)
Department of Mathematics
Faculty of Science
University of Tokyo
SHIGERU FURUYA
Department of Mathematics
Faculty of Science and Engineering
Aoyama-gakuln University
SIN HITOTUMATU (Executive Editor)
Research Institute forMathematical Sciences
Kyoto University
MASAYA YAMAGUTI
Department of Mathematics
Faculty of Science
Kyoto University

The MEMOIRS OF NUMERICAL MATHEMATICS will publish
original papers concerned with numerical mathematics, andhenceforth will be publishpd once a year.
All correspo~dence should be addressed to
Sin Hitotumatu
Research Institute for Mathematical Sciences
Kyoto University, Sakyo-ku, Kyoto
606, Japan

CONTENTS

Fumio KIKUCHI

On a Mixed Finite Element Scheme for
Buckling Analysis of Plates
1

Taketomo MITSUIThe Initial-Value Adjusting Method
for Solving Problems of the Least
Squares Type of Ordinary Differential
Equations21

Vitoriano RUAS B.SANTOS

A Direct Method for Solving Twodimensional One Phase Stefan Problems

Teruo USHIJIMA

Error Estimates for the Lumped Mass
Approximation of the Heat Equation

39

65

On a Mixed Finite Element Scheme for Buckling Analysis of Plates
Fumio KIKUCHI *1. Introduction
Our problem is expressed by the buckling equation of clamped
plates'\ A Zu+~u

2 T .. a 2u/ax.ax. = 0
Li,j=l 1J

1

in

J

n,

(1)

on an,

u = au/an = 0

where n c: RZ is a ~on-obtuse polygonal domain with boundary an,

6.

the Laplace operator, a/an differentiation in the outward normal
direction of an, x = (x 1 ,x Z) the independent variable of RZ, andTij (1 ~ i,j < Z) are sufficiently smooth given functions of x suchthat2

r aT. ./ax.J = 0
j =1 1J(Z)(i= 1,Z).

We are interested· in a non-trivial function u = u(x) and a (real)
number A that satisfy (1).
In our mixed method, we first decompose (1) into a system

-

6.u

=v

- A 6.v +2' 2
rT .. a u/ax.ax. = 0
. . 1
J

1, J

( 3)

=

1J

u = au/an = 0

1

in

n,

on an,* Institute of Space and Aeronautical Science, University of Tokyo,
Tokyo, 153, Japan.

1and then use appropriate finite element spaces and weak forms fordiscretization.

This type of mixed method has been proposed bymany investigators, see for example Kikuchi [1], and Ciarlet andGlowinski [2].
As for the error analysis of the mixed method for the above
eigenvalue problem, we can refer to Ishihara [3], where error estimates of approximate eigenpairs are derived for the piecewise
linear finite element model considered over nearly uniform triangulations.

The analysis is based on Miyoshi's results for the mixed

approximation of the corresponding boundary value problems [4J.
It also uses the Rayleigh and the min-max principles, and hence its
validity is restricted to equations such as A~2u + ~u = O.
The main objective of this note is to present error analysisof the above mixed method applied to more general cases, where theRayleigh and the min-max principles are not available, or the employed meshes are not necessarily uniform.

We will use the so-called

spectral projection to our end, see for example Grigorieff [5].

We

will also utilize the results obtained by Scholz [6J to obtain order
estimates of errors of approximate eigenpairs.

In our analysis,

it is essential to prove a certain compactness property of a familyof approximate operators.

The techniques developed in this note

wilL be available for higher order finite elements, and also forthe different type of mixed models where u itself as well as its
all second order derivatives are considered as independent functions
to be approximated (se for example Miyoshi [4]).
In this note, C will be used as a generic positive constant,which ~ay take different values when it appears in different places.

22. Preliminaries
We will use the Sobolev spaces Hm(O) and H~(n) equipped withthe same normII

11m (m = 1,2, ... }.

The inner product and the normof L2 (n) wi 11 be denoted by ( , ) andII II, respectively.

The problem (1) expressed in a weak form is to find a pair
{A,U} e Rl x H~(n) such that u ; 0 and

for all(4)2
u e HO (n)

,

where
2

(5)

b(u,u).. au/ax.1 ,au/ax.)= . r. 1 Cr 1)
J
1,)=

.

Here we have used the relation (2).
Another expression may be given to the same problem, if a weakform of (3) is employed: find a triplet {A, u, v} e Rl x H~ (n) x HI (n)such that u ; 0 and(6)(Vu, Vv) = (v, v)(

V-Iv e H (n) ) ,A(VV,VU) '"'b(u,u)

where (Vu,Vv) for example implies
2(7)

eVu,Vv)==r (au/ax. ,av/axi)
. 1
1=
1

In order to check the equivalence between (4) and (6), notice th3:t thesolution u of (4) is necessarily belongs to H4 (n) 'since n is a nonobtuse polygonal domain, see Mizutani [7].

Then the conditionv = - ~u e: Hlen) in (6) can be easily justified.

33. Finite Element Scheme
We triangulate n in the usual way.

We assume that the consi-

dered triangulation Th for h > 0 is K-regular in the following
sense, see Scholz [6]:(i) T e Th is a (closed) triangle, each side of which is either
a portion of an or a side of an adjacent triangle in T; ,(ii) there is a fixed constant K > 0 such that two circles Kl
and KZ exist with the properties
radius of Kl

radius of KZ = Kh

= h/K

for allandT e Th.

In the sequel, we will only consider a K-regular family of triangulations {Th} with a fixed K > 0, and the case when h ~o.

Let Xh = Xh(Th) be the space of continuous functions which arelinear polynomials in each T e: Th.
al subspace of Hl(n).

Then Xh is a finite dimension-

We also use the space x~ of all functions in

Xh that vanish on an, that is,x~ = xhA H~ (n)'.

Using the above two finite element spaces Xh and XhO' we canintroduce a mixed finite element approximation to (6): find a triplet
1{Ah,uh,v h } E: Rhhx Xo x X . such that u h

r 0 and(8)

For each u h E x~"

we can 'define . ~h u h E

4Xh by the relation

for all(9)

Then we can rewrite the system of equations (8) into a single onefor {Ah,u h }E: Rl x x~

(10)

for allXh 'u- h e O~huh.

where vh has been eliminated by the relation v h = 4. Properties of ~h2

In the space HO (n),II lm II can be regarded as a norm equiva-

II u 11 2 , and hence we have the following properties: Let{ui}~=l be a sequence in H~(n) such that II 6U i II ~ C for a certain

lent topositive constant C.

Then we can choose a sub-sequence, again

denoted by {ui}~=l for simplicity, such that for i ~ ~,

u.

(11)

1

~

U

2

1

weakly in HO(n), and 5trongly in HO(n),where u is a certain element in H~(n).

We will show that theoperator ~h' when considered over an appropriate family of X~,
possesses properties in a sense corresponding to the above.

Lemma 1

There exists a positive constant C = C(O) such that

(12)

where

II vU h II

5

Proof

The former part of (12) is nothing but the Poincare inequ-a1ity, while the latter is obtained by substituting v h = u h 6 xoh exh
and using the Schwarz inequality in (9). lUI

Lemma 2Let {Th(i)}~1.=1 be a (K-regu1ar) sequence of triangulations

such that lim h(i) = O.
isuch that

(13)

uh(i)

h (i)

e Xo

Consider a sequence of functions {uh(i)}~=l'

where C is a positive constant, and x~(i) and 6 h (i) are the finite
element space x~ and the operator 6 h associated with Th(i), respective1y.

00

Then we can choose a subsequence, again denoted by {u h (i)}i=l

for simplicity, such that for i~00

weakly in L2 (n) ,(14)strongly in' H~(n)where u is a certain e1emen~ of H~(n).
ProofIn this proof, we omit the index i, and h ~ 0 implies i

~

00.

o

1

Since 6 h u h is uniformly bounded in L2 (n), we can show the exi-

stence of a subsequence, again denoted by {uh }, such that for h ~ 0
weakly in L2 (fl)weakly in H~(fl) and strongly in L2 (fl) ,6

1
where wand u are certain elements of L2 (D) and HO(D),respectively,

and we have used the results of Lemma 1.
find for each v

E

On the other hand, we can

e HI (0) a sequence {vh } such that'1strongly in H
(D) as h ~ O.

xh, and

The existence of such an approximate sequence is assured by theusual approximation theory, see for example Ciarlet [81.

Takingthe limit of

we find(Vu, Vv)

(*)

== -

(w, v)

for allv e HI (0) .

Then we have for h ~ 0

1from which follows the strong convergence of u h to u in HO(D).
o2Choosing v from H~(D) in (*), we have

(Vu,Vv)

==

-

(w,v)

for all

Since D is a non-obtuse polygonal domain, we can assure that1·

2

u e HO (D) (\ H (n) and l1u = w, see Kondrat' ev [9].
Green's formula to the left-hand side of (*) yields-(l1u, v) +J~u v dy = - (w, v)

an n

7

Applying thewhere dy is the line element on an.
2• 0 on an, and hence u ~ HO(D).

We can now conclude that au/an/#/

5. Fundamental Properties of the Approximate Eigenvalue Problems
As before, we assume that the considered family of triangulations is K-regular.

Define a, a h , and aCE) by

a = set of all eigenvalues for the problem (4),0h = set of all eigenvalues for the problem (10),
0(£)

= U

]A-£,A+£[,

AEO

where £ is an arbitrary positive number.

Clearly a is a bounded

countable subset of Rl without any accumulation points except zero.
In particular, zero is either an eigenvalue of (4) or an accumulation point of o.

These properties directly follow from the theoryof compact operators.

We can also show that 0h is a finite bounded

subset of Rl, the boundedness being uniform with respect to h, as
follows from Lemma 1.

Notice that 0(£) is an open set.

In the sequel, we will only consider a non-ae~o eigenvalue Aoeo, which can be fixed arbitrarily.

The constants to appear in theforthcoming lemmas and theorems may differ with Ao ' but we will makeno special distinction among them.

DefineE = eigenspace for AO
m = dimension of E; 1 ~ m <00,(~$.,~$.)

a basis of E such that1

where 0 .. is the Kronecker delta.
1J

8

J

• 0 .. (1 ~ i,j < m) ,1J-

Let K bea positive constant such that

DefineEh = linear hull of eigenfunctions corresponding to all eigenvalues in <1h {\mh = dimension of Eh[A O - K,A O + K]
0 ~ mh < 00

.
mh
h
a
basis
of
E
such that each 'f'~h1. is an eigenfunction'l'hi i=l -{ ~}of (10) and satisfies (6 h $hi,6 h $hj)

= 0ij (i,j ~ mh ) ,Ahi = eigenvalue corresponding to $hi (i ~ mh ) .
The following lemmas may be derived from the results establishedin the preceding section and by reduction to absurdity.

Lemma 3For each E > 0, we can find a positive number hO= hO(E)such thatfor any(15)

That is, any element of <1h is close to a certain element of <1 when his sufficiently small.
Proof

Fix EO > 0, and consider the case E ~ EO'

ists a positive c~nstant L such that
<1h C[-L,L]aCE) C [-L,L]

9

Then there ex-Assume the contrary to (15).
{{ Ah (i) , uh (ii}

Then we can find a sequence:""1 such that

lim h(i)
i+oo

e

=0
.
.
f
.
Xh0 (i) 1Sof (10) corresponding toan e1genunct10nand satisfies Ilflh(i)uh(i) II=1

where {Xh(i)}~is a sequence of finite element spaces associatedo1=1with an appropriate sequence of triangulations {Th(i)}~=I' Hereafter,
we omit the index i, and h ~ 0 means i ~~.

Since IA h ' + II flhU h II is

uniformly bounded, we can apply Lemma 2 to show the existence of asubsequence, again denoted by HAh,uhH, such that for h ~ 0
Ah + A e Rl
uh + ustrongly in H~(O)flhU h ~ flu

weakly in L2 (0)2where U is an element of HO(Q).

,

Notice that A helongs to a closed

set [-L,L] \ 0(£), and in particular A fila O(E).

On the other hand,

we can find for each U E: H~(Q) an approximate sequence {uh } such-hthat u h E: Xo and for h ~ 0strongly in H~(n)strongly in L2 (0)10,Here we have used the results of Scholz [6] to assure the existenceof a sequence strongly convergent in the above sense.

Taking thelimit ofwe haveA(6u,6U)= b(u,u)for any

Since A -1= 0, the above relation implies u = O.

we have A = b (u, u) = o.
mulation point of o.
a contradiction.

Lemma 4Since A ttI=0,

u E H~ (n)Taking the limit ofthis means A = 0 is an accu-

But A also belongs to .[-L,L] \ o(€), and we havefifi

We can find a positive constant hO such that(16)

for anyh ~ hO .

That is, the dimension of Eh never exceeds that of E when h is sufficiently small.
Proof

We only sketch the proof.

Assume the contrary to (16).

Then we can find a sequence of triangulations (Th(i)]~=l such that

lim h(i)
i-+oo

=0

,

mh(i) ~ m + 1

11

where mh(i) is the dimension of Eh = Eh(i) associated with Th(i) .
In the sequel, we omit i as in the proof of Lemma 3.

We can find asequence {uh* } such that

for

where

mh

{~h·}· 1
J J=

1 ~ j

<

m ,

has been already defined in this section, and ~. forJ

I ~ j < m+l are real coefficients.

Notice that u h* is a linear combi-

mhnation of the first (m+l) functions of {~hj}j=l' and that we can always find u h* by an appropriate choice of a j .

As in the proof of thepreceding lemma, we can choose a suitable subsequence of {uh } toshow the existence of a non-zero function u*E:

H~(n) such thatfor all
for

1 ~ j

~

m

Here it is especially to be noted that the eigenvalues associatedwith Eh necessarily converge to AO as h ~ 0 due

to Lemma 3.

Theabove relations clearly contradict the fact that E is the eigenspacefor AO' and the proof is completed. lUI
6. Existence and Error Estimation of the Approximate Eigenpairs
In the preceding section, we have established some fundamental
properties of the approximate eigenpairs.

This section is devoted

to the proof that there actually exists a nice approximation to eacheigenpair of (4), except that for the zero eigenvalue.

12

Essentially,we will rely upon the standard technique, that is, the use of thespectral projection.

In our description, we will be mainly concernedwith the idea or the outline of analysis.
Let us consider, in the complex plane, a circumference r with
its center and direction taken as AO and anticlockwise, respectively.
The radius of r is fixed so small that there is no element of 0,
except AO' either on or inside r.
for each A e rand fE

Consider the following two problems :H~ (n), find u e H~ (n) and u h e x~ such thatv u- e HO(n)
2
,

(17)
and
(18)

respectively.

For each A e r and f E: H~(n), we can show that uhexists uniquely in HO(D), and that u h also does uniquely in Xo when2

h is small enough.

These results are derived by the use of theory ofcompact operators, or the approximate compactness properties established· in Lemma 2.

Although we have not stated explicitly, we haveconsidered real functions only.

Since A is a complex number, u andu h are in general -complex-valued functions even for real f, and hence
we are obliged to concern ourselves in complex-valued world.

How-

ever, the final results can be again expressed in real-valued world,
as will become clear later.
Now we can define the following two operators for each A E r:
Q(A) :H~(n)

(19)

-+-

H~(n)

Q(A)f = ue H~(n) of (17) for each13fe. H~(n). '(20)Qh(A)f = u h E: X~ of (18) for each

It is to be noted that Q(A)f and Qh(A)f (for h small enough) are
uniformly continuous in A E r for each f, the continuities being
those with respect to the metrics of H~(n) and x~, respectively.
Notice also Lemma 5 to be presented later.

( 21 )

Pf

Define= 2; i J Q (A) f drr

(22)

where f is an arbitrary function in H~(~), and i the imaginary unit.
We can show that P is actually a mapping of H~(n) into E such that

(23)

P u

=u

for anyu E

E •On the other hand, Ph is a mapping of H~(n) into Eh well defined forh small enough, but does not possess the property corresponding to(23) in the strict sense.
In the following, we will present three lemmas without proofs.
The first one may be proved with the aid of Lemma 2, while the second
and the last one are known results.

14

Lemma 5

Let A be an arbitrary number on'r.

For each f~

L2 (f2)H- I (fl) (= the dual of H~ (fl) ) , consider u h e. X~ such thatand g Ev u- (6 Xoh '

(24)

h

where for the expression (g,u h ) we have used ( , ) as the pairingon H-I(fl)x

H~(f2).

Then for h small enough, u h exists uniquely in

X~, and satisfies

where t'he positive constant e can be independent of A, h, f, and g.
Remark I

The corresponding properties hold for the continuousversion of (24).

Lemma 6For each A e rand f E H~(fl), Q(A)f(Mizutani [7])belongs to H4 (n) "~~(n) with the estimation

II Q(A)f "4 ~ ell f "2(26)

where e can be independent of A and f.
Remark 2,

2

Notice that f belongs to HO(O) and 0 is a non-obtuse

polygonal domain.

Lemma 7
Q(A)~

(27)

(Scholz [n])-For each A e rand f e H~ (n), consider

and Qh(A)f when h is small enough.

Then

II Qh(A)f - Q(A)f" + h l / 2 1ln hi" ~hQh(A)f -l1Q(A)f"~ e h lIn hi 211 f II 2 '
15

(28)

where C can be independent of A, f, and h.

Now we can follow the standard procedure.
4

Notice that ~i = p ~.1

2

E6 H (n) () HO(n), and consider(29)* and ~i satisfy the ·estimations corresponding
Then we can show that ~hi
to (27) and (28).

Therefore, for h small enough,linearly independent and belong to Eh.
= mh =< m from Lemma 4, and hence*~hl'.·.'

*

~hm

are

On the other hand, dim Ehm{~hi}i=lcan be regarded as a basisof Eh.
Notice that mh = m and consider the basis {~hi}~=l introduced
in section 5.

Then each ~hi can be uniquely expressed as

~hi'"

(30)

*
L
CL •• ~h .
j =1 1) )m(1 ~ i ~ m)by choosing the coefficients CL ij (1 ~ i,j ~ m) appropriately.

DefinellJ i bym(31)llJ·'"
1

L CL 1)• • ~.EE)

j=l(1 ~ iwhere CL •• are the same as those in (30).
1)

~ m)

Noting that CL •• are uniform~ly bounded with respect to i, j, and h, we have

16

1)

,

(33)

where C can be chosen independent of i and h.
The above estimations implicate that there exists a function
1jJ. e1

In general,

E sufficiently close to each <Ph"1

orthonormal in the sense
has been chosen such that(~1jJ.,~1jJ.)

1J= 0 ..•
1J(~hq,hi,6hq,hj){1jJ i }~=1

is not·However, the basis {q,hi};l

= 0ij' and hence we canfind an appropriate orthonormal basis of E, each basis function ofwhich is close to one of q,hi (1 < i < m) , cf. Chapter 6 of Strang==and Fix [10] . On the other. hand, the error of each Ahi can be evaluated by the use of the relation Ahi = b(q,hi,q,hi)/IIA~hill 2 .
The final results can be summarized as follows.

Theorem 1

For h small enough, we can choose a basis { q,i*}mi=l of E

such that**= 0 ..

(34)(~q,.,Aq,.)(35)

II Llh"'hi - Ll"'i II ~ Ch l / 2 1ln hi

J

1

II.

A..

1J

II.A..*,(36)

where 1 ~ i,j < m, {q,hi}~=l 'is the basis of Eh defined in section 5,
and C can be chosen independent of hand i.

Notice that the choiceof { q,i*}mi=l may differ with the triangulation Th.

As for thevalue approximation, we have

(37)

for 1

17

=< m .

eige~-

Remark 3

At present, it is not certain whether we may expect theestimations'" * II <II '"'I'hi - 'I'i

C h lIn h 12

for 1 =< i =< m .

7. Concluding Remarks
We have performed error analysis of a mixed finite element
model applied to linear buckling analysis of thin elastic plates.
The obtained order estimates of errors of the approximate eigenpairs
have been summarized as a theorem.

We have only sketched the proofs,

and the detailed analysis will be reported elsewhere.

The author

believes that the principles and techniques established in this note
will be available for the analysis of various finite element models.
This work is in part supported by the Grant-in Aid for the
Scientific Research from the Ministry Education.

The author would

like to express his deepest appreciation to Mr. T. Hanada andProfessor T. Ushijima of the University of Electro-Communications
for their valuable advices and discussions.

References[1] Kikuchi, F., A finite element method for plate bending analysis
by decomposition of differential operator, J. Nucl. Sci. Technol.,~(1971) 597-599.

[2] Ciar1et, P.G. and Glowinski, R., Dual iterative techniques forsolving a finite element approximation of the biharmonic equation, Comput. Methods Appl. Mech. Engng,~(1975) 277-295.

[3] Ishihara, K., The buckling of plates by the mixed finite element18

method, Memoirs of Numerical Mathematics,~(1978) 73-82.

[4] Miyoshi, T., A finite element method for the solution of fourthorder partial differential equations, Kumamoto J. Sci. (Math.),~(1973) 87-116.

[5] Grigorieff, R.D., Discrete Approximati~n von Eigenwertproblemen,
I. Qualitative Konvergenz, Numer. Math., 24 (1975) 355-374, II.
Konvergenzordnung, Ibid, 415-433, III. Asymptotische Entwicklungen, Numer. Math.,~(1975) 79-97.

[6] Scholz, R., A mixed m~thod for 4th order problems using linear
finite elements, R.A.I.R.O., Analyse numerique,~(1978) 85-90.

[7] Mizutani, A., On the regularities of solutions of biharmonicequations in domains with angular points -- notes on a paper ofKondrat'ev, KOkyu-roku of Reserach Institute of Mathematical
Sciences, University of Kyoto, No. 329 (1978) 2-9.
[8] Ciarlet, P.G., The Finite Element Method for Elliptic Problems,
North-Holland (1978)[9] Kondrat'ev, V.A., Boundary problems for elliptic equations indomains with conical or angular points, Trans. Moscow Math. Soc.,
16 (1967) 227-313.
[10} Strang, G. and Fix, G.J., An Analysis of the Finite Element
Method, Prentice-Hall (1973).
[11] Dym. C.L. and Shames, I.H., Solid Mechanics, A Variational
Approach, McGraw-Hill (1973).

l12J Canuto, C., Eigenvalue approximation by mixed methods, R.A.I.R.O.,
Analyse numerique, ~ (1978) 27-50.
[13] Kikuchi, F., A mixed finite element model for fourth ordereigenvalue problems, to appear.

19

I

The Initial-Value Adjusting MethodforSolving Problems of the Least Squares Typeof
Ordinary Differential Equations

By

Taketomo MITSUI*Abstracts
An application of the initial-value adjusting method forthe problem of the least squares type is considered.
gencepropert~of the method is discussed.

Conver-

An illustrative

numerical example is given.

§l.

The Initial-Value Adjusting Algorithm.

We are concerned with numerical procedures solving the
problems of the least squares type for ordinary differential
equations as following:
Find a solution of the differential equation~~ = X(x,t),(1.1)

a < t < bwhich minimizes the value

(1.2)
where

x

and

given points on

X

are real n-dimensional vectors,
I = [a,b],

tj

are

a=tl < t2 < ••• . < t N=b '*Research Institute for Mathematical Sciences, Kyoto University,
Kyoto, 606 Japan.

21

and

dj

are given nxn matrices and n-dimensional vectors

respectively.

As far as the equation (1.1) is stable,

problem is equivalent to find an initial value~hesolutionx(t)starting fromn

Let us define a linear operator

that

th~

J

(1.4)

such thatJ.

L : C(I) ~ mnN

suchx~C(I)

that for

Denoteminimizes

n

the

is

nN-dimensional vector

by

d.

Note

represented by

1 t
J = 2 {Lx - d}{Lx - ~}.

After the initial-value adjusting method proposed by
Ojika and Kasue [7J,

our algorithm can be expressed in thefollowing steps.
Step O.

Choose a suitable perturbation parameteran initial value
Step 1.

nand setnOEm ,

E

and

k=O.

Compute the numerical solutionxk(t)of (1.1)

xk(a)=n k , and obtain the resultIt
J k =2 {Lk-~}{Lk-~}·for the initial conditioning value

Lk=Lxkandvaries slightly in comparison

Step 2.

If the value'

Jk

J k_l

(i. e. Jk_l-J k

is sufficiently close to 0 by thewithcriterion given in advance),

terminate the

wise,

(If

go to the next step.

k=O,

iteration~

Other-

skip this step.)

Step 3.

Set

Step 4.

Compute the numerical solution

j=l.

(1.1) for the initial conditionmeans the j-th uint vector of

22

Yk(J)(a) = n k +~e
"- j.

mn.

of
Here

Step 5.
until

Replace

by

j

j+l,

and return to Step 4

j=n.

Step 6. Determine the nNxn matrix
matpi~)

such that(1.5)

S(e:;x k )S(€;x k ) (the adjusting= (-kn,y(l)_n,} ... ~n,y{j)-n, } •••
€
k
k'
'e:
k
k'
,
~ n,y ~ n) -n,k} ) .

Step 7. Determine the initial value

nk + lfor the next

iteration by(1.6)
Then replacekby, k+l,

and return to Step 1.

ComputationaZ Remapks. In the above process,ical

the numer-

integration of the differential equation (1.1) are

carried out by a 'suitable step-by-step method,
the Runge-Kutta method.
is an nxn symmetric,

Since the matrix

for example,

tS(€;xk)S(e:;x k )positive definite matrix,

the square-root-free Cholesky's method is preferable for the solutionof the linear, equation, in. (1.6).

We shall denote the Euclidean norm of an n-dimensional
vector

x

by

axJ.

C(I)

stands for the Banach space of

vector-valued continuous functions on

I,equipped with

the norm= maxI x(t H .
tEI
Cl(I)

means a subset of

C(I)
23

of continuously different i-able functions on

I.

mm ~ mm'

The norms for matrices

andother linear operators should be taken as the induced norms
bythe .corresponding vector norms.
Let2

be the domain of the tx-space bounded on x,

intercepted by two hyperplanes
points of ~on the hyperplanes

to be included in ~plane.

t=a

and

t=a

t=b.

The boundaryt=b

are supposed

andand to make an open set on each hyper-

Put

D = {XECl(l); (t,x(t)for tel},D'- { x4l!C(l); (t,x(t)for tel}.

We shall analyse the initial-value adjusting algorithm
within the above framework.

Because of space limitations,

we shall omit all the proofs of the statements,
be shown inwhich willother papers [5J,[6J.

§2.

Some Preliminaries.

We shall adopt the following assumptions to the problem( 1 • 1) and (1. 2) .
Assumption 1.

X(x,t)

is defined and suffiaiently aon-tinuousZy differentiabZe ~ith respeat toits

de~ivatives

a~eaontinuous~ith

Let us consider an operator ~tion ofCl(l)on~ . . X

respect tomapping

mn

t

on

and
I.

to a func-

along the flow generated by the differential

equation (1.1):~t Q'in = X[~n(t),tJ,
~n(a)

x

a<

= n.

Define a function

J(n)

by

24

t<

b,

I t{LOfn -ti}{:n,~n-iIl}.
ItI"'C.

J(n) = '2

(2.1)
Let

~(t;x)

stand for the matrizant of the linear

homogeneous matrix differential equationa<:L~(·;~n)

We shall call the matrix
it by

t<:

lJ.

the G-matri~ and denote

G (n) •

Assumption 2.

(a1n*)E-J;)1

e~istsTheren*a vectorsuch that"n*ED,(2.2)

anddet tG(n*)G(n*)~o

(2.3)

hoLd.

Furthermore,there

0 <: 1£,1 ~ £'such that for

41i..
dt -- X ( y,t),

a <:

t <:

e~istsa positive constant-£the initial value problemb,yea) = n*+f:.e.
J

has a unique soLution in

point of

J(n).

L

D.

of

on

Let

LO

For the matrizant

for everybe the bound of the operator-norm~(t;x)(2.4)

is evident.

the estimation
for

Since

j.

the e~act (locaL) isolated minimaln*

We shall callD

D is open in
25

XED'

CI(I),

for a positiveconstantin thewe can take a positive numberMIll-neighbourhood

Bll

of

II

such thatn*Un-n*U <= ll},.

~(n)

the inverse of compound G-matrix defined byt= G(n)G(n~,exists and have the estimation

ReplaceIIby a suitable value less than II .if necessary,then it is possible that" ~ n- ~n*1/ C~ <5 0n~Bllholds fornumberand small positive

We shall fix the

00'

ll.

Lemma 1.

Forn€Bll~

putx(t)=~n(t).

y(j)(t)in

Step 4 of the algorithm has the following expression.
y(j)(t) = X(t)+E{p(j)(t)+v(j)(t)}

(2.6)

where

Cj'(j) (t)onis the j-th column vector ofv(j) (t)satisfies the differential equation(2.'1)

dv(j) (t)
dt

1

(.)

('J

= ~ X[x(tJ+cl<J' J (t)+v JI~~(t;x)

and(tJ}, t]-X[x(tJ ~ tJ}a <t< bsubject to the initial condition( 2. 8)v(j)(a) = O.

Moreover,for arbitrary small

constantC*such thatv(j)(tJ

26

lEI

there exists a positive

satisfies(2.9)

Throughthe equation

S(E;~n) -G(n) = L(v(l) , ••• ,v(j) , ••• ,v(n»,Lemma 1 immediately implies that there exists a small positive
numberEO

0 ~. I E I ~ EO

such that for

the estimation(2.10)
holds in

BLl.

Furthermore,

for

0 <

I E I ~. EO

andn4iB l\tS(E;i1n)S(E;~n)

the inverse of the compound adju8ting matrix
exists and(2.11)
l)olds.
Corresponding to our iterative process,mappingmn

an operator~

into itself is defined by the following:

The domain ofA

is identical to

for

BA.

neBLl"The iterative process is simply represented by(2.13)

k=0,1,2,···.

Thus our aim is concentr'ated on the analysis of the operator

4.
§3.

Fixed Point of the Operator~.

It is noteworthy that the exact minimal point

27

n*

is

~,

not a fixed point of

because a fixed point of~must satisfy the equationi.e.
(3.1) ,

from (2.2)

On the other hand,n*satisfiestG(n*){L~n*-~}=O.
S(E;~n)is surely an approximation ofidentical to that as far as
Therefore,A

point ofexists.

Define a function. (n, E)such that

Lemma 2.

Define( J. 3)but is not

E~O.

we should answer the question whether a fixed
We shall show tl1at it actually exists

in some neighbourhood of

for

G(n),

~(n~ 0),.

n*by the implicit function theorem.

~(n,E)

byand~(n~ 0)by

lim J(n~ £) ~
£+0

then
(3.4)

hoLds.

The following two lemmas say how the adjusting matrix

28

e:

depends on

n

or

when

they make a small variation.

(3.5)

holds.

n, e: 1 and e: 2 .

Here the constant does not depend on

Lemma

4.

3

For sufficiently small numbers

and

E,

the

equation(3.6)

holds i fn, n+~EB3

le:I ~E.

andlinear mapping ~hich maps

~E8n

P(t;~n,e:)Hereto an nxn matri~~ith components of continuous functions

is aP(t;~n,e:)~I.

on

Lemmas 2 and 3 and the implicit function theorem imply
the desired result.

Theorem 1.

e~ist

There

e:

positive numbers

such that,for anyJrn, e:) = 0

has a unique solution

lin - nltu ~ 1::. 1 'andle:I ~e:l'satisfying,..
n -+ nitas,..

n = n(e:)

E

1::.1

.to-

.

and

e: 1

the equation

l.n the ballo.

We shall denote the fixed point ofAby~ = ~ (E)whose existence has been guaranteed by Theorem 1 and call it
the appro~imate minimal point ofapproximation of

n*

J(n)

and tends to it as

because it is ane: -+ O.

§4.
Since

holds if

Convergence of the Iterative Process.
is defined by

Jln-n*1I <== 1::..

An - n = n - '"n~n(n,e:)

Here
to

n.

~('n, e:) =0,

the equation

Then we have{ t S ( e: ; ~ n ) S ( e: ; «f n) }-1 J( n , e: )~n (n,e:)(n-n)}+o(IIn-n. fl ).

stands for the Jacobian of~with respect

Utilizing the results of the previous section,

we

can investigate the estimation for the termtS(e:; ~n)S(e:;Y;n) (n-rD- ~(n,e:) (n-n).
Then,with the estimation (2.11),

the contraction mapping

principle brings the following

Theorem 2.

There exists a positive number

1::.2

suah thatthe iterative proaes8 (2.13) starting from anynin thel::.2-neighbourhood ofRemark.

,...

n.

converges to

The results in the previous section also gives

the explicit formula for the Frechet derivative ofBy

the consideration o.f

At (~) ,~at~.

we can conclude that theconvergence of the iteration would not be expected to bequadratic in the neighbourhood ofchosen to be apart from zero.

30~

while a fixede:isNext,""nfor

we shall investigate the approximation order of
The following improves the statement of Lemma 3.

n*.

(4. 2)

hoZds,~hepep(t;~n,E) is an nxn matpi~ ~ith componentsof continuous functions onI.

By virtue of Lemma 5 we obtain informations about thederivatives ofJ(n,E)

bourhood of (n*,O),

Theorem 3.
{(n,E);

with respect tOo:E

in some neigh-which implies the following

In the neighbourhood

Un-n*\I ~ 6. 3

andlEI ~'E3},nee) ~hose e~i8tence is guaranteed in Theorem 1 has theestimation such as
(4.3)

§5.

An Illustrative Example.

Let us consider the following problem originally mentioned
in [4],
reacto~which occurs relating to the tubular flow chemicalwith axial mixing.

The differential equation is
'J

(5.1)

d'-x
d
2
6~ - l2x = 0
dT2 dT'

By the transformation

31

o < T < 1.

(~.2)= 1 - 21',t

the equation (5.1) can be reduced to the following differential equation:2dX
d x
2 + 3dt

2 · 0,- 3x

dt-1 < ,t < 1.

The constraining condition of least squares type is asfollows:
10

(5.4)

r

{x(t j )j=O-

y j ,}2. minimum,

where.
(5.5)(j • 0,1,··· ,10)0.2j - 1.0

tj •(j = 0,1,·· ·,10)

and

are given in Table 1.

Table 1.
j

0

1

2

3

4

5

Yj

0.387270.39476

0.41305

0.43862

O~470l7

0.50764

6

7

8

9

10

0.55172

0.60372

0.665590.740120.83129

This problem is equivalent to the example in [2].
We shall rewrite the problem into the vector form
as (1.1) and (1.2).
dX 1

dt
(5.6)

The equation (5.3) is transformed into

= x2 ,-1 < t < 1.

dX 2
2
- = 3xi - 3x 2 ,dt32x = [:~1,then (5.6)

l'-1 < t < 1

and the functional valueJ

to be minimized isHereare the same as in (5.5),Letxbe the vector

is equivalentto~~ = [ 3X2 -2 3x( 5 · 6) ,Xltj (j=O,l,···,lO)

matrices

Lj

2(J=O,l,· _. ,10)

LO=L l =-· ·=L IO =and the vectors

the

are

[~ ~) ~d j (j=O,l,···,lO)

are given by

d = [Yj0.0 1 .
j

The results of the numerical computation carried out by
FACOM M-190 in the Data Processing Center, Kyoto University
are shown in Tables 2~4.

For the numerical integration ofordinary differential equations the Runge-Kutta-Gill method
programed by T.Ojika was used. All the calculations werecarried out in the double precision arithmetic.

. 33

Table 2. Results obtained by the iteration.
itcration
times

1

0

x2

xl

2

3

4

xl

x2

xl

x2

xl

x2

xl

x2

t o"'-1. 0 0.5

0.0

0.37806769956

0.10806510612

0.38702086627

0.00261251861

0.38727177735-0.00004288586

0.38727191330-0.00004431630

t 1=-0.8 0.51252090985

0.11507113359

0.40171204036

0.12821189883

0.39490532424

0.07015010707

0.39476040385

0.06869739433

0.39476032659

0.06869661205

t 2=-0.6 0.54334736769

0.18909213769

0.42937707807

0.14869156970

0.41341909766

0.11229585764

0.41304637216

0.11139645273

r,.4130461722'1 '0.11139596872

t 3"'-0.4 0.58721750119

0.24863472740

0.46135283076

0.17162147360

0.43913778607

0.14393243601

0.43861177501

0.14326111462

0.43861149266

0.14326075367

t 4"'-0.2

0.64281390289

0.30821629103

0.49830337268

0.19871943583

0.47081948324

0.17287829974

0.47016701758

0.17226427328

0.47016666726

0.17226394342

t 5=0.0

0.71115159268

0.3774847.7983

0.54123626991

0.23176626442

0.50841638291

0.20369262002

0.50763820223

0.20303615170

0.50763778440

0.20303579926

t 6-0.2

0.79503126629

0.46516837677

0.59155143668

0.27296261250

0.55263785706

0.23962085391

0.55171808024

0.23885031164

0.55171758641

0.23884989812

t 7=0.4

0.89914650517

0.58192807415

0.65116185255

0.32531012072

0.60480645049

0.28369352839

0.60371555428

0.28274087195

O.6017l49hR6J

0.28274036083

t S"'0.6

1.0307690460

0.74350525298

0.72270288155

0.39313865121

0.66689448238

0.33949736231

0.66558820885

0.33828049978

0.66558750766

0.33S27984706

t 9=0.8

1. 2012163950

0.97571861039

0.80987497209

0.48294851916

0.74170948R34

0.41193878825

0.74012428084

0.41034340656

0.74012343005

0.41034255099

tlO=1.0 1.4287209566

1.3241632086

0.91800946931

0.60486691312

0.83325492973

0.50828488609

0.83129911340

0.50613875295

0.83129806389

0.50613760232

Jki

1. 16 X l0- 2

4.93 xlO- 1

5.81 x10- 6

l.12 x lO- IO

1.10 X l0- 10

i
I
I
I

I
convergence criterion-7

1.0x10step size for the Runge-Kutta method
E •l.Oxl0-S

0.003125

Table 3. The results for different starting values (I)(E"';1.OXIO- 8 )
(a)

(b)

(c)

Starting { xl (-1.0)

0.5

0.4

0.1

Values

0.0

0.0

0.0

4

3

5

X2 (-1.0)

Iteration timesJO

4.93 x lO- l

4.60XlO- 3

1.16

Jl

1.16xlO- 2

1.54 xl 0--6,

3.2l x lO- l

J2

5.8lxlO- 6

1.11xlO- lO 3.45 x lO- 3-

J3

1.12XlO- lO 1.10xlO- lO 5.70 x IO- 7

J4

1.10xlO- lCf

J5'~::::==-

1.10xlO- to
1.10xlO- lO

For the above three cases all the converged
values at each tj coincide to elevenfigures.
The values are shown in Table 2.

convergence criterion

1.OxlO- 7

step-size for the Runge-Kutta method

35

0.003125

Table 4.
The results for different starting values (II) (E=1.0 X10
(b)
(a)
(c)

Starting
values{Xl(-1. 0)

0.5

0.4

0.1

x 2 (-1.0)

0.0

0.0

0.0

~6

4

3
4.60x10- 3

1.15

J1

4.93 x10- 1
1.16x10- 2

1. 54 x10- 6

3.21 x10- 1

J2

5.81 xl0- 6

1.11X10- 1O

3.45 x10- 3

J3

1.12 xIO- 1O

1.10x10- 1O

5.70 x10- 7

Iteration times
Jo

J4
J5

5

--------------1.10 x10-10-'---------'--)

1.10x10- 1O
1.10x10- 1Oy)

Por the above three cases all the converged
values at each tj coincide tb elevenfigures.
The values are following:~'---..
j,
-1. 0

-0.8
-0.6
-0.4
-0.2
0.0
0.2
0.4
0.6
0.8
1.0

xl

x2

0.38727191327
0.39476032658
0.41304617226
0.43861149265
0.47016666726
0.50763778440
0.55171758641
0.60371496863
0.66558750766
0.74012343005
0.83129806389-0.00004431619
0.06869661211
0.11139596875
0.14326075369
0.17226394343
0.20303579926
0.23884989812
0.28274036083
0.33827984706
0.41034255099
0.50613760232

convergence criterion

1.0X10- 7

step-size for the Runge-Kutta method

36

0.003125References[lJ

Banks, H. T.

& Groome, Jr., G. M., Convergence theorems

for parameter estimation by quasi1inearization, J. Math.
Anal. App1., ~(1973), 91-109.
[2]

Fujii, M., Numerical solutions to problems of the leastsquares type for ordinary differential equations by theuse of Chebyshev series, Bull. Fukuoka Univ. Education,
27, Prat III(1978), 15-26.

[3]

HpacHoc8nbcKHH, M.A., BaHHHHKo, r.M. H ~p.,
P8W8HH8[4]

on8paTopH~X

ypasH8HHH,

npH~nm8HHo8"HaYKa", MocKBa, 1969.

Lee, E. S., Quasi1inearization and invariant imbedding,.Academic P., New York, 1968.

[5]

Mitsui, T., On the convergence of the initial-value
adjusting method for nonlinear boundary value problems,
submitted to Publ. RIMS, Kyoto Univ.

[6J

Mitsui, T., The initial-value adjusting method forproblems of the least squares type of ordinary differential equations, in preparat~on.

[7J

Ojika, T. & Kasue, Y., Initial-value adjusting method
for the solution of nonlinear multipoint boundary-value
problems, to appear in J. Math. Anal. App1.

[8]

Urabe, M., The Newton method and its application toboundary value problems with nonlinear boundary conditions,U. S. -Japan Sem. Different. Funct. Eqs., Benjamin,
New York, 1967, pp383-4l0.

E9]

Urabe, M., On the Newton method to solve problems of theleast squares type for ordinary differential equations,Mem. Fac. Sci., Kyushu Univ., Sere A,

37~(1975),173-183.

A Direct Method for Solving Two-dimensional
One Phase

Ste~anProblems

By
Vitoriano RUAS B. SANTOS*

Sutnrnary
A new algorithm is proposed for solving the Stefan problem
for the case when the initial region occupied by the medium ina given phase is a starshaped

two-dimensional domain.

The evolution or this region as time increases is determinedby plotting the rree boundary directly.

At each time step this1s approximated by a polygon in the interior of which the heat
equation is solved with piecewise linear finite elements.
NUmerical experiments indicate that the method is stable
and that convergence is to be expected, under suitable assumptions
on the regularity or"both the boundary and initial data of the
problem.

1.

Introduction

We shall be concerned here with an algorithm for solving
one-phase Steran problems in two-dimension space.

The basic

tool ror the algorithm is a method or automatic generation oftriangular rinite element meshes that the author initially proposed
for solving boundary value problems derined on stationary star-* Computer Science Department, Pontif!cia Universidade Cat61ica
do Rio de Janeiro, Brazil
andDepartment of Mathematics, University or Tokyo.

39

shaped domains [8].

In that work many technical details of its

implementation are given, while regularity properties of the sogenerated mesh have been analysed in [9J.
This automatic triangulation process is used here both forplotting the free boundary and for adjusting the mesh, so thatthe heat equation can be conveniently solved just in the domain
occupied by one of the phases of a certain medium undergoing achange of phase.

This is what we call a direct method of solutionof a Stefan problem.
Actually, our algorithm generalizes the one proposed byMori [6] for the one-dimensional case.

In his work the freeboundary is determined after each increase in time of a fixed
value

6t,

while the mesh containing a fixed number of intervals

moves according to the new shape of the domain.

As it was well

remarked by Prof. Fujita [3J, our triangulation process could be
the appropriate means of doing the same thing in the two-dimensional case, provided that the outer boundary of the initial
domain is starshaped and that it remains starshaped as it evolves
in time.

Of course, since the domain is in itself one of theunknowns of the problem, the latter condition cannot be satisfied
a priori.

However, one might expect that, in most practicalcases, the domain becomes more and more regular as time increases,in the sense that it actually remains starshaped, althoughrigorous proofs are not,available yet.
Let us also say that Mori applied with success a slight
generalization of his one phase algorithm to a two-dimensional
strip domain [7].

Later Bonnerot and Jamet proposed the direct

solution of one phase two-dimensional Stefan problems by time-10

1

space finite elements [5] .

Though their method works well inmany cases, as shown by their examples, they did not propose any
solution to theincreases.

proble~of dealing with the spatial mesh as time

We believe this may cause either numerical incon-

venience or difficulties of implementation for practical cases,
unless a suitable solution of this problem is a part-of thealgorithm itself (see comments in Section 5).
Finally we note that a number of reasonable algorithms for
indjrect solution or the two-dimensional two phase problem areavailable at present.

As a signiricant example or those wemention the work or Ciavaldini [1J, that is based on a suitabletransrormation or the heat equation into a semilinear problem.
This is solved on the rixed bounded domain occupied by themedium in both phases, and the position of the free boundary isthen determined with the help or a runction or the solution of
the transformed equation.

By introducing appropriate modifica-

tions, this approach can also be used ror solving the pure one
phase problem.

However, in this case, although one can linearize

the discrete problem, it is necessary to perform the calculations
ror a domain much larger than the one occupied by the only phase
undergoing changes or temperature.

This can be impractical in

the case where the domain evolves to very large regions comparedto the initial one.

2.

The Continuous., Problem

As we will comment in Section 5, with the eventual introduction of

simpl~I Actually,

modifications, our algorithm can be applied totheir algorithm gives approximate problems very

close to Mori's in the cases he.cbrtsiders.
41the solution of a very wide class of Stefan problems.

However,

we shall confine ourselves here to the solution of the special
case of continuous problem defined on domains whose growth is
not at all limited a priori, as described below.

Our algorithm

seems to attain the best of its efficiency compared to others
particularly in the case where the domain may become much larger
than the initial one.
In the sequel we introduce the formulation of the one phase
Stefan proble~ we shall consider here, as a special case of thegeneral multidimensional one treated by Friedman [2J.
Let a certain medium exist in two different phases

1

and

2.

Phase 1 occupies initially a certain bounded region

the

x-y

plane, whose boundary consists of two disjoint curves

rO

and

r*.

t

rO

= ° and r*

applied.

is the interface separating both phases at time
is a fixed curve where heat sources are to be

We assume thatr*lies in the interior of

we denote its own interior byn*.
nO u n*rO

and

n* are starshaped
domains such that there exists one point 0 E n* for which thesegments joining it to any other point either of n* or nOun*
We assume also that bothandlie completely in' those sets, respectively.

n*

and

nO

An example of setssatisfying this assumption is illustrated in

Figure 1.

Figure 1
4L

As time increases a change of phase 2 into phase 1 will
occur in the medium lying in the region beyondrD.

So, thecorresponding problem will be defined on the region

R2 - n* .

As a practical example of this situation we have the problem
of melting ice by heating some already melted region

a pipeline represented by
We shall denote by

nO

with

r* (Figure 1).
net)

and

ret)

the domain occupiedby the medium in phase 1 and its outer boundary respectively, at
time

t,

that

ret)

~ O.

Clearly we haveu r*is the boundary of

t

n(O)

= nO,

reO)

=rO

and

net).

Let

(1)

cj)(x, y, t) = 0

be the equation of the free boundaryis such thatcj)(x, y, t) > 0

cJ>(x, y, t) <u n*.

the heating process takes place up to timeture

We assume that

Now, supposing that'T,for the tempera--uof the medium in phase 1 we have the heat equation:(2)

au = 6u
atin

Vt E [0, TJ

net)

with the initial condition:(3)u(x, y, 0) • uo(x, y)and the boundary conditions:

(4)

.....

.....

u(x, t) • g(x, t)or(v

denotes the unit outernormal vector offorx E r*,

cj)° for (x, y) E n(t).U n* and(x, y) ~ n(t)

for

ret).

Vt E [0, T],
43

r*),

and
-+

Also for each timeon-+

x ~ r (t ) ,foru(x, t) = 0~t~ [0. TJ.

t, the so-called Stefan condition holds

ret)

(6)

( V <P ,Vu)

1-+

XE r(t)

where

a<I>
= Kat·

is a positive constant andKproduct in

(.,

.)

denotes the scalarl1 2 •

Thus the problem we want to solve is finding

@ and

u

to sat i s fy (2) -- (6).
As we have already mentioned, the algorithm we shall employ
for solving this problem, is constructed upon a process of
automatic triangulation.

This in turn is based on thA represen-tation of the free boundary by an equation in polar coordinates(p, e).

So it will be useful to take$to be an exprAssionof the form:o-s(8, t).

<:> ...

In order to do so we must choose first of all a suitable origin of
coordinates ~ying inn*.

According to the assumptions that wehave made it is possible to choose such an origin so that both
r*

and

ret)

can be represented in polar coordinates 2 •it makes sense to write

~

I

Thus

in form (7), and we have:dS

V<I> = (cos e + S 19 sin e,
and

2 Strictly speaking
shaped.

For

th~nO U 0*should be "non-singular" star-corresponding definition and othersee [)J.
44

detai~5andas-= - ft"Thus (6) becomes:

(8)
where

u

uJ

I asa
[a
=ap-sae~ l(p,8)£r(t)andare, respectively, derivatives of

two perpendicular directions (see Figure 2).

uin

The meaning of thefirst derivative is clear whereas the second one is simply
I au

p ae'
On the other hand, according to (5), the tangential derivaau
tive aT on r(t) must vanish.
So denoting byva

the angle between the polar radius and by

the OUtL t' normal on

we have:r( t) ,

0....:;.---'----------FIgure ;2(9)~u sina +au cosa -= O.

apaa

Using (8) and (9) and taking into account thattan a -=

as
s1 as
'

we get:45as)~Jaul
asl]r -I- (1S ae
ap p =s ( e , t) = - Kat·

( "]C) )(10) is the expression of the Stefan condition whenuisgiven in polar coordinateR.

3.

The Algorithm

We shall now describe the algorithm that we propose forsolving (2) ~ (6) numerically.
First of all thp discretization of (2) -(5) is performed

by standard methods applied to the corresponding variational
formulation, namely:() u v dx = - /
IT

!

(11)

net)

\lu\lv dx

net)

if the Dirichlet-type boundary condition (4) holds, or(ll)~ /~~VdX = -/net)

\lu\lvdx+!

net)

r*

g'vds

VVEV', uEV'

in case the Neumann-type boundary condition (4)' holds.
andV~

are the subspaces ofeverywhere on

u r*

ret)

is the linear variety of
w =g

that

on

Hl[~(t)Jand on

V

ret),

Here

V

of functions that vanish
respectively, andconsisting of· functions

w

V

such

r*.

For the discretization in space of (11) and (11)' we use
a triangulation of
spacesginterpolate
lying onV'hor

Functions

r* ,

n(t)

and
gh

of piecewise linear continuous functions.
g'andwill be replaced by their piecewise linear
gh

which coincide with them at the nodes

respectively.

approximate solutionofon which we construct finite elementuhThus in case (4) holds, thewill belong to-, i.e., the linear variety of

V.

Vh ,

discrete analogue

such thatimplies

r*.

= gh

wh

on the inner polygonal boundary approximating

If (4)' holds, integration will be performed along thispolygon instead ofr*.

For the discretization in time we employ standard schemes.
This means that once obtained a linear system of ordinary differential equations after discretization in space:

(12)

o=u

we start with

uh

U

h

o

being the Vh-interpolate ofhu h1 ' u h2 ' ... , u hn 'and calculate

Ln = wt

where

o '

n

•.. , approximations of+ (l-w)t n- l' . andw E [O,lJuc'u h ( t)is the scheme

parameter.
As a matter of fact, we take a constant increment of time
fit,

so thattn = nflt, and approximaten (t n ), say byIntegrations are performed on a weighted domain given by

Wg~+ (l-w)n~-l.
t,M~Matricesandfor the approximate domain

argument

LnA~

do not really depend onchanges discretely." Thus· theabove only accounts for this interpolation.

Before describing how we determine

we should give ashort account of the triangulation method that we use.
Letp

= sO (6 )coordinates ofNow, given
and

p

and

rO
e: > 0pand= s * (e )r*,be the equations in polar

respectively.

(possibly small), we choose integers

such that:

47

m

I
max
ee [ 27T ]

°,

s*(6)

-

p

SO(6)-s*(6)1< E
mand we define
M = 8(p +m).

nO

Let the boundary ofandr~whose vertices are the intersections of

e = 6j ,••• ,8p-l,= 0,1, ... ,M-l,

rO

and

J t3

and

e = e~,
0* = i8*,ibe the domain bounded by

rO

and

with the lines

8*

be approximated by polygons

jrespectively, where

ej

=

and

i

r~,

and

= 0,1,
and8

being given by

8 = 7T/4(p+m)

(lIn
and

8* = 7T/4p .

(14)'

p =

Let

nO

S~(6)

be the equation ofhPkt

Now the verticeshrO

in polar coordinates.
hof the triangulation ofaredefined as follows:

(15)

k

:0:

p, p+ 1, ..• , p+m;

R,:o:1, 2, ... , 8k .

An illustration of the so obtained triangulation isin Figure 6 (see Sectio·n 4) for

p

=4

and

m = 2. 3~iven

For otherdetails see [9J.
The meaning ofEis establishing a relation between

3 This triangulation is slightly different of the one we
consider in [9J. The essential difference is due to the
exclusion of the triangles lying in n*.
48

p

and

m

so that the length of all the edges of the mesh remainof the same order as we refine the mesh.

In the ideal caseo(m- 2 ).

In this way we canshould be chosen to be at least

immediately conclude that the spacial step sizeasO(m - l )

cases

£

O(p-l).

or equivalently as
is justO(h)· O(m- l )h£is defi.ned

But in most practicaland still the above equiva-

lence could be verified.

r~Notice that by this process,is a polygon ofvertices and that the domain is divided inton/4,

of amplitude

basic sectors

where the same kind of partitinn 1s performed.

Also, in each sector we have
p+m

8

M

p+m

r~

triangles adjacent toboundary triangles, i.e.,along an edge; between neighbor-ing pairs of those boundary triangles we havetriangles, intersecting

r~

p+m-l

interposedby one vertex.

Let us now introduce the process for determining theposition of the free boundary

1, 2, •.. ,and consequently

First we denote bydirections

8j ,

j

0

Sj

n

rn
h

f2 h'

at the n-th time step,
the interior of

the polar radii ofr~=r~ur~.

r hO in the= 0,1, ..., M-l.

Now we assume that during the fixeci ttme j.ncrement

r11-1-moves to
h

n

6t,

in the way shown in Figure 3, which repre-

sents a part of the domain.

Figure 3
4yIIbeing the displacement of theboundaryr~nthe polar radii

0,1, ••• ,M-l,bound~ryin the direction

defining the newSj

are given bynn-l

Sj = S,j

The increments

n

AS j

Stefan condition (10)n+6s j ,

n = 1,2, . . . .

are calculated by discretizing theas follows:

a n-l
Uh

(16)

ap

J

n-l n-l 2(P )[1 + (Sj+l-Sj_i)j2Bs~-1 )n-l

is given by (14) anda,u n-lh(P. )are the vertices ofas shown in Figure 3.
dO
Jcalcqlated in two different ways according to the position of= sM_l'

where

1.

If

j

= k(p+m)~k

B= 0,1, .•. ,7

calculated along the edge

then it is simply

QjP j

= 6j

lying on the lipe 6(see Fig. 4a),Fig. 4a~lig.

4bJi'tgure 4
2.

Ifj;J k(p+m),k = 0,1, •.. ,7,

calculated along the segment

QjP j

then

is

(see Figure 4b)

lying in the interposed triangle between the boundary
triangles having

Pj

as a common vertex.

50

We note tha't due to the definition of our triangulation
process, the latter operation always makes sense.
Finally, keeping the number of triangles constant, we simply
define the new mesh by the automatic triangulation method (15)by replacing

o

sh

n

by

sh'

n > 0,

forn

P = sheS)

being ther~.

equation of the polygon

We should note that in our method each basis function
depends on time,
time step.

for the corresponding node

Denoting by

the coefficient of
node

Pi'

N(h)

uh(t)

Pi

moves at eachthe number of nodes and by~i'

associated with

~j

ui(t)

that is, with

we haveh• Nfh) (dU

aU

at

i=1

~) •

j (t)

9'j + at u jdt

Now recalling (12) we have:
Ah(t)

lit

{o(t)
r [ , if1
V"1 V~

a C'fj+'1 at] dX}

andThe mass matrix

Mh(t)is obviously symmetric butAh(t)is non-symmetric due to the second term in the integrand. Actually the non-symmetric component ofAh(t)is matrix

a fj

Cf'i at dx}
that is called by Mori the velocity matrix, since it accounts
for the effect of displacement of the nodes with respect to time.
Another point that is worth a remark is the following.
Since for this algorithm we do not increase the number of triangles, and since tpe edges of51

r~. have fixed length at every

r~time step, triangles closer toincreases.

will tend to elongate asn

Although from the accuracy viewpoint this may not be

a desirable situation, as far as the rate of convergence isconcerned, this violation of the classical angle condition for
unboundedndoes not yield any disadvantage.

Indeed, as it

has been proved by Jamet [4], the essential condition for maintaining the optimal rate of convergence of triangular finite
element approximate solutions is not the lower boundedness of
the angles, but the fact that no angle approaches
mesh is refined.

180 0

as theSuch a situation is not occurring here.

4.

Numerical Results

We have tested the numerical viability of the algorithm
proposed here which generally presented a good performance.

Somesingificant exa:mples of the obtained results are given below.
All the calculations were done in double preciSion on the
HITAC/8800-8700

of the University of Tokyo Computer Center.

For the solution of the systems of linear equations we used theGaussian method for band matrices.
Example 1:

For the given initial domain

nn = {(x, y)/l < p <. 2}
we take as an exact Rolution the function1
2
u = -(t + 2 - -p- )2

t+2that satisfies a non-homogeneous equation of the form:
au _ flu = fat.

The equation of the free boundary isp= t + 2. We preferred tochoose (4)' as a boundary condition, so that we can compare52computed values ofuat the stationary nodal points of ther*.

grid, i.e., those lying on
We calculate up to time

T = 1

so that the final domain is

Q(T) = {(x, y)/l<p<3L
Due to the symmetry of the problem we have performed thecalculations only for the sector
h = 11m

We define
We also takew - 1

0 <8 < n/4.

and we takeIn this way we choose

Table 1: Computed va;t.ues of

0

.u

for

~t

= h/5.

p = 1.

ht = 0.1

t = 0.5

t = 1.0

1/2

0.81940.8111

1. 0528

1.31751.3068n/80e: = 0).

as the scheme parameter so that we obtaina usual implicit scheme.

8

p = m (so, in (13)1.04240.81540.8114

1. 04421.32061.3161

1/8

0.81320.81171.0483
1.04671.32551.32370.81230.8118

1. 0488

1/16

1.0483

1.32891.3283

Exact Value

0.81191.05001.3333

1/4

n/8
0

niB
0

niB

1. 0487

Table 2: Maximal absolute errors of the computed

solutio~.

11

t = 0.1

t =0.5t = 1.0

1/2
1/4
1/8
1/1f)

0.01470.00700.00330.0016

0.06610.03190.01560.00770.13100.06490.03220.016053

Remarks:

1.

Those maximal absolute errors occur on thecomputed free boundary where the computed solution
vanishes.

2.

From Table 2 one can see that the observed rateof convergence in the maximum-nol'ill is one for
each

t.

Table 3: Position of the free boundary for
8 = 0 and computer time.
h1/2

t = 0.1

1/4
1/8

t = 1.0

Camp. Time2.0875t = 0.5
2.44232.88074.019sec.

2.09402.47132.48572.49295.963sec.
25.284sec.

1/16

2.09712.09862.93962.9696·2.9847

Exact Value

2.10002.50003.0000.~.--Remarks: 1.

254.324sec.

For other values of

-.

~-.

8

~-

the computed values

are nearly the same which means that the freeboundary is stably plotted.
2.

One can observe linear

convergen~efor theposition of the free boundary.
Example ?:one.

We take an example similar to the preceding

Only this timer*

is reduceo to a point, namely theeentre of the circular doma.ill:
gO = {(x, y)/O<p<l}•'1'11e exact Rnllltlon tn this case is chosen to be:so that the equation of the frpf> bnunoary is

p = t + 1.

In this case it only makes sense choosing boundary condition( 4) •

54

h = 11m,

Again we takethis time= O.

po ~ e .:: n/4,

w = 1

h = 6t/5.

andObviously

We perform the calculations only in the sector
T = 1,

up to time

O(T)

===

i. e., the final domain is given{(x, y) 10 < p < 2}.

Table 4: Maximal absolute errors of the computed solution
h

t

= 0.1

t = 0.5

t = 1.0

1/2

0.02990.13170.2516

1/4

0.0141

0.06430.1271

1/8

0.00670.03150.0635

1/16

0.00320.0156

0.0316

Remarks: 1.

The basic convergence properties shown in thepreceding example can be observed here also,
although convergence itself seems to become
slower.

2.

For the position of the free boundary thepreceding remark also applies.

Example 3.

edirections

III

We take data symmetric with respect to thein/4,

i . 0, 1,2, 3

4

and we solve the followin[

problem:
au = 6uatinnet) x [0, 1 ] ,plus Stefan condition onVor*

is given bys*(9) = 1rO

is given bysO(9) = 2 _ cos 49
10

aU(t

~) = 1av '-+u(t, x)

===

0

for
for

-+

x 6 r*-+x IS r(t)
55

and
and

t;

E [0, 1]t

6 [0, IJ.

r (t) •

= 11m and we take p = m (so that £ < 0 . 1m-1 ) .
We use again w = 1 but this time we take ~t = h/2.
For any quantity 1r, we denote by 1rh its approximation
for a certain value of h, and we define oirh by:
We definehe = o.

Table 5: P0sition of the free boundary for
t+

0.50

0.25

h

sh(O,t)

1/2
1/4
1/8
1/16
1/32

2.15002.1153
2.09602.08592.0807

Remark:aShsh(O,t)

347
193
101

2.31102.25852.23032.21552.2080

52

1.00

0.75
oSh

sh(O,t)sh(O,t)

aSh

350
181

2.55332.47682.43662.4158

765
402
208

92

2.4053105

aSh

2.4407
525
282
148
75

2.37502.34002.32192.3127657

Again the observed rate of convergence is one forthe position of the free boundary.
Table 6: Computed values of
t-+-

0.25

h

Uh (l,O,t)

1/2
1/4
1/8
1/16
1/32

0.87450.86500.85410.84690.8428Remarks: 1.

u

0.50

for

p =1

e = o.

and

0.75

1.00

oU h u h (l,O,t) QU h Uh (l,O,t)oU h Uh (l,O,t)

95
109
72
41

56

0.85220.84260.82910.820153

0.81510.85450.84240.82790.8183121
145

0.813096135900.86050.85450.84270.834660
118

50

0.830145~t/h

for theIt seems that the 'optima1 rate

purpose of economy is far less than

81

1/2.

we observed that for the same value of56

QU h .

hIndeed
the

l8a

\

\

\

/

\

/

\

\
\
\
\
I

,,".2.5~>:'1.t}

Figure 5= '.0'1

Figure 657

18b

Figure 758

approximate values obtained with

ht = hiSaremuch closer to the estimated solution (by extrapolation to the limit) than those in Table 5.
So, in practice

htshould be" small even for

w = 1.

2.

The number of operations necessary for solving
the system of linear equations iso(m 4 )

whereasthe composition of matrices at each time stepis onlyO(m 2 ).

So we suggest also the following

procedure:
Choosew = 0,

i.e., an explicit scheme, withlumped mass system [6J.

In this case there isno need to solve a system of linear equations at
each time step.

Although for explicit schemes

we must take= O(h 2 )ht

this could be advanta-

geous in some cases (for instance, whensufficiently large).

misThis statement is basedon the fact that the number of operations withthis modified scheme becomesO(m 5 ),o(m 4 )

instead of

and is also supported by the preceding

remark.
In Figure 5 we shaw the evolution of the free boundary forincreasing values oftcircles with centerO.

Example 4:andh · 1/32.

The dotted lines show

We have tested our algorithm to a whole domain

by solving a non-symmetric problem similar to the one of thepreceding example.

We have observed basically the same prop-

erties, although in this case we could not use too small meshsizes because the computer time increases sharply.
59

So we prefershowing the results obtained for the following problem in which
the domain evolves up to large values of- rO

is given by a set of

96

- r*

rO

So we consider

is defined by= 0.8.

pis the function such that

-

Uo

-

au(-+x 'av'd t E [0 , TJ.

t) = 1epoints equally spaced in

1. aand with polar radii. ranging from
Figure 6.

T:

and

1.4, as shown intoto be a polygon.

au .'
ap

LI -+

_0 = 1-+x E r

a, x En,= 0"*.

The aspect of the initial triangulation is shown in Figure

6 for

p

= 4,

m = 2.

m = 4 (so that= h,i. e. ,

E6t

we show in Figure
increases.

p

We actually calculate with-< o. 05).

= 0.1.

Definingh= O.8/p

=8

and
6t

we choose

We calculate up to timeT

= 10.0, and7 the evolution of the free boundary as time

In this way the computation lasts about

5.

3

minutes.

Concluding Remarks

Although we did not consider explicitly other cases, it
should be clear that everyone-phase Stefan problem in which thedomain evolves as a starshaped one can be treated in a similar
way to thatdes~ribed

in Sections 2 and 3.

In particular we

mention the case of bounded domains such as the regionnin Figure 8.

we canBy introducing appropriate

modif~cations,showntake into account the gradual transformation of edges of theapproximate free boundary into edges of the fixed external boundaryf,and then solve the problem similarly (see

60

[5J).

an = fur*

Figure 8
Unfortunately generalization of our algorithm to some casessuch as that of the two phase problem, seems not to be sostraightforward.

However, as we have mentioned other algorithms

based on the direct determination of the free boundary could be
efficiently used in such cases, provided that good methods foradjusting the spacial mesh step by step are available.

Indeed,one of the main features of our triangulation method is generating triangles whose angles remain reasonably bounded away fromzero, or, in other words, approximateiy equal.

In so doing thenumber of nodes of the mesh necessary to attain a given precision
can be minimized and computer time saved.
So, generally speaking, we think that our work could be astarting point for research on the application of automatic
discretization processes to the direct solution of two or threedimensional free boundary problems.

This is because we alRo

believe that with this work we have helped to ccntrad1ct thelong-prevailing opinion, that this approach is inadequate to thenumerical

solution of such problems.

6J

Acknowledgements
The author wishes to express his gratitude to Prof. H.
Fujita of the Department of Mathematics of the University ofTokyo for extending a warm welcome to him and for giving valuablesuggestions which were essential for the development of theseresults.
Thanks should be extended to Prof. T. Ushijima of the
Department of Information Mathematics of the University of ElectroCommunications, for his kind assistance in various academic
matters during the author's study in Japan.
This work was supported by the Japanese Ministry of Education(Mombusho) .

Refoerences[lJ

Ciavaldini, J. F.,,

Analyse Num~rique d'un Probleme de~,Stefan a Deux Phases par une Methode d'Elements Finis,
SIAM Journal on Numerical Analysis, Vol.12, 3 (1975), 464487.
[2J

Friedman, A.,

The Stefan Problem in Several Space Variables,

Trans. Amer, Math. Soc.,[3J

Fujita, H.,

132 (1968), 51-87,

Personal Communication,

Department of Mathe-

matics, University of Tokyo.
[4J

Jamet, ,P.,

Estimations d'Erreur pour des El~ments Finis

Droits Presque D~g'n'r~s,

R.A.I.R.O., Vol.lO, 3(1976),

43-fi1.
[5J

Jamet, P. and Bonnerot, R.,

Numerical Computation of theFree Boundary for the Two-dimensional Stefan Problem by
Time-Space Finite Elements,
2 (1977), 163-181.
62

J. of Compo Phys., Vol.25,

[6]Mori, M.,

Stability and Convergence of a Finite Elemertt

Method for Solving the Stefan Problem,

Publ. RIMS, Kyoto

Univ., 12 (1976), 539-563.

[7]

Mori, M.,

A Finite Element Method for Solving Moving

Boundary Problems,

Preprints of IFI'P Working Conference

on Modelling of Environmental Systems, Tokyo, April 1976,
167-171.

[8J

Ruas B. Santos, V.,

Sur l'Application de Quelques Erements

,

~

\

Finis non Conformes a la Resolution de Problemes Biharmoniques,,,

These, Universite Paris VI, April 1976.

[9J

Ruas B. Santos, V., Automatic Generation'of Triangular Finite
Element Meshes (to appear in Computer and Mathematics with
Applications).

63

Error Estimates for the Lumped Mass Approximationof the Heat Equation

By Teruo USHIJIMA
January 1979

Introduction.
In this paper error estimates for the lumped mass approximationof the inhomogeneous heat equation with zero Dirichlet boundary
condition are considered.

Generally the error is dominated by

the term caused by the lumping effect.

It is, however, bounded by

the term having the same order as in the case of consistent mass
approximation if the triangulation poseses the local symmetricproperty, which is a kind of regularity of triangulation, being
defined in §l of this paper.

Here we restrict our consideration

to a semi-discrete approximation scheme.

Namely we adopt a systemof ordinary differential equations as an approximate equation.
It may be well-known that the error is O(h) in the lumped mass
approximation for the heat equation.

And classical analysis teaches

us that the error is O(h2) in the usual difference scheme using
five point differnce approximation of A, whereas

this difference

scheme can be regarded as a special case of the lumped mass
approximation scheme, as was shown by courant [4].

Original

motivation of this work is to give a persuasible explanation ofthis discrepancy.
An idea of division of error term into two terms is essential.
The first term represents the error of Ritz projection, which
already appeared in the consistent mass approximation.
65

The second

term is proper to the lumped mass approximation, which is dominant,in general, being O(h) with respect to mesh size h.

In the case of

local symmetric triangulation, the estimate for the second term isimproved up to 0(h 2 ) , so that the first term is principal.
In §l, after giving formulation of the problem, we state main
result concerning L2-estimate in Theorem 2.

An abstract routine

to treat the error estimation is constructed in §2.
to prove Theorem 2.

§3 is devoted

In §4, LOO-estimates are driven as applicationsof our ahstract routine.
As for LOO-estimates, similar results to us have been obtainedby Tabata [8] (see also [9]).

He regards the lumped mass scheme

as a finite difference scheme defined on an irregular mesh.
In comparison with his method, this paper may be considered todescribe an operator theoretical approach of the lumping method.
The author expresses his sincere thanks to Professor Kikuchi
of Institute of Space and Aeronautical Science, University of Tokyofor his valuable advices and discussions during the preparation ofthis paper.

§l.

Statement of the problem and L2-estimates.

We consider the following continuous problem
au(t,x)_at-flu(t,x) - f(t,x),(~ )u(t,x)==0,u(O,x)

==

a(x),

(~).

(t,x)E(O,T]xn,(t , x) E (0 , T] x r ,x f S1,

where n is a hounded convex polygonal domain in n-dimensional
Euclidean space Rn with its boundary r, and f(t,x) anda(x) are given

known functions.
Let us denote the inner product of L2 (n) by ( , ) .

=

(u,v)

fnu(x)v(x)dx

forNamelyu,vEL 2 (n).

Define the inner product a( , ) of V = H~(n) by the formula:a(u,v) =

~

av)

(au

for u,v£V.

ai"j 'axj

j =1

Then there is a unique selfadjoint operator A in L2(n) with its domain
D(A) satisfying the following properties from (1.1) to (1.4):(1.1)

Av = -6.v(1.2)

for vED(A),(1.3)

a(u,v)(1.4)= (Au,v),

for uED(A), VEV.

As is well known,
D (A)= H~ (n)n H2 (n)holds since n is a convex polygonal domain (see Kadlec [6J).
As usual, 'we seek the solution u(t_,x) of (€) as the solutionof the following variational formulation (II) of(~),

or operational

formulation (E) of (t).
gt (u(t),v)+a(u(t),v) • (f(t),v),(II)u(t)£V,
u(O)

(E)

{

O<t~T, v£V,

O<t~T,= a.

d~~t) + Au(t)= f (t) ,t>O,uCO) = a.
Following Fujii [5J, let us formulate the lumped mass
approximation method.

Let positive numbers h be indices.

Assumethat there is a triangulation :1h = {T: simplex, diam(T) ~ h}of n,

67and that the family of triangulations {Jh:h>.O} is regular in thesense of Ciarlet-Raviart [2].

Namely there is a positive constanta satisfying that~h(T)/p(T)for T£~, h>O,a\'1here h (T) is a .. diam.eter of T and p (T) is a diameter of the inscribed
sphere of T.
Now we set our approximation space Vh of V as follows.
Vh = {vhEC(ITh ) : vhlr = 0,
vhl T £ PI

for any TE Jh}'

where C(n) denotes the space of continuous functions on IT, and vis
denotes the restriction of the function v to the set S, and Pk denotes
the totality of polynomials with degree at most k.
is said to be a nodal point.

A vertex of T£:Th

Let us count the interior and theboundary nodal points of n as b l , b 2 , ... , b N, and bN+ I , b N+ 2 ,b N+ M, respectively.
{w. :J

l~j~N+M}
W.

E

J

with the properties:
C(n h ),

Wj IT E PI
w.(b k )Jfor any T£

7h ,= o.ikol~j~N}

Clearly {w j :There exists uniquely a set of functions

forms a basis of Vh .

Let S. be the support ofJ

wj ' and let B j be the lumped mass region corresponding to the nodal
point b ..
J(For the definition of lumped mass region, see, forexample, p479 of Ushijima [10].)The characteristic function of B.

J

Let Vh be the linear space spanned by the set ofis denoted by w .•J.. l~j~N}.
J

functions {w.

The lumping operator Lh from Vh onto Vh

is defined as the linear mapping naturally generated by thecorrespondence: w.--+ w .•J

J

Namely we have68N

for .v h = J=
.L l V.w.£Vh..
J J
inve~seLet Kh be theof Lh .

The spaces Vh and Vh being closed

subspaces of L2 (n), we can define adjoint operators Lh*' and Kh*'from Vh onto Vh' and from Vh onto Vh' respectively, inner productsof these spaces being induced from L2 (n).

It holds thatand that(K h*vh,u- h ) = (vh,Khu
.
- h)

-h·
vh£V h , -uh£V

Let Ah be a linear operator acting on Vh defined by the formula:(AhUh,v h )= a(uh,v h ),uh,VhEV h ,We say that Ah is the Galerkin approximation of A.

Clearly Ah isa bounded positive self adjoint operator in the space Vh.

...

Definea bounded linear operator Ah acting on Vh by the formula:-

Ah = KhKh*A h ·Let us denote by TIh the linear interpolation
onto Vh"

operator from C(Q)Let P lh be the orthogonal projection from V onto Vh withrespect to the inner product a( , ), which is called Ritz projection,
sometimes, and le·t P Oh' and Jfh , be the orthogonal proj ection fromL2 (n) onto Vh' and Vh' respectively, with respect to the innerproduct ( , ).
Our approximate problem .Cnh) can be written as follows.
dat(Lhuh(t),Lhvh ) + a(uh(t),v h ) = (Lhfh(t),LhV h ),O<t~T,

uh(t) E Vh'uh(O)VhEV h ,0< t~T ,= ah ,

where fh(t) is a Vh-valued function, and a h is an element of Vh"
This problem (TI h ) is equivalent to the following Vh-valued evolution
69

equation (E h )(E h )~ uh(O) = a h ·du~~t)+

Ahuh(t) • fh(t),Defini tion 1.

A triangulationt~O,'7h is said to be locallysymmetric if S. is symmetric with respect to b. for any interiorJ

nodal point b jJ(l~j~N),- (S . - b .) = S. - b . ,J J
J J

namely, if it holds that

1 ~ j ~N •It is noted that wj(x+b j ) is an even function of x if
the triangulation is locally symmetric.
Two typical examples of locally symmetric triangulation in
2-dimensional case are illustraled in Fig. 1 and 2.

Fig. 1

Fig. 2

In both examples, we obtain the stiffness matrix corresponding tousual five points difference formula for -6.

As for mass matrices,

the matrix element of the point P is hk, where as that of Q, and R,
are 4/3hk, and 2/3hk, respectively.
Theorem 2.
where n~3.

Let g be a bounded convex polygonal domain in R nAssume that {Jh:h>O} is regular.
10

Then there is aconstant C, which may depend on T, with the property that the estimate

O~~~T II u h (t) -u(t) II L2 (0){m max= C h (O't~T<

II u(t)Hm(O)+O~t~T
II
max II dU II
max
d·t Hm(O)+O~t~T+ II a h - P 1 h a II L 2 (0) +0 ~~~ T II f h (t) - f ( t) II L 2 (0) }holds for the solution u(t) of (E) and the solution uh(t) of eEh)
.
l+m/2
au
. m/2 ), and that f(t) is so
provIded that u(t)EO(A), dU (t)ED(Asmooth that the quantity of the right hand side may be meaningful.
Here we can put m=2 if:rh is locally symmetric for any h, whereas

m=l in general case.

§2:

An abstract theory for error estimation.

Now we consider the following five conditions.
Condition I.

There is a Banach space X satisfying that

Condi tion II .

There is a Banach space Y contained in X andV, and a scalar function E(h) such thatII PI h v - v II X ~ E( h) II v II Y'
Condi tion ill.

v£ Y.

There is a Banach space Z contained in X; anda scalar function.o(h) such that

II Kh Kh Po h v - v II X ~ 0 (h) II v II Z'
111:

Condi tion N •-There is a constant MO satisfying

II e - tAhv h II X ~ MO II Vh II X'
Condition V.

v£ Z •

O~t~T,

vhEV h .

The solution u(t) of (E) satisfies

u(t) E YAD(A) , Au(t) E Z, andduet) £ Y
dt
for any tE: [0, T] •71

Under conditions from I to V, we have the following

Theorem 2.

estimate with a suitable constant C.

o~~~TII u h (t) -u(t) II X

~ C{e:(h) (o~~~TII u(t) II y+o~~~TII ~~II y)+6(h)o~~~TII Au(t) II z+1I ah-Plhall X+o~~~TII fh-fll Xl.
Proof.
~(2.1)

Ahr h = KhKh*POhAu,

since rh=A h
(2.2)

-1

d

POhAu.

~

Let gh=(ff(uh-rh)+Ah(uh-rh)' then we haveauaugh = (fh-f)-(P lhat - at) - (KhKh*POhAu-Au).

-

In fact, by (E h ) and (2.1),d~gh = fh-dfrh-Ahrhau= fh-Plha[ - KhKh*POhAu
auau= fh-Plhat + at + Au - f - KhKh*POhAu.
Therefore the 'function e h "'u h -r h is the solution of the following
evolution problem (2.3).
(2. 3)d[~~h + Ahe h "" gh

O~t~T,

eh(O) = a h - Plha.
By Duhamel's principle
(2.4)

-

we have~eh(t) = e-tAheh(O) +Ii e-(t-S)Ahgh(s) ds.

Substituting condition IT, TIl, V into (2.2), we obtain

II gh(t) II X ~ II fh-fll x+E(h) II ~~II y+6(h) II Aull z'
Therefore (2.4) and C"ondition N

imply the following estimate

II eh(t) II X ~ Moll ah-Pl h all X"'MOTCII fh-fll X+e:(h) II ~~II y+6(h) II Aull Z),

72

On the other hand, Condition n implies
II r h (t) -u(t) II X = II Plhu(t) -u(t) II X
~ E(h) II u(t) II y'Thus the conclusion of Theorem 3 follows from the triangular inequality:

Corollary 4.

Adding conditions from I to V, assume further

that Y is a closed subspace of Z.

Then we have the following-

estimate with E=max(E,&).

o~~~TII uh(t) -u(t) II X

~ C{E(h) (o~~~TII u(t) II z+o~~~TII :~II z+o~~~TII fll Z)+11 ah-Plhall x+o~~~TII fh-fll Xl.
Proof.

§3

Substitute Au=f-if into

the estimate in Theorem 3.

Proof of Theorem 2.

First we prepare some Propositions.
Proposition S~
Kh*P Oh can be considered as operators from LP(Q) to LP(Q) with theirnorms not greater than 1.
Proof.

See Proposition 2.1 of Ushijima [10].

Proposition 6.

For any p£[l,~], there is a constant L=L

independent of h such that

Proof.

See Proposition 2.2 of Ushijima [10].

Propos i tion 7.

. For aE:C (IT); we have

KhKh *p Oh a - Hha73p,n= j~l{mtB.) fs.(a(x)-a(bj))wj(x)dxw j },J

J

where m(Bj)=/BjdX.
Proof.

Noticing the orthogonality relation:

we have

On the otherhand,Nilha = J.Ela(b.)w.
=.
J J

=.~
1
I Sa
(b)
. t..l (B). w. d xw.,
J=

m j

j

J

J

J

where we use the equality:m(B J.) = IS w.(x)dx.
jJTherefore we have the conclusion.

The proof of Theorem 2 for. m=l is summarized in the following
Proposition.
Proposition 8.
MO=L 2 ,n'

Let X be L2(n), and let Y and Z be V.

Set

Then the conditions from I to N hold with suitable scalarfunctions E(h) and 6(h) which behave O(h) as h tends to O.
Proof.

Condition I

i~txivially valid.

By a standard error

estimation method for elliptic problem, we have

v€.v,

(3.1)

which assures condition IT (see Ciar1et [1]). Set vh=P1hv

for v£V.

By definition of P lh ,II vhll V ~ II vII V'
Since II vII V= II vvll L2 (n) is equivalent to II vII Hl (n)' (3.1) implies
(3.2)

(3.3)
74

KhlePOhVh-LhVh= (by Proposition 7)N
1=j~l{m(B.) IS.(Vh(X)-Vh(bj))Wj(x)dXWj}]
JN
1=j~l{m(Bj) ISj(X-bj'~Vh(X))nwj(x)dxwj}'

for vh(x) is piece wise linear
inner product of a,bdRn .

on Sj' where (c,b)ri means the Euchidean

Therefore

II Kh leP Oh v h -L h vh II 2
N

1

= J=
.t l m (B j ) lIs j (x-bJ"~VhCX))nw.(x)dxl
J

2

2 N
1
2
~ h j~lm(B.) IS. I~vhl dxmCS j )
J

J

222

~ h (n+l)

II VVh II L2 (fl)•Hence (3.2) implies( 3 • 4)

II Kh Ie P0 hv h - Lh v h II L2 (fl) ~ (n + 1 ) h II v II V •Using Proposition 5, we have

II Kh Kh Ie Po hv - v II L2 (fl)~ II (KhKh lep Oh -1) (v-v h ) II L2 (fl) + II KhKh lep Dh v h -vh II L2 (O)~ 2 II v - v hilL 2 (0) + II Kh Ie Po hv h - Lh v h II L2 (fl)~(by (3.2) and (3.4))~

Ch II v II V.

Thus Condition III is established.

Now let Ah=KhleAhKh.

a positive self adjoint operator acting in Vh.
lIe-tAhvhIlL2(O) ~II vh Il L2 (O)'

Since e -tAh=K h e-tAhL h'-

II e - tAhvh " L2 (fl)7S

Then Ah is

So we have~(Proposition 5)~

II Lhvhll L2 (n)~ L2,nll vhll L2(n)'Therfore Condition N holds with MO=L 2 ,n'The following Proposition completes the proof of Theorem 2 form=2.
Proposition 9.

Assume:1h be locally symmetric. Let X be L2 (n),12and let Y and Z be HO (n)n H (n). Then Condi tions IT and ill holdwith some scalar functions E(h) and o(h) which behave O(h 2 ) as htends to O.
Proof.

By a standard argument, we have

(3.5)

(see Ciarlet[l]).

This gives Condition IT.

To prove Condition ill,

first we note that IThv is well defined for v£H 2 (n), for H2cn)Cc(n)holds by Sobolev imbedding Theorem since we assumed n~3.
And we have
(3.6)

(see Ciarlet [1]).

Now we admit

the following estimate (3.7)

for a while.
(3.7)

II Kh Kh * P0 h v - v II L2 (n)<II Khll L(L 2 (n)) II Kh*POhv-LhIThvll L2 (n)+11 IThv-vll L2 (n)~(by Proposition 5, (3.7), and (3.0))~ Ch 2 11 vii H2 (n)'7bThis gives Condition Ill.
Now we proceed to establish (3.7).
(3.8)

Proposition 7 implies

II Kh*POhV-LhTIhVII2L2(0)N

12"'j~l m(B j ) I/Sj(u(X)~U(bj))Wj(x)dxl .

Let ~=x-bj'

Uj(~)=u(~+bj)' w(~)"'wj(~+bj)

and S=Sj-b j .

u(x)-u(b j ) = Uj(~)-Uj(O)= (VUj(O)'~)n+~l t(D2Uj(t~)~'~)ndt,
where D2Uj(~) the following n~n matrix:au
D2Uj(~) = (axjaxk (~+bj))l~j,k~n .
2Since S is symmetric with respect to the origin, we have
ISj(U(X)-U(bj))Wj(X)dX= IS(Uj(~)-Uj(O))W(~)d~... I s(/01t((D 2u j ) (t~)~,~)dt)w(~)d~1
2
2~ 10t(/slI D Uj(t~)II1Rnll ~111Rn d~)dt~(by Schwartz inequality and~ (diam(S) /2) 2

II ~II ~diam(s)/2)

Ilt II D2Uj (t~) II L2 (S) dt m(S) 1/2... h 2m(s) 1/2 101 t 1-n/2 11 D2u j (~) II L2 (tS) dt~ h 2m(S) 1/2 Ilt 1 - n / 2 11 D2Uj (~) II L2 (S)dt... h2m(S) 1/2 4-n
2 II D2U II L2(S.)'

J
Substituting this estimate into (3.8), we have

II Kh*POhV-LhIThVII £2(Q)
2 2 4 N m(S')
22'"' (4-n) h j~l m(B~) II Dull 1. 2 (Sj)<

2 2 411 22...< (n+l) 2 (4-n)
h
D ull L2(0)77Then we haveThis assures the validity of (3.7).

§4.

Loo-estimates for nonnegative triangulation.

Since we have prepared our abstract theory aiming a general00•

purpose routine, we can automatically obtain some L -est1matesgathering already es~ablished results if the triangulation isrestricted to be nonnegative.

Following Ciarlet-Raviart [3],the triangulation :1h is said to be nonnegative if and only if it
holds
(4.1)

a(w. ,w.)
J

1

~

f or

0,

'L'
l<'<N
1rJ,
=1=, l<'<N=J= +M•

In 2-dimensional problem, (4.1) is equivalet to the requirementthat all the angles of triangles T of ~h are not greater than n/2.
Following Fujii [5], this triangulation is said to be of acute type.
Throughout this §, triangulations are assumed to be nonnegative.
00

200

100

Let X=L (n), Y=W ' (OlAV and Z=W '(O)nV,
trivially satisfied.

Then Condition I is

Due to Ciarlet-Raviart [3], the following

estimate (4.2) holds.
(4.2)
This implies the condition IT with e(h)=O(h).

Moreover(4.4)

Proposition 7 assuresmean value theorem implies

II IIhv-vl! Loo(O) ~ hI! vII Wl,OO(n)'Condition III with o(h)=O(h) follows from (4.3) and (4.4).
Fujii [5], we have
(4.5)

II e- tAh vhll Loo(O) ~ II vhll Loo(n) ,78

Due to(See p487 of Ushijima [10].

It is noted that the nonnegativityof triangulation is essential to establish (4.2) and (4.5).)
Since e-tAh=Khe-tAhLht-

II e-tAhvhll Loo(O) ~ (Proposition 5 and (4.5))~

II Lhvh" Loo(fn= II vh" Loo(n)Thus Conditions W with MO=l holds good.

holds for vh£V h .

Hence Theorem 3 implies the following result.
Theorem 10.

Let 0 be a bounded convex polygonal domain in Rn

with arbitrary n.

Assume that the triangulation ~h is nonnegative

for any h.

Then there is a constant C with the property that theestimate

O~~~T II uh (t) -u(t) " L2 (n)~ C{h(o~~~TII u(t) " W2 '~(O) +o~~~TII ~~II W2 ,00(0)+O~~fTII f(t) II WI ,00(0))+11 ah-Plhall LW(O)+ 0 ~~~ T

II f h - f" L00 ( 0) }holds for the solution u(t) of (E) and the solution uh(t) of (E h )provided that U(t)£H~(0)nW2,00(0), AU(t)£H~(O)nWl,oo(O) andA result of O(h)-convergence of this problem is already obtainedby Tabata [8].
It is also possible to utilize the recent result of Loo-estimate
for the stationary problem due to Nitsche [7].

79

We restrict ourproblem to the 2-dimensional case, and assume further that our familyof triangulation {:rh:h>O} satisfies the inverse assumption.
Namely there is a positive constant v satisfying thath~vhCT)

C4.6)

Then Nitsche proves the following estimate C4.7):
C4 . 7)

II PI h v - v II L00 Cfl)

~ Ch 2 110 g h I II v II w2 ,00 Cfl)provided that v£VnW2,ooCfl).
Theorem 11.

Let fl be a bounded convex polygonal domain in R2.

Assume that the family of regular nonnegative triangulation{:Th : h>O} of n satisfies the inverse assumption C4.6), and that ~h
is locally symmetric for any h>O.

Then there is a constant C withthe property that the estimateo ~~~ T II u h (t) - u ( t) II L 2 (n)~ C{h 2 llog h I (o~~~TII u(t) II W2 ,oo(S"n +O~~~T II ffll W2 ,00(fl)+o~~~TII f(t) II W2JOD cn))+1I ah-Plhall LODcn)+o~~~TII fh-fll LOO(fl)}'holds for the solution u(t) of (E) and the solution uh(t) of (E h )provided that u(t), AuCt), }fCt)£H5(fl)nW2,OO(fl), O~t~T.
Proof.

Let X=LOD(n), and let Y and Z be D(A)nW2,00(fl)=H6cn)nW2,OOCfl).

Then Condition I is trivial.

Condition IT follows

from (4.7) and Condition N with MO=l is already shown.
Condition N is also proven by a step by step modification of theproof given in Proposition 9.

Namely instead of C3.6) and C3.7),

we have(4.8)

80

(4.9)
Corollary 5 implies the conclusion of Theorem.

References[1]

Ciarlet, P.G.,

The finite element method for elliptic

problems, (North Holland, Amsterdam, 1977).
[21

Ciarlet, P.G.,

Raviart, P.A.,

General Lagrange and Hermite

interpolation in Rn with applications to finite element methods,
Arch. Rational Mech. Anal., 46, 177-199(1972).
[3]

Ciarlet, P.G.,

Raviart, P.A.,

Maximum principle anduniform convergence for the finite element method, Comput.
Meth. Appl.
[4]

Engrg., ~, 17-31(1973).

Mech~

Courant, R.,

Variational methods for the solution ofproblems of equilibrium and vibrations, Bull. Amer. Math. Soc.,
49,1-23(1943).
[5]

Fujii, H.,

Some remarks on finite element analysis oftime-dependent field problems, IN Theory and practice infinite element structural analysis, 91-106(Proceedings of1973 Tokyo seminar of finite element analysis, Univ. Tokyo
Press, Tokyo, 1973).
[6]

Kadlec, J.,

On the regularity of the solution of the

Poisson problem on a domain with boundary locally similar tothe boundary of a covex open set (in Russian), Czech. Mat. J.,!!, 386-393(1964).
[7]

Nitsch, J.A.,

L~-convergence of finite element approximation,(Second conference on finite elements, Rennei, 1975).

[8]Tabata, M.,

Uniform convergence of the upwind finite

element approximation for semilinear parabolic problems, J. Math.
Kyoto Univ.,[9]~,

Tabata, M.,

327-351(1978).
LOO-analysis of the finite element method,

In Numerical analysis of evolution equations, 25-62(Lecture notes in numerical and applied analysis ~, Kinokuniya,
Tokyo, 1979).
[10]

Ushijima, T.,

On the uniform convergence for the lumped

mass approximation of the heat equation, Jour. Fac. Sci. Univ.
Tokyo, Sec lA, 24, 477-490(1977).

Department of Information Mathematics
The University of Electro-Communications
1-5-1, Chofugaoka, Chofu-shi
Tokyo, 184, Japan.

82