MEMOIRSOFNUMERICAL
MATHEMATICS
No.1

1974

-:<

N: ~~ k rq:
2289815
Memoirs of Numerical Mathematics No.1

1974

Masatake MORl
Yoshitane SHINOHARA
Hideo KAWARADA
and Makoto NATORI~

fkJf.#r~ff~fIJt

Contents
Tetsuro YAMAMOTO

r1J

Notes on the rate of convergence of theRayleigh-Ritz and Weinstein-Bazley methods

1

On the superiority of the trapezoidal rule forthe integration of periodic analytic functions

11

Numerical investigation of 1/2-subharmonic
solutions to Duffing's equation

21

On numerical solutions of Stefan Problem I.

43Managing Editors:
Minoru URABE

Department of Mathematics, Kyushu Univ.,
Fukuoka, Japan-8l2

Sin HITOTUMATU

Research Institute for Mathematical Sciences,
Kyoto Univ., Kyoto, Japan-606

with the cooperation of
Shigeru FURUYA
Department of Mathematics,. Univ. of Tokyo,
Tokyo, Japan-113

Hiroshi FUJITA

Department of Mathematics, Univ. of Tokyo,
Tokyo, Japan-ll3

Masaya YAMAGUTI

Department of Mathematics, Kyoto Univ.,
Ky .. to, Japan-606.

All correspondence should be sent to the Managing Editor.

i1Notes on the Rate of Convergence of the Rayleigh-Ritz
and Weinstein-Bazley Methods
Tetsuro

YAMAMOTO(Received 19 January, 1974)
1. Introduction
It is of considerably difficult problem to estimate therate of convergence of the Rayleigh-Ritz(RR) and theWeinstein-Bazley(WB) methods applied to general self-adjoint
operators.

For a certain class of operators, however,it canbe easily estimated by comparing the corresponding matrices.
The purpose of this note is to show this for the finite Hill
operator and give some observations foranother operator .

Let A be a self-adjoint operator in a real Hilbert spacehaving the inner product ( , ) and the norm Ilull=/(u,u).

Weassume that A is bounded below and that its spectrum consistsof isolated eigenvalues
plicity.

Al~A2~ .. 'each having finite

multi-

Further we assume that A can be represented as A(O)+A', where A(O) is a self-adjoint operator with known
crete eigenvalues

dis-1.1°) (AiO)~A2~ •.. ) and known orthonormal

eigenfunctions ulO). We denote by Aln)(Ain)~A~n)~ ••• ~A~n»(n)_t( (n)(n)
and r i - r il , ••• ,rin ) be the eigenvalues and the orthonormal eigenvectors of the nxn RR matrix Rn=(AulO),uJO»
based upon the test functions u l(0) , ••• ,u n(0) •
<A(n»
-n

and wln)=t(wl~), ••• ,wl~»
- 1 -

(n)

Ai

(n)

(AI

~

•••

stand for the eigenvalues

2and the orthonormal eigenvectors of the WB matrix Wn=(AiO)Oij)+([Pi,PjJ)-l respectively, where Pi= A,-luiO)(BaZley's special,

choice) and [u,vJ=(A'u,v) for U,VED(A'), the domain of A . It
is known that

if A(n)<A(O) A(n)<A(n+l)<A <A(n+l)<A(n) for,n - n+l' i ~ i- i- i= ieach i<n. u(n)=r(n)u (0) and u(n)=w(n)u(O) give the RR=ij =1 ij jij =1 ij j
and WB approximations for the orthonormal eigenfunction u iII

corresponding to Ai' respectively.
2. The Matrices Rn and Wn for a Certain Class of Operators
We first state the following theorem.
Theorem 1.

Let, (0)_ i+L
(0)
_
(1)
A u i - l 8ijUj(i-l,2, ... ,n)
j=l,
(n)
( )
for some scalar 8ij . Then Rn-Wn=(£ij ) where EiJ =0 1f i~n-Lor j<n-L.
-=
,

(0)

Proof. We first note that 8ij =(A u i

(0)

,u j

)=B ji and thatthe cond1t1oris (1) are equ1valent to, (0)_ i+L(0)
A u1 - l
81kuk(1=1,2, ... ,n)
k=1-Lor(0)_1+L
u1 - l
81k Pk'k=1-L(0) _i+L_ (0) (0) _ 1+L
Hence we have 0ij-(u i ,u j )- l8ik (Pk,u j )- l
8ik [Pk,Pj,J.
k=1-L
k=1-L
Th1s implies that, for each 1 such that i~n-L, the vector
i-L-l
n-i-L.
(6, ... ,0,81i_L, ... ,8i1+L'0, ... ,0) (or its transpose) g1ves the~

1-th row (or column) of the nxn matrix (bij)=([Pi,Pj])-l. In
order to determine the other elements bOT' we put bOT=BoT-E~~)(o,T=n-L+l, .• ;,n).

Then, subst1tuting the relat10n
- 2 -

3

into
n-L

n
(n) _
r
[Pa,PkJl3kT+ r[Pa,PkJ(l3kT-EkT )-oaT'
k=T-L
k=n-L+l
we now obtainn(n)_ T+L_(2)[p a' Pk J Ek T - [p a' Pk J 13 k T (a-I, 2 , .. • , n) .
k=n-L+lk=n+lClearl~ the system (2) has the unique solution 'Ej~)(n-L<j~n)r·r

for each T(n-L<T~n), and (2) is equivalent to
n
()
T +L
(3)
r
[Pa'PkJEk~ = - r[Pa,PkJ13kT (a=n-L+l, •.. ,n).
k=n-L+lk=n+lThis completes the proof.
Q.E.D.

3. The Rate of Convergence of the Methods Applied to
the Finite Hill Equation
In this section weshallapply Theorem 1 to estimate the rateof convergence of the RR and WB methods applied to the finite
Hill equation.

Let A be defined by
L

Au = -u"+( I 2c 2kcos 2kx)u, c 2L#o,k=lwhere tne Qomain of A consists of the functions u such thatu(O)=u(~)=o, u(x+~)=u(x), u' is absolutely continuous and u"E

L2(O,~).

We define the operators A(O) and A' as
A(O)u= -u"+au,

A'=A-A(O)L

respectively, where a is a constant such that a>M: I21c2kl.
k=lThenAiO)= 4i 2-a,uiO)=

J! sin 2ix (i~l).

Therefore the conditions (1) are satisfied with
L

l3 ij

= aO ij +. k=l
I c 2k (oi_J k-oi+j k)'
'
,
- 3 -

4

Lemma 1. Let A be the operator defined as above. Then there
exists a positive constant E independent on n such that IEly) I~E(n-L<i,j~n)for all n.

Proof. We rewrite (3) in the matrix form([ Pn-L+i,Pn-L+j]) «n)_
En-L+i:n-L+j)- -([Pn-L+i,Pn+j])C( 4)(i,j=1,2, •• ~ ,L),
where
Sn+l,n-L+l' "Sn+l,n
C =

=

o

o

As is easily seen, we have
[Pn-L+i,Pn-L+J] = 'lT2 J'IT ~(x)sin 2(n-L+i)x sin 2(n-L+j)xdx
o·-+- 'lTl J'IT ~(x)cos 2(i-j)x dx =Pij
.

0

and
1 J'IT ~(x)cos 2(j-i+L)x dx =qij-+-'IT

as n-+-oo , where ~(x)= (a+20

r c 2kcos 2kx)-.1 Further we can show

L

k=l
LM2
that, if a is sufficiently large (e.g.,a>M+2Ic2LI)' the LxL .

matrix P=(Pij) is strictly diagonally dom1nant, and nonsingular.
Therefore we conclude from (4) that(as n-+-oo(n)En-L+i,n-L+j)-+- ~P-1

QC, where Q=(qij) (i,j=1,2, ••• ,L).

Q.E.D.

From Lemma l and the relation W W(n)=A(n)W(n) we haven iii •- 4-

5such that n-mL +as nLemma 2. Let k>l be fixed and m be any positive integer
~

w(

n)

kj

0 ({

=

+

Then we have

00.

}m)

MIL

4(n-mL)2(j=n-L+l,n-L+2, •.• ,n),k(n)wn-ok , j

= O( n-[L J )where
M = M + 2(L-l)£.
Therefore, by the same way as in the previous paper[16J,
we can apply Wilkinson's result[15;pp.172-l73] to estimate therate of convergence of the RR and WB methods.

The results arestated without proofs as follows:
Theorem 2. Under the same assumption as in Lemma 2, we haveMIE

= 0( {}2m)

4 (n-mL) 2= O( n-2E~])!" )2 ~ 2~JCorollary 1. A(n) _ A(n) = o( MIL
k
k
n2

).

Theorem 3. Let (u~n)'Uk)~O for fixed k. Then, as n+~, wehave/lu(n) -k.

Uk II = o( {

"

MIL
}
4(n-mL)2- 5 -m) (k>l,-n-mL+~)

6

lI u (n)

U ~k

_ u

II = O( n~k-EkJ-IL) (k>O).

"[2£J

Corollary. 2.llu~n) - uk II = o( { M~}).

and r(n) so as to satisTake the vectors wen)ii

Theorem 4.
n

fy

-

l

O','I:=n-L+l
we have, as n-+oo,

w(n)r,(n»Ofor all n.
iO' 'i'l:Then... for fixed k,(k~l,n-mL-+oo)-(k>O) •

Corollary 3.

Remark. We can prove that~L

E(n)W(n)w(n»O for any
0''1:
iO' i'l:O','I:=n-L+I
nand i such that n>L and l~i~n. Therefore the assumption ofTheorem 4 will be satisfied for sufficiently large n.

4. Observation
Considerfor the Other Equation

an operator A defined by
Au= -u" + q(x)u,u(x)=u(-x), U(X)EL 2 (_oo,oo)

where q(x) is a polynomial of degree 2L such that q(x»x 2 -afor some constant a>O and q(-x)=q(x). (The case of q(x)=x 2 +ax 4 (O<a<l) has been treated by Bazley and Fox[2].) Then
A(O)u= -u"+(x 2 -a)u and A'=A-A(O) are suitable decomposition.
We have- 6 -7,(0) __' 4
I\ii- 3- a.,where,

m x 2 dm

H (x)=(-l) e - e
m .
dxm

_x 2

By the reccurence relation(the m-th Hermite polynomial).
2xHm = 2mHm_ l +Hm+l or

we see that the conditions (1) are satisfied.

Further wehave= 0(1)sincex2
I Hm (x) I~ e 2" /2 m+1 1IDi.

On the other hand, the relation (5) shows that Bij=O(i L ).
Hence, if n is large,

we may roughly conclude from (4) that

E(n)= O(nL ) (n-L+l<i,j<n).
Therefore, even if the methodsij=
are convergent, the rate of convergence in this case may berather slow, compared with the case of the finite Hill operator.

The more precise information will be given

elsewhere.

- 7 -

8

References[lJN.W.Bazley, Lower bounds for eigenvalues, J.Math.Mech.,lO(196l),pp.289-307.
[2JN.W.Bazley and D.W.Fox, Lower bounds for e genvalues ofsc~rgdinger's equation, Phys.Review, l24(196l),pp.

483-492.
[3JG.Birkhoff and G.Fix, Rayleigh-Ritz approximation bytrigonometric polynomials, Indian,· J.Math.,9(1967), pp.269-277.
[4J __________ , Accurate eigenvalue computations for elliptic
problems, in SIAM-AMS Proceedings, Numerical SO-'lution of Field Problems in Continum Physics,
G.Birkhoff and R.S.Varga, ed., A.M.S., 1970,pp.lll151.
[5JG.Fix, Order of convergence of the Rayleigh-Ritz andWeinstein-Bazley methods, Proc.Nat.Acad.Sci. U.S.A.,
6l(1968),pp.12l9-l223.
[6]D.W.Fox and W.C.Rheinboldt, Computational methods fordetermining lower bounds for eigenvalues of operators in Hilbert space, SIAM Review,8(1966),pp.427462.
[7]S.H.Gould, Variational Methods for Eigenvalue Problems,
2nd ed., Oxford University Press, 1966.
[8]A.J.Hoffman and H.W.Wielandt, The variation of thespectrum of a normal matrix, Duke Math.J.,20(1953),
pp.37-39.
- 8 -

9

[9J W.Magnus and S.Winkler, Hill's Equation, Interscience
Publishers,1966.
[lOJL.T.Poznyak,

Estimati~nof the rate of convergence ofa variant of the method of intermediate problems,
U.S.S.R.Comp.Math.Math.Phys.,8(1968),pp.246-260.
[llJ, The convergence of the Bazley-Fox Process
and the evaluation of this convergence, U.S.S.R.
Comp.Math.Math.Phys.,9(1969),pp.167-l84.

[12J, The rate of convergence of the Bazley-Fox
method of intermediate problems in the generalized
eigenvalue problem of the form AU=ACu, U.S.S.R.
Comp.Math.Math.Phys.,lO(1970), pp.56-74.

[13]A.Weinstein, Some numerical results in intermediateproblems for eigenvalues, in Numerical Solution ofPartial Differential Equations, J.H.Bramble,ed.,Academic Press, 1966.
[14JA.Weinstein and W.Stenger, Methods of Intermediate
Problems for Eigenvalues, Theory and Ramifications,
Academic Press, 1972.
[lS]J.H.Wilkinson, The Algebraic Eigenvalue Problem, Oxford
University Press,1965.
[16]T.Yamamoto, The Rayleigh-Ritz and Weinstein-Bazleymethods applied to a class of ordinary differential
equations of the second order, to appear.
[17]Urs.W.Hochstrasser, Orthogonal polynomials, in Handbookof Mathematical Functions, M.Abramowitz and I.A.Segun,

- 9 -

10

ed., DoverPublications,Inc.,1965.

Department of Mathematics
Faculty of Science
Ehime University
Matsuyama, Japan- 10 -JJ

On the Superiority of the Trapezoidal Rule for the
Integration of Periodic Analytic Functions
Masatake Mori
(Received 22 January, 1974)§lIntroduction
It is· well-known that the trapezoidal rule with anequal mesh size. yields a result with remarkably high accuracy
when applied to the integration of a periodic analytic
function over its period.

Frequently this rule is evencalled the "best" rule with rather vague definition of theword "best".

Several attempts seem to have been made toshow the "best" property of the trapezoidal rule (e.g.[ll)by comparing the results with other various well-known
formulas such as Gauss rules.
In the present paper we will define the best rule fornumerical integration of periodic analytic functions overits period in terms of the asymptotic decay rate of the
exponent of thecharacteristi~

function of the error [3],and show that the trapezoidal rule is the best under thisdefinition in the similar way as in the preceding paper [3].

§2

The complex integral representation of the error
Let f(x) be analytic over [a,b] and periodic with period

b-a.

I =

We consider the integral of f(x) over its period:

f:

(2.1)

f(x)dx- 1 -J?

We write an approximate formula for (2.1) as follows:(2.2)

where a k and Ak are the sampling points and the weights,
The error 6I n =I-I n of (2.2) can be expressed

respectively.

in the form of a complex integral:
Theorem 1

Let f(x) be analytic over [a,b] and periodic withperiod b-a.
61

Then

= Sba f(x)dx -

n

=~
t;1T1

1.

C

f Akf(a k)k=l'¥(z)f(z)dz - ~
t;1T1

J rC

(2.3)

where
'¥(z)

1m z > 0. { -7Ti

(2.4)

7Ti

00

'¥

1m z < 0n2: {I:

1n (z) =m=-oo k=l Ak ( z-ak-m(b-a)

1+ak+m(b-a)

) }

(2.5)

The path C of the integral (2.3) consists of two line segments
C1 and C3 parallel to the real axis as shown in Fig.l and istaken in such a way that there exists no singular pOint of.

fez) between the two line segments.
a+in
a:i

,

a-in':-----~~~-----: b-in'
C3

Fig.l

The path C of the integral (2.3)
- 2 -

13

Proof

The first term of (2.3) is evident since, if we

deform the path by letting n~O and n'~O, then we have
'Y ( z) f ( z ) d z =Jb f (x) dx

(Z.6)

a

It can be seen that the right hand side of (Z.5)
converges uniformly on compact sets not containing any ofthe poles and defines a meromorphic function.

Since fez)and 'Y n (z) are both periodic, 'Y n (z)f(z) is also periodic withperiod b-a so that we can join the lines CI and C3 with nochange of ~In by two line segments C2 and C4 parallel to theimaginary axis and form a rectangular contour C' withpositive orientation.

If a sampling point is located at

the end point x=a or x=b, Cz and C4 are deformed slightly
in such a way that they do not pass through the sampling
point and that they have the same shape(Fig.2).

Thenfrom the residue theorem,(2. 7)

and we have (2.3) from (2.6), which completes the proof.

Fig.Z

The path C'

We call
(2.8)

- 3 -

14
the aharaateristia funation of the error [3].

In case of

the trapezoidal rule with an equal mesh size of h=(b-a)/n,
we have

4>

(z)

n

={-lTi -lTcotKz =- 21T il-exp(-21Tiz!h)

1m z > 0
(2.9)

+2lTi
l-exp(+2lTiz/h)lTi -1m z < 0with the aid of the partial fraction expansion of cotangent
function.

Note that ~n(z) of the trapezoidal rule has nozeros on the finite z-plane.

Asymptotic decay rate of the exponent of the characteristic§3function and the best formula
When 1m z»h in the upper half-plane,~n(z)of (2.9)can be approximated by
I~

n

(z)1

~

2lTexp(-2lTy/h), y=Im z

and hence the quantity~

n2lT ' z=x+iy(3.2)

can be regarded to define the decay rate of the exponent ofI~n(z)1 in the upper half-plane.

Since .n~ (z)=~(z), wenwill consider the decay rate only in the upper half-plane.
It would intuitively be clear that ~n(z) of a useful formula
must converge to +lTi as 1m Z+±oo with all its derivatives andthat generally a formula having a larger decay rate of the
exponent of the characteristic function results in a smaller
- 4 -

'5
error when applied to the integration of a certain fixed
function.

And hence the decay rate can be considered to

be a criterion to define the best formula.

Accordingly

we define the average decay rate of the exponent of I~ n (z)/at a distance d from the real axis in the upper half-plane by

J ab+idI= - D7a

red)

.
a+,1d

d

y

logl~ (z)/dz ,

(3.3)

n

and define the asymptotic decay rate byr = lim
red)
oo

(3.4)

d-+

Now we call a formula having the largest asymptotic
decay rate the best formuZa.
Theorem 2

The trapezoidal rule with an equal mesh size isthe best formula for the integration of a periodic analytic
function over its period b-a among formulas having the same
number of sampling points on [a,b).
ProofConsider an integral

(3.5)

where

~

n

The path r is taken to(z) is defined by (2.5).

be a rectanglar contour with positive orientation with.comers
at b+id, a+id, a-id' and b-id' as shown in Fig.3.
a+id
Fig.3

The path rof (3.5)

rl

~:

b+id

b~i r 4

r 21,-_a___ ___a-id'

When ab-id'

16

sampling point lies on x=a or on x=b, r is modified as in
the proof of Theorem 1.
points.
TIi+~n

Now n is the number of the sampling

Note that n is nothing but the number of poles of(z) located insider.

Then from the principle of theargument we have
(3.6)

where n z is the number of zeros of TIi+~n(z) inside~~(z)/{TIi+~n(z)}

r.

Since

is periodic with period b-a, the integrals

along r 2 and r 4 cancel each other for any value of d and d'.
Since ~ n (z) must converge, on the other hand, to the constantfunction TIi as

d'~oo,

the integral along r3 can be made as

small as one wishes by letting d'~.

Hence we have

J

b+id

J

= 2;i

as

d'~~.

(3.7)

a+id

Since log{ni+~n(z)} is analytic along r l , we canreplace d/dz by partial derivative a/ax so that

Jb+id
2;i
Jb+id
= 2;i
Jb+id .J logl~n(z)ldz

J =

a+id

a+id

=

2;i

a+id

1 Jb+id_ ~~ya logl~n(z)ldz-'on.d rJV

(3.8)

a+lin view of Cauchy-Riemann equation.

- 6 -We know that J is real

17

from (3.6) and hence the first term of the right hand sideof (3.8) is zero.

Accordingly, by letting d~oo, we obtainan important relation between n-nz

and the asymptotic decayrate r:(3.9)

This shows that the asymptotic decay rate r cannot be larger
than ~~nd the maximum value ~is attained by suchv-ao-athat satisfies n =0 i.e. that has no zeros in thez 'We have, on thefinite z-p1ane, if such -~i-~n(z) exists.
-~i-~n(z)other hand, already seen that -vi-inez) of the trapezoidal
rule has no zeros in the finite z-plane, and hence thetrapezoidal rule is the best formula, in which r_2~n_2~-o::a-n'§4 Discussions
The characteristic function of the error of Simpson's
formulaIs = ~{f(a)+4f(a+h)+2f(a+2h)+4f(a+3h)+'"+f(b)}= jh{jf(a)+f(a+h)+f(a+2h)+"'+i f (b)}-~(2h) {~f(a)+f(a+2h)+f(a+4h)+" '+jf(b)} , h=b~a(4.1)

is given by
Im z > 0(4.2)
1m z < 0which can be approximated as
- 7 -

,

J8

I <f>sCz)I ~ ;nexpC-nlyl/h)

11m z I» h

(4.3)

Comparing this with (3.1) we see that the decay rate of the
exponent of (4.3) is half of that of (3.1).
situation also from the point of (3.9).

We see thisIndeed, if we put<f>s(z)=O in (4.2) for 1m z>O, we find an infinite array ofzeros of~s(z)in the upper half-plane arranged with

distance of 2h=2(b-a)/n at
z = i~log 3 + 2jh,j=O,±1,±2~ .•.

(4.4)

This means n z =n/2, and hence the asymptotic decay rate for
Simpson's formula is
r

nn
1T
= b-a
=n

(4.5)

from (3.9).
It should be remarked that there may happen to be acase where an approximate formula other than the trapezoidal
rule gives better result for certain analytic periodic
functions.

In fact Simpson's formula (4.1) would yield anexact result except the round-off error when applied to ameromorphic function having simple poles the location ofwhich coincides with that of zeros (4.4) of <f>s(z).

Sucha ~ase, however, is rather exceptional and , since theasymptotic approximation of

I <f>nCz)I such as (3.1) or (4.3)

becomes more precise as h is made smaller, the trapezoidal
rule would become superior to any other formulas when ahigh precision is required so that the number of the sampling- 8 -points is sufficiently large.
Finally we show a numerical example that shows thesuperiority of the trapezoidal rule.

To the integral

representation of Bessel function(4.6)

we applied the trapezoidal rule, Simpson's formula, Fi10n's
formula [2J and Legendre-Gauss rule.

The number of thesampling points of each formula is chosen to be 16.

Thecomputed absolute errors are as follows:Integration f~rmula

Absolute error

Trapezoidal rule

3.7X10- 19

Simpson's formula2.5x10- S

Fi10n's formula2.4x10- 3

Legendre-Gauss rule

6.1x10 -8 .

References
1.

P.J.Davis, On the numerical integration of periodic
analytic functions, On Numerical Approximation, R.Langer,ed., Univ. of Wisconsin Press, 1959, p.45-59.

2.

L.N.G.Ptlon, On a quadrature formula for trigonometricintegrals, Proc. Roy. Soc. Edinburgh, 49, 38-47 (1928).

3.

H.Takahasi and M.Mori, Error estimation in the numericalintegration of analytic functions, Report of the Computer
Centre, Univ. of Tokyo, 1, 41-108 (1970).
Research Institute for Mathematical Sciences,
Kyoto University, Kitasirakawa, Kyoto, Japan
- 9 -

21

Numerical Investigation of 1/2-Subharmonic Solutions
to Duffing's Equation
By

Yoshitane SHINOHARA(Received Jan. 1~, 1974, Revised Feb. IS, 1974)
1. Introduction
The present paper is concerned with 1/2-subharmonic
solutions to Duffing's equation
d2x
d tZ

(1 • i )dx
3+a . crt +bx+cx =e· cos wt,

where b>O, c>O and e>O.
The mathematical proof for the existence of a 1/2-subharmonic solution to (1.1) with damping absent, that is, with
a=O, has been given by C. T. Taam (6), T. Shimizu OJ and W.S.
Loud (2).

But, as far as the author is aware, good approxi-

mations to 1/2-subharmonic solutions are not yet found andthe mathematical proof for the existence of a 1/2-subharmonic
solution is not yet given to Duffing's equation (l.l) with
damping present, that is, with a*O.
In the present paper, good Galerkin approximations tothe subharmonic solutions in question will be given by meansof Galerkin's procedure established by M. Urabe(7J and of theauthor's program C5J of finding solutions to systems of nonlinear equations, and the mathematical guarantee for the
existence of the corresponding exact subharmonic solutions
will be also given by the use of Urabe' s existence theorem (7].
Replacing wt by t, we rewrite (l.l) as follows:
(1.2)

d.
d

J

x
t~

+ ~. dxw dt+ ~'(bx+cx3)= ~ cos t.
w
w

To a 1/2.subharmonic solution to (1~1), corresponds asolution to (1.2) of the form

-1-

22

(1.3)
Hence, replacing t by 2t in (1.2) and (1.3), we Can
reduce the problem to the one to find a solution of the form
(1 .4)

to the equation
(1.5)
where

The equation (1.5) can be rewritten in the following form

(1.6 )
where y=/Yx. The equation (1.5) oan be also rewritten in astand.ard form as follows:

(1.7)

dz
3
+0' . dT- +z+£· Z =oos

VT

,

where

Ii3 ·t=T , -t-x=z,

la =0' ,=£ , - -2= v

18

As for 1/2-subharmonio solutions to (1.5) with damping
present, C. Hayashi ell suggests some interesting properties
from his phase-plane analysis, but his assertions are madefrom rough approximate solutions of the formp
- 3·00S

2t

without givlng the proof for the existence of 1/2-subharmonic

-2-

?3

solutions.

Hence, his assertions do not seem to be of enoughmathematical confidence.
In the present paper, at first, we compute Galerkinapproximate solutions of the formby the use of the author's geometric method

C4J,l5)and then, starting from these rough approximate solutions, we computethe approximations to the solutions of the form (1.4) by theuse of Galerkin's procedure established by M. Urabe (8).
After having found an approximate solution by the above
procedure, it is necessary for the ~ompletion of the process
to verify the existence of an exact solution near the approximate solution obtained. and to find an error bound for theapproximate solution obtained.. Using the method developed by
M. Urabe (7) , we have checked the existence of the exact
subharmonic solutions and we have calculated the error bounds
for the approximate solutions. In Tables 2--8, c5 shows theerror bound obtained in this way. In Tables 1.... 9, the stabi,",lity of the subharmonic solutions are also shown.
The computations in the present paper have been carriedout by the use of TOSBAC 3400 at the Computing Institution of
the Reseach Institute for Mathematical Sciences, Kyoto
University.
The author expresses his hearty gratitude to Professor
M. Urabe for his constant advice and encouragement.
The author also wishes to acknowledge the .assistance of Miss.
T. Takahashi, Mr. Y. Kono and Mr. Y. Nakatani at the Computing Institution of the Research Institute for Mathematical
Sciences, Kyoto University.

-3-

24

2. Numerical Computation
In his paper (8), M. Urabe established a Galerkinprocedure to compute numerically the solutions of the form
(1.4) to the equation (1.5) by solving the so-called determining equation by the Newton method. In the present paper,
we shall use his procedure. According to Urabe's procedure,
we rewrite equation (1.5) in the first order system.

(2.1 ). =y,X(. =d Id t){

y=X(x,y,t)

where X(x,y,t)=-6X-YX 3 -ay+p.cos 2t.
For (2.1), a Galerkin approximation of order m is of the form

(2.2)

{

X(t)=c 1 +I. (c 2n ,sin nt +c 2n+1· cos nt),y(t)= t(-n.c 2n +1· sin nt +n·c 2n ·cos nt).

'>\. .. ,Hence, for (2.2), the determining equation for the Galerkinapproximations of order m can be reduced to the following
equationf 1 (C)=

~1f

51'"

X (x (t ) ,y (t ) , t 1 d t

=0,

o

(2.3)

1 (llt2
f 2n (C)=-;-) X(x(t),y(t),t)-sin nt dt+n 'c 2n=0,o
1.1t

f 2n +1 (C)=; ~ X(x(t) ,y(t) ,tl.. cos nt dt+n~c2n+l=0o(n=1,2,where C=(c 1 ,c 2 , - ~ -,m),,c2m,c2m+l).

In order to apply the Newton method to the system (2.3),
it is necessary to find the starting approximate s91utions to

-4-

25

the determining equations (2.3). For this purpose we.oonsider a Galerkin approximation (2.2) with m=2, that is,(2.4)
Then the determining equations for (2.4) are as follows:
)

)

2

2

f1 (01,02,03,04,05)=~01+y(01+1.50102+0.75°305+1.5°1C)
2 2 2+1.50104+1.50105+1.5020)04-0.750205)=0,)2

f2(01,02,03,04,05)=(a-1)02-~03+y(0.7502+30102

222+0~75°203+1.5°204+1.50205-3010205+3010304)=0,

(2.5)

3

2

f3(01,02,03,04,05)=~02+(a-l)03+y(0~7503+30103

22·2+0.750203+1.50304+1.50305+3010204+3010305)=0,
f 4 (01 ,0 2 ,0 3 ,04,0 5 )=( a-4)04-2Cl05+Y(0. 7504+30i04
222+1.50204+1.50304+0.750405+3010203)=0,

f5(01,02'03,04,05)=2a04+(a-4)05-P+Y(0.750~+30i05
2

2'

2

2

2

+1.50205-1.50102+1.50305+1.50103+0.750405 )=0.
For a system of the above determining equations witha=o • 1 , B =0. 3 , Y =0. 7, P=2. 0
that is,0=0.18), E=10).704, v=3.651,
we have oomputed the eleven solutions shown i:t:J. Table 1 1n thebounded region 10i'~) (i=1,2,3,4,5), by the use of thefortran program oonstruoted in the previous paper (5).
Now from the form of (1.5), we oan easily see,.that if

-5-

26

x(t) is a solution to (1.5), then x(t+rr), -x(t+{'31T/2)J and-X(t+{1T /2)1 are also solutions to (1.5), and that if x(t) is
a Galerkin approximation of 2~-periodic solution to (1.5),
then x(t+7r ), -X(t+('31T /2)) and -X(t+(1T /2)J are also Galerkinapproximations of 2TI-periodic solutions to (1.5) with thesame order as x(t).
For the solutions to (2.5) shown in Table 1, we readily
see that the Galerkin approximations x 2 (t), x 5 (t) andx6(t)corresponding to the 2nd, 5th and 6th solutions in Table 1
are equal respectively to xl (t+1T), -xlCt+('31T/2)1 and -x l Ct+(1T/2») , where x 1 (t) is the Galerkin approximation corresponding to the first solution in Table 1.
Likewise we readily see that the Galerkin approximations x4(t), x 7 (t) and x8 (t) corresponding to the 4th, 7th and. 8th
solutions in Table 1 are equal respectively to X3 (t+1T), -x 3 C
t+ (J1T /2 )land -x3 Ct+ (1T /2)1, where x3 (t) is the Galerkinapproximation corresponding to the '3rd solution in Table 1.
The Galerkin approximations x 9 (t), x l0 (t) and ll (t)corresponding to the 9th, 10th and 11th solutions in Table 1
will be supposed to be the Galerkin approximations ofharmonic solutions to (1.2).
Starting from the solutions to (2.5) shown in Table 1,by the use of the techniques described in the paper ( fa) wehave computed the Galerkin approximations of high order for
subharmonic solutions and, harmonic solutions. However, by
the reason mentioned above, we have not carried out the com:putations starting from the 2nd, 4th, 5th, 6th, 7th and 8th
solutions in Table 1. The results are shown in Tables 2 and3. In these tables, for each approximate solution is given
an error bound G such thatx

I

C\ x(t) _ X(t)\2 + ri(t) _ ~(t)A,l ]-r~ <5,where • =d/dt and, x (t) is an exact solution corresponding to

-6-

?7

the approximate solution ~(t).
3. The Effect of Nonlinear Term on the Subharmonic Solutions
In order to consider the subharmonic solutions to weakly
nonlinear equation (1.5), that is, equation (1.5) with IYK<l,we pursue a behavior of the stable subharmonic solution x 3 (t)obtained. in the previous section as Y ~ O.
From the last equation of (2.S) we have
1.S yc~(cS-cl)+l.SY

2

c~(cS+cl)=P-2 aC4+(4-8)cS- YC S (0.75 c~

2+3 c 1 +O.7S c 4)·
Taking account of the numerical results, we may consider thats is close to -P{3 which does not vanish forthe value of C
P=2.0.
(1-

Hence, we have:~ ) c~ + (1 + :~ ) oj =2P-4ac4+(8-28)C S

Jyc"Letting a~ 0 and. y-t 0, we see from the results of the numerical computations thatIc 1 /c S 1<1 and. IC4'« 1.
Hence, for small values of y and lal , we may estimate theamplitude

Jc~ +c~of the principal part of 1/2-subharmonic

solution as follows:(J .1 )

j c~ +c~

The value of the right side of (J.l) tends to infinity as Y
approaches zero. Hence, the amplitude of the subharmonic
solution increases as the value of y decreases.

-7-

28

Fig. 2 shows this phenomenon by waveforms of the stablesubharmonic solutionXJ(t) for various values of y •These solutions are shown in Tables 4",6.
These results show that 1/2-subharmonic solutions to
nuffing's equation do not exist in a neighborhood of periodic
solutions to the linearized equation. The fact tells us thatthe methods of perturbation which seek for solutions in aneighborhood of solutions to the linearized equation are notadequate for the computation of 1/2-subharmonic solutions
to nuffing's equation. This may be the main reason why 1/2subharmonic solutions to nuffing's equation have been nevercomputed so far.
M. Urabe (9] proves that for sufficiently small valuesof £ and 0(>0) there exists only the lj3-subharmonic
solution to (1.7), except for the solutions z(t) such thatI z (t)l-+ 00 or \dz /d t \ ~ 00 as £ ~ O.
Our results illustrate affirmatively his result.
From (J.l) the amplitude /c~ +cj of the principal partIof 1/2-subharmonic solutions tends to infinity with order y-~
as y approaches zero. But the value of IY. Ic~ +cj
isnearly constant as y approaches zero. The fact tells usthat it.is better to compute 1/2-subharmonic solutions to(1.6) instead. of (1.5) for small values of (l and y •Taking account of the fact, we have computed. the subharmonic
solutions to (1.6) with a=O.OOOOl, 8=0.3, Y=0.0000675 and
P=2.0. The results are shown in 'l'able 9.
Tables 2", 9 show that in the Fourier series of thesubharmonic solution in question, 'the t'irst five terms domi...,nate the remaining ones in strongly nonlinear cases, but thefirst seven terms dominate the remaining ones in weakly nonlinear cases. The fact tells us that one can know thequalitative character of l/c-subharmonic solutions to strongly nonlinear nuffing's equation by investigating the

-8-

29

character of the Galerkin approximations of the formx(t)=c 1 +c 2 ·sin t +c 3 ·cos t +c4·sin 2t +c 5 ·cos 2t,
but in weakly nonlinear cases one must take the more accurate
Galerkin approximations for the same purpose.
Remark. The effect of the damping term in the stable
1/2-subharmonic solution x3 (t) of the form
X3 (t )=c 1 + 'f c 2 . sin nt +
"":1
n

00

Ec2

'1'1=\

n +1' cos nt

00

appears strongly in the terms E c -sin nt"'-\ 2n•
When the damping is absent, we observe in Tables 7 and 8 that'f c -sin nt ~O,"".. , 2n
but when the d.amping is present t we observe in Tables 4 and 6
that
00

Ec

VI-' 2n. sin nt.

is not small and it increases as the value of y decreases.
In this case, however,!. ,c 2n+1 " cos nt
is effected very little by the damping term.
References(1)

Hayashi, C., Nonlinear oscillations in physical systems,
McGraw-Hill,New York, 1964.

(2)

Loud, W. S., Subharmonic solutions of second orderequations arising near harmonic solutions, J. Differential Equations Vol. 11(1972), 628-660.

[)]

Shimizu, T., Theory of nonlinear osci,llations, Baifukan,
Japan, 1965.

-9-

30

(4)

Shinohara, Y., A geometric method of numerical solutionof nonlinear equations and error estimation by Urabe's
proposition, Publ. RIMS. Kyoto Univ., Vol. 5(1969),1-9.

(5)Shinohara, Y., A geometric method of numerical solutionof nonlinear equations and its application to nonlinear
oscillations, PUbl. RIMS. Kyoto Univ., Vol. 8(1972), 13
-42.

(6]

Taam, C.T., The solutions of nonlinear differential
equationsN, Duke Mathematical Journal Vol. 24(1957),

511-519.
(7)

Urabe, M., Galerkin's procedure for nonlinear periodic
systems, Arch. Rational Mech. Anal 20(1965),120-152.

(8)

Urabe, M., Numerical investigation of subharmonic
solutions to Duffing's equation, Publ. RIMS. Kyoto Univ.
, Vol. 5(1969),72-112.

(9)Urabe, M., Subharmonic solutions to weakly nonlinear
Duffing's equation, On global theories of differential
equations, Lecture Note of RIMS. Kyoto Univ., Vol. 87,
April(1970), 61-91.
Yoshitane SHINOHARA
Department of Mathematics
Faculty of Engineering
Tokushima University
Minami Jyosanzima
Tokushima, Japan

-10-

31

01
°2
°3
°4
------------1: 0.0923726693 -0.2752659728 -0~6784945626 0.0597774582
2: 0.0923726692 0.2752659730 0.6784945625 0.0597774582

3:
4:0.1285965000 -0.1825706906 -0.7583156503 0.0632102469
0.1285964999 0.1825706907 0.7583156503 0.0632102469

5: -0.0923726692 0.6784945625 -0.2752659730 0.0597774582
6: -0.0923726693 -0.6784945626 0.2752659728 0.0597774582
7: -0.1285965000 0~75831565d5 -0~1825706905 0~0632102469
8: -0.1285964999 -0.7583156502
0.1825706909
0.0632102469
.
.
9:
10:
11:

0.0000000000 0.0000000000
0.0000000000 .0.0000000000
0.0000000000 0.0000000000

O~OOOOOOOOOO 0~0319575661
O~OOOOOOOOOO

0.8312796602
0.0000000000 0.5462865832

Os
1 : -0.6783580768
2: -0.6783580768
3: -0.6899371440
4: -0.6899371440
5: -0.6783580768
6: -0.6783580768
7: -0.6899371441
8: -0.6899371440
9: -0.5644062145
10: 2.760755463
11: -2.272539725
Table 1-11-

32

Periodic solutions to (1.5) with a=0.1,B=0.3,Y=0.7,P=2.0:
1) Xl (t)=0.0928604913-0.2729725113 sin t-0.6622053451 cos t+0.0592095043 sin 2t ~0.6775492475 cos 2t+0.0081759893 sin 3t -0.0104856750 cos 3t-0~0081146901 sin 4t -0.0070125277 cos 4t-0~0015293305 sin 5t -0.0075428283 cos 5t+0.0005026926 sin 6t -0.0017641013 cos 6t-0.0000735301 sin 7t -0.0001068038 cos 7t-0.0000611871 sin 8t -0.0000950798 cos 8t-0.0000011729 sin 9t -0~0000382199 cos 9t+0.0000008986 sin 10t-0.0000054232 cos lOt-0.0000010457 sin l1t-0.0000011578 cos llt-0~0000002502 sin 12t-0~0000006431 cos 12t+0~0000000026 sin 13t-0~0000001539 cos 13t-0~0000000098 sin 14t-0~0000000231 cos 14t-0~0000000066 sin 15t-0~0000000092 cos 15t-0.0000000009 sin 16t-0~0000000033 cos 16t-0.0000000001 sin 17t-0.0000000006 cos 17t-0.0000000001 sin 18t-0.000000000l cos 18t-0.0000000001 cos 19t,o=1.4xl0-7, Stability: unstable.
-0.7455563965 cos t
2)X3(t)=0.1325574730-0.1765810119 sin t,+0.0627505778 sin 2t -0.6894792968 cos 2t+0.0065184359 sin 3t -0.0122670665 cos 3t-0.0057077000 sin 4t -0.0096779756 cos 4t-0~0002877279 sin 5t -0.0085009813 cos 5t+0.0005555277 sin 6t -0.0019074934 cos 6t-0.0000489424 sin 7t -0 '. 0001903'122 cos 7t-0~0000306596 sin 8t -0 '.0001267016 cos 8t+0~0000061376 sin 9t -0.0000442724 cos 9t+0~0000015217 sin 10t-0.0000075120 cos lot-0~0000005539 sin 11t-0.0000020591 cos l1t-0.0000000308 sin 12t-0.0000008552 cos 12t+0~0000000388 sin 13t-0.0000002050 cos 13t-0.0000000016 sin 14t-0.0000000445 cos 14t-0.0000000021 sin 15t-0.0000000156 cos 15t+0.0000000004 sin 16t-0.0000000047 cos 16t+0.0000000001 sin 17t-O'.0000000011 cos 17t-0.0000000003 cos 18t-0.0000000001 cos 19t,o=1.8xl0- 7 , Stability: stable.
Table 2

-12-

33

2

'X"C1.)

Ol--~

____

ol--I-------rj=-~--~I

-/

-2
1t

"

t_

o

11

t-

2n

t-

Fig. 1

Waveforms of the 1/2-subharmonic solutions
given in Table 2.

-13-

34

Periodic solutions to (1.2) w1tha =0.1~ =O.l,y =0.7,P=2.0:
11 x9{t)=0.0000000001+0.0319635845 sin t -0.5644528144 cos t+0.0001659289 sin 3t -0.0008785586 cos lt+0.0000004612 sin 5t -0.0000014402 cos 5t+0.0000000011 sin 7t -0.0000000023 cos 7t.
6 =4. 7 10- 9 , Stability: stable.
l(2') x 10 ( t )=0.0000000002+0.808599 8 729 sin t +2.6920234635 cos t+0.1077254 917 sin 3t +0.0940200536 cos 3t+0.0067635484 sin 5t +0.0009948792 cos 5t+0.0002990601 sin 7t -0.0001350083 cos 7t+0.0000085908 sin 9t -0~0000132058 Cos 9t+0.0000000051 sin llt-0.0000007564 cos lit-0 '.0000000194 sin 13t-0~0000000307 cos 13t-0 '. 0000000016 sin 15t-0.0000000007 Cos 15t-0.0000000001 sin 17t,·-8
6 =9.11(10, Stability: stable.
3) x,,(t)=0.0000000001+0.530508 56 0 8 sin t -2.2314158468 cos t-0.0000000001 sin 2t+0.0462811064 sin 3t -0.0527460092 cos 3t+0.0019772093 sin 5t -0'.0007726566 cos 5t+0 '.0000635606 sin 7t +0.0000067958 cos 7t+0.0000016033 sin 9t +0'.0000010635 cos 9t+0.0000000281 sin 11 t+O '.0000000506 cos llt+0.0000000017 cos 13t,
6 =3 .4/C 10-8 , Stability: unstable.
Table 3-14-

35

Periodic solutions to (1.5) witha =0.001,s=0.3,P=2.0:-11.4799570291 cos t
X3(t)=0.2930123755-1.2590140124 sin t+0.0341190409 sin 2t -0.6997695126 cos 2t-0.1130152934 sin 3t -0.3336762141 cos 3t-0.0052523794 sin 4t -0.0329736966 cos 4t-0.0054219996 sin 5t -0.0101881061 cos 5t-0.0005610178 sin 6t -0.0014392693 cos 6t-0.0002349807 sin 7t -0.0003180846 cos 7t-0.0000342579 sin 8t -0.0000535047 cos 8t-0.0000098891 sin 9t -0.0000100401 cos 9t-0.0000017126 sin 10t-0~0000018143 cos lOt-0.0000004099 sin llt-0.0000003113 cos llt-0.0000000771 sin 12t-0~0000000569 cos 12t-0.0000000167 sin 13t-0~0000000092 cos 13t-0.0000000033 sin 14t-0~0000000016 cos 14t-0.0000000007 sin 15t-0.000000OO02 cos 15t-0.0000000001 sin 16t,
0=5.6 x10- 6 , Stability: stable, where y=0.00675)(1.0.
-5.12751979 4 2 cos t
X3(t)=0.2961798507-0.1083332858 sin t.
+0.0071156489 sin 2t -0.6962044099 cos 2t-0.0093081401 sin 3t -0.1496922067 cos 3t-0.0009722017 sin 4t -0.0322234685 cos 4t-0.0004426696 sin 5t -0.0070009504 cos 5t-0.0000999061 sin 6t -0.0015101723 cos 6t-0.0000260509 sin 7t -0.0003249848 cos 7t-0.0000063492 sin 8t -0.0000699167 cos 8t-0.0000015387 sin 9t -0.0000150223 cos 9t-0.0000003674 sin 10t-0.00oo032274 cos lot-0.0000000869 sin llt-O.0000006931 cos 11 t
-0.0000000204 sin 12t-0.0000001488 Cos 12t-0.00000000J..J.7 sin 13t-0.0000000319 Cos 13t-0.0000000011 sin 14t-0.0000000069 cos 14t-0.0000000003 sin 15t-0.0000000015 cos 15t-0.0000000001 sin 16t-0.0000000003 cos 16t-0.0000000001 cos 17t,
0= 5 .. 3 IC 1 0 -7 , Stability: stable, where y =0.00675 l( 5.0.
Table '4

-15-

36

Periodic solutions to (1.5) with a=0.001,8=0.3,P=2.0:
X3(t)=0.2914436944-0.0378554054 sin t-3.5922663954 cos t+0.0037432939 sin 2t -0~7002181244 Cos 2t-0.0030002747 sin 3t -0.0995051261 cos 3t-0.0004751640 sin 4t -0.0313425705 cos 4t-0.0001337269 sin 5t -0.0065898220 Cos 5t-0.0000441812 sin 6t -0.0015028645 cos 6t-0.0000107546 sin 7t -0.0003512176 cos 7t-0.0000028383 sin 8t -0.0000793848 cos 8t-0.0000007390 sin 9t -0.0000181357 cos 9t-0.0000001854 sin 10t-0~0000041449 Cos lOt-0~0000000467 sin l1t-0~0000009443 cos l1t-0.0000000116 sin 12t-0~0000002154 cos 12t-0.0000000029 sin 13t-0~0000000491 cos 13t-0.0000000007 sin 14t-0~0000000112 cos 14t-0.0000000002 sin 15t-0~0000000026 cos 15t-0~0000000006 cos 16t-0.0000000001 cos 17t,
6=2.2)(10- 7 , Stability: stable, where y=0.00675l(10.0.
x3(t)=0.1729240185-0.0013345231 sin t-0.8434488108 cos t+0.0006816817 sin 2t -0~7085991147 cos 2t+0.0000573890 sin 3t -0.0152417827 cos 3t-0.0000449412 sin 4t -0.0121875975 cos 4t+0.0000049194 sin 5t -0.0094382956 cos 5t+0.0000061983 sin 6t -0~0021498124 cos 6t-0.0000002845 sin 7t -0~0002818746 cos 7t-0.D000000987 sin 8t -0.0001522820 cos 8t+0~0000001164 sin 9t -0.00005206 6 9 cos 9t+0~0000000236 sin 10t-0~0000105388 cos lot-0.0000000010 sin llt-0.000002890 4 cos llt+0.0000000014 sin 12t-0.0000010626 cos 12t+0~0000000008 sin lJt-0.00000028J3 cos lJt+0.0000000001 sin 14t-0.0000000699 cos 14t-0.0000000218 cos 15t-0.0000000066 cos 16t-0.0000000017 cos 17t-0.0000000005 cos 18t-0.0000000001 Cos 19t,
86=7.0",10- , Stability: stable, wherey=0.00675x 100.0.
Table 5-16-

37

Periodic solutions to (1.5) with CI. =0.001 J3 =0.3 ,P=2.0:-0.5514182480 cos t
X3(t}=0.1113354509-0.0010744200 sin t+0.0005308837 sin 2t -0.6737278095 cos 2t+0.0000737828 sin 3t -0.0136047713 cos' 3t-0.0000344817 sin 4t -0.0055800058 cos 4t-0.0000018926 sin 5t -0.0076511441 cos 5t+0.0000064399 sin 6t -0.0023812607 cos 6t+0.0000002092 sin 7t -0.0001772272 cos 7t-0.0000002303 sin 8t -0.0000978628 cos 8t+0.0000000720 sin 9t -0.0000496638 cos 9t+0.0000000385 sin 10t-0.0000107395 cos lOt-0.0000000005 sin llt-0.0000018419 cos llt-0.0000000001 sin 12t-0.0000008334 cos 12t+0.0000000008 sin 13t-0.0000002810 cos 13t+0.0000000002 sin 14t-0.0000000598 cos 14t-0.0000000156 cos 15t-0.0000000057 cos 16t-0.0000000016 cos 17t-0.0000000004 cos 18t-0.0000000001 cos 19t?
6=7. 6 )C 10 -8 ,Stabili ty: stable, where y=O. 00675 K 140.0.
-0.2491740565 cos
X3(t}=0.0341643187-0.0014957675 sin t+0.0003984503 sin 2t -0.6185496392 cos+0.0000860099 sin 3t -0.0092860164 cos-0.0000264466 sin 4t -0.0009737959 cos-0.0000170626 sin 5t -0.0036648354 cos+0.0000047346 sin 6t -0.0021630661 cos+0.0000006181 sin 7t -0.0000844859 cos-0.0000003103 sin 8t -0.0000226179 cos-0.0000000801 sin 9t -0.0000241399 cos+0.0000000300 sin 10t-0.0000084863 cos+0.0000000020 sin l1t-0.0000005800 cos-0.0000000021 sin 12t-0.0000002061 cos-0.0000000003 sin 13t-0.0000001293 cos+0.0000000002 sin 14t-0.0000000348 Cos-0.0000000037 cos-0.0000000014 cos-0.0000000006 cos-0.0000000002 cos

t
2t
3t
4t
5t
6t
7t
8t
9t
lOt
11 t
12t
13t
14t
15t
16t
17t
18t,

\5 =1.2)( 10- 7 , Stability: stable, where y =0.00675)( 180.0.

Table 6

-17-

38
x

·.1

...
a=O.OOl/3=0.3 '.
P=2.0.

y=o.oo675&... 10I. 0 . -_ _ . ......
_ ••• .._ .._ .•.• _ _~.".H. ______ ._.;-___.. --. ..-,.-.

tr. F1g. 2

39

Periodic solutions to (1.5) with a=0,B=0.3,P=2.0:

X3 (t)=0.3002627 61 J+0.0000000306 sin t-11.5505768414 cos t-0.0000000011 sin 2t -0.6923730347 cos 2t+0.0000000027 sin 3t -0.3526089531 cos 3t+0.0000000001 sin 4t -0.0328902369 cos 4t+0.0000000001 sin 5t -0.0116843449 cos 5t-0.0015252088 cos 6t-0.0004043236 cos 7t-0.0000630335 cos 8t-0.0000144947 cos 9t-0.0000024908 cos lot-0.0000005309 cos lit-0.0000000963 cos 12t-0.0000000197 cos 13t-0.0000000037 cos 14t-0.0000000007 cos 15t-0.0000000001 cos 16t,
0= 5 .0)( 10 -6 , Stability: neutral, where y =0.00 675)( 1.0.

X3(t)=0.2964666764+0.0000000032 sin t -5.1288186026 cos t-0.0000000003 sin 2t -0.6959386175 cos 2t+0.0000000002 sin 3t -0.1500079715 cos 3t-0.0322212179 cos 4t-0.0070228217 cos 5t-0.0015132373 cos 6t-0.0003263119 cos 7t-0.0000702336 cos .8t-0.0000151121 cos 9t-0.0000032506 cos lot-0.0000006991 cos l1t-0.0000001503 cos 12t-0.0000000323 cos 13t-0.0000000070 cos 14t-0.0000000015 cos 15t-0.0000000003 cos 16t-0.0000000001 cos 17t,o =.5.1)( 10- 7 t Stability: neutral, where y=0.00675 x 5.0.

Table 7-19-

40

Periodic solutions to (1.5) with a=0,8=0.3,P=2.0:-0.5514227536 cos t-0.6737293984 cos 2t-0.0136049668 cos 3t-0.0055800896 cos 4t-0.0076512213 cos 5t-0.0023812868 cos 6t-0.0001772325 cos 7t-0.0000978644 Cos 8t-0.0000496646 cos 9t-0.0000107398 cos lot-0.0000018420 cos 11 t
-0.0000008334 cos 12t-0.0000002811 cos 13t-0.0000000598 cos 14t-0.0000000156 cos 15t-0.0000000057 cos 16t-0.0000000016 cos 17t.-0.0000000004 cos 18t-0.0000000001 cos 19 t ,. &=7. 6 x 10 -8 ,Stability: neutral, where y=O. 00675 x 140.0 •X3{t)=0.1113376437+0.0000000002 sin t-0.2492037267 cos t-0.6185554553 cos 2t-0.0000000001 sin 2t -0.0092870763 cos 3t-0.0009740593 cos .4t-0.0036653121 cos 5t-0.0021631498 cos 6t-0.0000845008 cos 7t-0.0000226241 cos 8t-0.0000241435 cos 9t-0.0000084870 cos lot-0.0000005801 cos l1t-0.0000002062 cos 12t-0.0000001294 cos 13t-0.0000000348 cos 14t-0.0000000037 cos .l.5t-0.0000000014 cos 16t-0.0000000006 cos 17t-0.0000000002 Cos 18t,

X3{t)=0.0341731575+0.0000000007 sin t

6 =1.2)( 10-7 ,Stability: neutral, where 'l-O.0067.5x 180.0.

Table 8-20-

41

Periodic solutions to (1.6) with a=10- 5 ,B=0.3,
y=O • 00 6 75 x 10 -2 , P= 2.0 :Y1 (t)= .ry,x1 (t)-0.7414373463 cos t=0.0005363391-0.5948614354 sin t+0.0002773728 sin 2t -0.0069298418 cos 2t-0.0263042508 sin 3t +0.0129549704 cos 3t-0.0003350052 sin 4t -0.0000915269 cos 4t+0.0002071635 sin 5t +0.0008519793 cos 5t-0.0000074282 sin 6t +0.0000141233 cos 6t+0.0000261856 sin 7t -0.0000004987 cos 7t+0.0000004872 sin 8t +0'.0000004127 cos 8t+0.0000001559 sin 9t -0.0000007657 cos 9t+0.0000000190 sin 10t~0~0000000144 cos lOt-0.0000000212 sin llt-0~0000000095 Cos 11 t
-0~0000000b04 sin 12t-0~0000000008 Cos 12t-0.0000000004 sin l3t+0.0000000006 cos 13t,
Stability: unstable.

Y3 (t )=IY. x3 (t)=0.0024156208-0.1036331599 sin t-0~9449355963 cos t. +0.0002773708 sin 2t -0~0057425427 cos .2t-0.0094390560 sin 3t -0~0277658105 Cos 3t-0.0000433228 sin 4t ~0.0002723404 cos 4t-0.0004562977 sin 5t -0~0007517027 cos 5t-0.0000047033 sin 6t -0~0000118494 cos 6t-0.0000182337 sin 7t -0~0000190670 cos 7t'-0.0000002870 sin 8t -0.0000004228 cos .8t-0.0000006579 sin 9t -0~0000004411 cos 9t-0.0000000139 sin 10t-0.0000000131 Cos lOt-0~0000000221 sin llt-0.0000000087 cos 11 t
-0~00ooooo006 sin l2t-0~0000000004 cos l2t-0.0000000007 sin lJt-O.OOOOOOOOOl Cos 13t,
Stability: stable.
Table 9-21-

43

On Numerical Solutions of Stefan Problem I.
Hideo KAWARADA and Makoto NATORI(Received Jan. 19, 1974; Revised Feb. 15, 1974)§l.

Introduction
Various works have been published on Stefan-type freeboundary problems for heat equations.(1)~(4)Stefan problem is to seeku=u(x,t)ands=s(t)

which

satisfy the following equations (FBP.I):( 1.1)

u t = u xx ' O<x<s(t),(1.2). u(O,t) = f(t) > 0 ,(1.3)

u(s(t),t) = 0 ,(1.4)

u(x,O) =(1.5a)

ds(t)/dt = -ux(s(t) ,t),(1.5b)

s(O)

where

O<T<+oo.

~(x)

O<t~T
0~t2'l'

O~t~T

~

,
,

,

0 , O<x<b ,~(b)=O, b>O ,

O<t<T,=b ,The caseb=O

will be discussed in the

succeeding paper.
The assumption required on the Stefan data is as follows:
Assumption (A):
(ii)
(1.6)

(i)

f

~

and

are nonnegative and continuous.

there exists a positive constanto < ~(x) < D(b-x)

,

for

D

such that0 < x < b .

====Recently, NOgi(5) presented a difference scheme for FBP.I
with a term of artificial heat flow and proved the unique existence and the convergence of the solution of his scheme.
We propose a new difference scheme by the penalty method.
Application of the penalty method to initial-boundary value
problems in a noncylindrical domain (i.e., moving boundary
prQblems) was done earlier by FUjita.(7)In order to apply the penalty method to ~tefan problems,
an integral representation of the free boundary-1-

x=s(t)

is44
indispensable.
Therefore, let us use the following reformulation(6)of Stefan condition (1.5):

Io

s (t)

s(t)2 = F(t) - 2(1.7a)~u(~,t)d~ ,O<t<T ,
= =where

21b~q,(~)d~F(t) = b 2 + 21tf(T)dT +(1.7b)o0

We shall denote the differential system(1.l)~(1.4)

and (1.7) by FBP.K.
From (1.7) we sees(t)2 =< F(t)
since

f(t»O
(O<t<T)
-=e::

,
and

iii::

u(x,t»O(O<x<s(t),0<t<T),which===e::: -=_ _follows from maximum principle.
the domain

Choosing

X

is included in

O<x<s(t), O<t<Tsuch that

X>F(t),=

O=[O,Xjx[O,T].

We consider the following equations with the penalty term(?)
in the region

O=[O,X]x[O,T],(1.8)

u t = u xx - KX(x,t)u(1.9)

u (0, t) = f(t)(1.10)u(X,t) = 0(l.ll)

u(x,O) = t~X)

,

,

in

0

X(x,t) =

O<t<T

== =

,

{ ,
1

s(t)2 = F (t) -(1.13 )

,

O<t<T
= =

,

O<x<b
b<x~X

(1.12)

0,

O<x<s= =:z (t) , O<t<T= =,s(t)<x<X,
O<t<T
=
= =

,

21X~u(~,t)d~ ,
0

where

K

is a positive number.

- 2 -The system

O<t<T
- ==

,

{l.8}~(1.13}

45
approximates well FBP.E if

K

is sufficiently large.

A

merit to solve the above-mentioned system in place of FBP.E

is that the calculation is markedly simplified by the replacement of the boundary condition (1.3) on
condition (1.10) on
Now,

~'le

u

by the boundary

x=X (which is independent of t).

introduce the following difference scheme (FBP .JI!) for

the system(1.14)

x=s(t)(1.8)~(1.13):

m,n

-

u

k

u
- 2u
+ u
m,n-I
m+l,n-l _ KXum,n-l = m-l,n-l
h2
m,n m,nl<m<M-l, l<n<N= ===== ===(1.15)(1.16 )
(1.17)

u O,n = f n ,u M,n = 0 ,O<n<Nu m,O =f:

Xm,n ={ ,- 2h2 m=ll m u m,n ,m

0(1.18)

O~n<N

I

,

l<m~.rb/h],[b/h]+I<m<M

,

l~[s

- -

n !h],O<n<N= =[Sn!h]+l~m.91-l

M-l, O<n~N(1.19a)

s

2
F
n = n

(1.l9b)Fn-l[bi h ]
b 2 + 2k I f n + 2h2=
m<l>m •
n
m=l
n=O

O<n<N
= =

where
h = mesh size in the x direction,k = mesh size in the t direction,
T

= Nk,

X = Mh,um,n = u(rnh,nk) ,

- 3 -

I

46
sn = s (nk..1,<Pm = <P (mh)fn

= flnk).

and raj denotes the greatest integer not" exceedinga.

The objective of the present paper is to construct anapproximate solution of FBP.TIIby iteration.
In the succeeding paper, we will discuss
existence of the solution of FBP.TII;

(ii)

(i)

the uniquethe convergence ofthe solution of FBP.DI to the solution of the difference scheme
for FBP. IT

asK -+ 00;(iii)

ence solution of FBP.IT
system FBP. IT

asthe convergence of the differ-to the solution of the differentialh, k -+ O.

In §3, we propose our algorithm.
numerical example.

In §4, we show a§2 is devoted to the preparation.

.,.. 4 -47§2.

Some preliminaries
Asa preparation we state the following lemmas.

sn

In FBP.ill, the free boundary
But ifis determined by (1.19).

s

is previously given and K is sufficiently large,n
we can consider (1.14)~(1.18) as the system which approximatesa moving boundary problem.
system by MBP.I.

Lemma 1.
and~We shall denote the above-mentioned

Let the solution of MBP.I be

In MBP.I, suppose thatO<k/h2~1/2

satisfy the assumption (A).

Then

(K)

um,n·

and that

u(K)
m,n

f

satisfies(2.1 )
where

C = max{ max f n ,
max
~}
O<n<N
l<m~Lb/h]m= =Proof.

From (1.14) we have

(2.2)

u(K) = _ _l--:=-c-{).u (K)+ (1-2),) u (K)+ Au (K)}m,nl+kKX(K)m-1.n-1m,n~1
m+1,n-1m,nwhere

By the assumption we see that the coefficients in the bracesof the right hand side of (2.2) are all positive and the sumof them equals
0 =< fn =<unity.

Making use of the facts

max f
n
O<n<N
= =

and
0 =< ~m =<max~l<~[b/hJ m

- 5-

48
we easily obtain (2.1).

Lemma 2.

and

u"
m,n
O<s'<s"<X=. n= n=s' and S" be given and u'm,nnn
be the corresponding solutions. Suppose that
In MBP.I, letin

O<n<N.

Then we haveu'm,n =< u"m,n(2.3)

Proof.

Put.

X~,n =!10 ,l<m<[s~/h],O~n~N[s~/h]+l<m~-l,,

O<n<N ,

and

X;,n

=f: :

l<m<[s~/h],

O~n~N

Ls~/h]+l~m~M-l,

We can prove by induction.

,

O<n<N •First we claim thatu'm,l'"< u"m,l •In fact, there holds

U~,l = l+kK~'m, 1 {AU~_l,O + (1-2A)u~,0 + AU~+l,O} ~

since

,
> X"
Xm,l
m,l
U'
<utI
we getm,n-l= m,n-l 'Urn' ,0 = u".m,O (ini tial condi tion) and
Generally if there holds

1u"m,n -u'm,n = l+kKX"m,n{A(u~_l n-l-u~-l n_l)+(l-2A) (u"m n-l-u~ n-l)""1

- l+kKX')m,n~{AUm',. -.,n- 1+(l.-2A)u' m,n- l+Au' m+1 ,n-. l} -> 0 •
-

6 -

49
§3.

Algorithm{u(R.)}:m,n )\,=0

We shall construct the approximating sequences
and{s~R.)};=O

for FBP.n by the solution of the following

difference scheme:u(R.) _u(R.)
m,n

(3.1)

m,n-l =u(R.)

m-1,n-l -hl<~M-l,

(R.)

= f

l~~N,

O<n<N
= =

(3.2)

u

(3.3)

u M,n = 0

,

{

<Pm

,

l~Lb/hJ

0,[b/hJ+l~M{0l<~[s (.q /hj , O<n<N= ....
- n1[s(R.)/hl+l~.91-l,
O<n<N
n
--

O,n

n

(R.)

(3.4)

u

(R.)

m,O =

(R.)

Xm,n =(3.5)(3.6a)(3.6b)

O<n<N= ={s(R.)}2 = F - 2h2
n
n

M-l

(R.-l)

rn um,n,O~n~N{s(R.)}oo
R.=O
n

and

{u(R.)}oom,n R.=O

l

m=l

,

{s (0) }2 = Fnn

The sequences

satisfy-(3.7)< s(3)<···<s(2R.+1)<···<s(2R.)<···<s(2) < s (0)0 =< s (1) ==
n""
==n
==
==n
....
"""n
n
n
'

(3.B)

(1) < u(3)<···<u(2R.+1)<···<u(2R.)<···<u(2) < u(O)o < u m,n== m,n=== m,n....
.... m,n"".... m,n .... m,n

and

-

(R.)
sn

is monotone increasing in

-

7 -

n

for anyt.

•50This fact can be proved as follows.
First we showGn-

O<S(l)== n •{ (1)}2 = b 2 + 2k
Sn

n-l

I

f

n=O

Put
M-l
h
+ 2h2 [bf ]m. . - 2h 2 I m u (0)
n
m,n
m=l
m
in=l

Then we have
GO = b 2 => 0.

andGn - G
n-l

= 2kf

n-l

-

M-l

I

2h2

m=O

m(u(O)
m,n

u (0)
)
m,n-l

M-l

= 2kf n- 1 - 2k m=O
I m(u(O)-2u(0)+u(O)m-l,n-lm,n-l m+l,m-l

which implies that
Therefore

GnGn>O (O<n<N).

is/monotone increasing inn.

Let us set

By Lemma 1 we haveu(l) < u{O).
m,n = m,nfrom which and Lemma 2 follows

Therefore we haveand{s(2)}2n{s(0)}2 = -2h 2
n

M-l

I

m=l

which implies
(3.9)

8. -

m u(l) ~ 0
m,n

I

51
Using (3.9) a.nd

~eIlUl\a2( we ha.ve

Repeating the similar argument as above, we get (3.7) and (3.8).
If we use (3.7) and (3.8), the monotonicity ofns~) (ti2}

is shown similarly as
If[s (R.) ~ s CR.-I) Inn

becomes sufficiently sma.ll r(S~R.), u~~~'} will be considered an approximate solution of
FBP.m.

... 9 -

in

52
§4.

Numerical Example.
We show the result of numerical solution of FBP.n in the

case

b=l.

We choose the following

functio~sas the boundary

and initial conditions,
t

0<t<2
=== =

,

f (t)
2~t~2.5=T<P(x)=1 - x ,

O<x<l
== ==

In this case, F(t) is calculated by (1.7b) as

F(t) =

Thus we takeK{X

is taken to be

• 7T
8
!+
s~n4 t
7T
3

!+
-7T8
3

0<t<2
== =

,

2<t<2.5== ==equal to 2.5.
220.

,

.

The coefficient of penalty term

The mesh sizes

hand

k

are

= 1/210,k

The results of calculation are shown in Fig. 1 and 2.
fifth iteration of (3.6), we get
(4.1)

max
n

..... 10~After theReferences
1)

J. Douglas & T. M. Gallie: On the Numerical Integration of
a Parabolic Differential Equation Subject to a Moving Boundary
Condition, Duke Math. J., 22 (1955) 557-571.

2)

A. Friedman: Free Boundary Problems for Parabolic Equations I.
Melting of Solids, J. Math. and Mech., 8 (1959) 499-518.

3)

O. Oleinik: On the Stefan-Type Free Boundary Problems for
Parabolic Equations, Seminari 1962-1963 Anal. Alg. Geom.
e ,Topol.,. Vol. 1. 1st. Naz. Acta. Mat., 388-403. Ed. Cremonese.

4)

L. I. Rubinstein:

5)

T. Nogi: On Numerical Solutions of Stefan Problem,(to appear).

6)

J. R. Cannon & C. D. Hill: Existence, Uniqueness, Stability,
and Monotone Dependence in a Stefan Problem for the Heat Equation,
J. Math. Mech., 17 (1967) 1-20.

7)

H. Fujita & N. Sauer: Construction of Weak Solutions of theNavier-Stokes Equations in a Non-Cylindrical Domain, Bull.
Am. Math. Soc., 75 (1969) 465.

The Stefan Problem, Am. Math. Soc., (1971).

Hideo KAWARADA
Department of Applied Physics,
Faculty of Engineering,
University of Tokyo
Makoto Natori
Computer Centre,
University of Tokyo- 11 -

t
2.5
uex,t)
2.0

1.0

1.5

1.0

0.5

0.5

t=2.0
0.0~c-'-_---.L----------

1.0

1.5

~(t)

Fig. 1

2.0

0.0 "----~--.II---~---""--"-O
0.5
1.0 . X 1.5"
2.

Fi &,.2