NUMERICAL SOLUTION OF LINEAR SECOND ORDER FAlRABOLIC EQUATIONS BY THE METHOD OF COLLOCATION WITH CUBIC SPLINES. by 1USEBIUS JACOBUS. DOEDEL B. Sc., Un i v e r s i t y of B r i t i s h Columbia, 1971 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOE THE DEGREE OF MASTER OF SCIENCE i n the Department of Ma&hematics We accept t h i s thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA A p r i l , 1973 In presenting t h i s t h e s i s i n p a r t i a l f u l f i l m e n t of the requirements f o r an advanced degree at the U n i v e r s i t y of B r i t i s h Columbia, I agree that the L i b r a r y s h a l l make i t f r e e l y a v a i l a b l e f o r reference and study. I f u r t h e r agree that permission f o r extensive copying of t h i s t h e s i s f o r s c h o l a r l y purposes may be granted by the Head of my Department or by h i s r e p r e s e n t a t i v e s . I t i s understood that copying or p u b l i c a t i o n of t h i s t h e s i s f o r f i n a n c i a l gain s h a l l not be allowed without my w r i t t e n permission. Department of Mathematics The U n i v e r s i t y of B r i t i s h Columbia Vancouver 8, Canada Date A p r i l , 1973. i i Abstract C o l l o c a t i o n with cubic splines i s used as a method for solving Linear second order parabolic p a r t i a l d i f f e r e n t i a l equations. The c o l l o c a t i o n method i s shown to be equivalent to a f i n i t e difference method that i s consistent with the d i f f e r e n t i a l equation and stable i n the sense of Von Neumann. Results of numerical computations are given, as w e l l as an a p p l i c a t i o n of the method to a moving boundary problem f or the heat equation. i i i Table of Contents INTRODUCTION 1 PART I. 1.1 Description of the C o l l o c a t i o n Scheme .... 4 1.2 The F i n i t e Difference Scheme 9 1.3 Consistency of the Difference Scheme ..... 19 1.4 S t a b i l i t y of the Difference Scheme 21 1.5 Equivalence of the F i n i t e Difference Scheme and the C o l l o c a t i o n Scheme... 24 1.6 The Algorithm f o r Solving the C o l l o c a t i o n Equation 35 PART I I . Results of Numerical Computations 43 PART I I I . A p p l i c a t i o n to the Stefan problem 55 111.1 Introduction 55 111.2 Description of the Scheme 57 111.3 Numerical Results 68 Bibliography 74 i v Acknowledgement The author s i n c e r e l y wishes to thank Professor James M. Varah f o r generous assistance during the course of the research. He wishes to thank a l s o the National Research Council of Canada f o r bursary support and the U n i v e r s i t y of B r i t i s h Columbia for teaching assistanceships. He i s also thankful to Mrs. Magdalena Alemparte f o r bhe typing of his manuscript. -1-INTHODUCTION. Co l l o c a t i o n with piecewise polynomial functions has been developed as a method for solving two-point boundary value problems f o r ordinary d i f f e r e n t i a l equations by f o r example B a s s e l l and S'hampine {* 13J , and de Boor and Swartz (3j • Recently Douglas and Dupont£5Jhave proposed a c o l l o c a t i o n method with Hermite piecewise - cubic poly-nomials f o r solving non l i n e a r parabolic equations. In t h i s thesis an e f f i c i e n t method f o r solving l i n e a r second order parabolic equations by the method of c o l l o c a t i o n with cubic splines w i l l be presented. In part I the c o l l o c a t i o n procedure w i l l be introduced and shown to be equivalent to a f i n i t e difference scheme that is second order accurate. With a uniform mesh and constant c o e f f i c i e n t s the f i n i t e difference scheme w i l l be shown to be stable i n the sense of Von Neumann, so that f o r t h i s type of problem the sol u t i o n 6f the f i n i t e d i f -ference equations converges to the sol u t i o n of the con-tinuous problem with second order accuracy. -2-The scheme has the advantage that l i k e K e l l e r ' s Box scheme i t does not require a uniform net. Ap-proximate values of f i r s t and second derivatives at the net points can be obtained e x p l i c i t l y from the so l u t i o n of the f i n i t e difference equations. Convergence of these to the corresponding derivatives of the continuous so l u t i o n w i l l not be shown; however, numerical computations seem to indi c a t e that they are second order approximations, which is i n agreement with the results obtained by K u s s e l l and Shampine. F i n a l l y the matrix equations r e s u l t i n g from the f i n i t e d ifference equations involve only t r i d i a g o n a l matrices that are of the same dimension as those of the corresponding Crank Nicholson scheme. In Part II a number of problems w i l l be solved using the scheme and the r e s u l t s w i l l be compared to those obtained by the Box scheme and the Crank Nicholson scheme. In Part III we s h a l l i n d i c a t e how the scheme can be used to obtain an approximate solution to a Stefan problem and some numerical r e s u l t s w i l l be given. It i s not neccessary to r e s t r i c t to cubic splines i n the c o l l o c a t i o n scheme. Good results have been obtained using Hermite piece-wise q u i n t i c polynomials for example. -3-The cubic spline c o l l o c a t i o n scheme wins, however, as f a r as s i m p l i c i t y i s concerned. -4-PART I. 1.1 Description of the C o l l o c a t i o n Scheme. Consider the l i n e a r second order parabolic P. D. E. ( p y ) + = y + qy + r y + f ( 0 $ x $ l ) ^ XX X with p, q, r and f functions of x and t . Upon setting z = e .y we f i n d that s a t i s f i e s + e . f We assume that the i n d e f i n i t e i n t e g r a l s J Aqdx and J q^dx can be evaluated, so we may r e s t r i c t our att e n t i o n to equations of the form (1.1) (1.2) ( p y ) t = y x x + qy • f subject to y(x.o) = g ( x ) y(o.t) = g (t) o y ( l , t ) = 6 l ( t ) ( 0 N < x * l ) (0 n =>*(x..t ) j J n n (>*! ) =. >^ (x , t ) where the subscript z z j j n represents a f i r s t or higher order d e r i v a t i v e with respect to x or t . J J n ;. For functions defined on £ 0,1 J x £ t _7 • ' y n = />(x,t ) n ( y 3 , ) n , y>. (x ,t ) z j z j n ; ) n = (x,t ) z z- n Given a p a r t i t i o n 0 = x < x , < ^ - x . - l 0 1 J of the i n t e r v a l £ 0,1 ~J , a spline of degree m defined on -7-Qo,l3 i s a function S 6 C \_ 0,\~] that i s given by some polynomial of degree m or less i n each subinterval [x.,x. . S\ 3 3 (0 4- 3 4 J - ^ l ) . The points fx 1 J are c a l l e d the 1 J J j = 0 knots or j o i n t s of the spline f u n c t i o n . In p a r t i c u l a r i f S i s a cubic spline defined on the unit i n t e r v a l then 2 S e C f O - l l and i n each subinterval Tx , x 1(0 i i S j - l ) L j j + r S i s given by a polynomial of degree 3 or l e s s . For an introduction to spline functions see for example [*9l • A more rigorous treatment of spl'.nes can be found i n f 1~\ . The c o l l o c a t i o n procedure i s now defined as follows. At each time l e v e l t ( l £ n ^ N ) we look f o r an approximate n n solution S to the problem (1.1, 1.2) s a t i s f y i n g (1.3) i s a cubic spline defined on [ o , l l x \. n^~\ (1.4) The boundary conditions: 8 J • 6 i ( V (1.5) The c o l l o c a t i o n conditions: f s n - f P " 1 S*" 1 1 f J j * J j r k n n n n n n-1 n-1 n-1 n_i) (S ) . + q.s. + f . + (S )• + q 4 S. + f . (0< H J) -8-0 Vftiere we l e t S be the cubic sp l i n e interpolant of the i n i t i a l data g(x), ( i . e . s g(x.) ( 0 ^ j ^ J ) ) s a t i s f y i n g the endpoint conditions: 0 and (S ) - g (1) XX J XX In the special case of the heat equation t XX the c o l l o c a t i o n equation (1.5) becomes 1 n 2 «. (s ) n + (s f - 1 v xx'j v xx'j which resembles the Crank-Nicolson difference approximation n n-1 IT - U. 1 .i J 211* • U n ' u n - ^ u , - if\'\ - 2irn"1 • u1^; J - l J J * l J " 1 J J * i - l ^ T T n - l k 2 n Except that only the time v a r i a b l e has been d i s c r e t i z e d . Usually the term c o l l o c a t i o n i s meant to imply n that the approximate s o l u t i o n S s a t i s f i e s the d i f f e r e n -t i a l equation exactly at a c e r t a i n number of points. Because of the d i s c r e t i z a t i o n in time t h i s i s a c t u a l l y not true f o r condition (1.5). However, we s h a l l r e f e r to condition (1.5) a s the c o l l o c a t i o n condition. -9-1.2 The F i n i t e Difference Scheme. Since at each time l e v e l t the approximate so l u t i o n n S should be a cubic s p l i n e , i t i s determined uniquely by four c o e f f i c i e n t s per subinterval. That i s we could represent S n by S n(x) = q . ^ • b . > n x * c . ^ x * d . ^ ( x . . ^ x ^ ) ( U J 4 J ) This would lead to 4«J unknown c o e f f i c i e n t s per time l e v e l . These c o e f f i c i e n t s then could be obtained from the 2 boundary conditi->ns ( 1 . 4 ) , the J + 1 c o l l o c a t i o n conditions (1.5) a n-d 3(J-1) continuity conditions. Such a procedure would give large systems to be solved, and i s therefore not advisable. Another approach would be to express S n i n terms of basis-functions having the property that a l i n e a r combination of these automatically s a t i s f i e s the con-t i n u i t y conditions. Such basis-functions have been constructed, see f o r example ^14^ • Vfe s h a l l not use either of the above methods; instead the c o l l o c a t i o n procedure w i l l be shown to be equivalent to a f i n i t e d ifference method. -10-F i r s t t h i s f i n i t e difference scheme w i l l be stated and s u f f i c i e n t conditions w i l l be given f o r the existence of a unique sol u t i o n to the r e s u l t i n g matrix equation. In section 1.3 we s h a l l show that the f i n i t e difference scheme i s consistent with the continuous problem. (1.1), and i n section 1.4 the scheme w i l l be shown to be stable i n the case of constant c o e f f i c i e n t s and a uniform net. The equivalence of the f i n i t e difference scheme and the c o l l o c a t i o n scheme w i l l be s e t t l e d by theorem 5»3 cf section 1.5• The difference scheme w i l l then serve as the basis f o r the algorithm f o r solving the c o l -l o c a t i o n scheme. The algorithm w i l l be given i n section 1.6. Consider the difference scheme defined by J - l , 1$ n$N) (2.1) U n = g (t ) 0 0 n (1$ n « N) (1* n$N) U (0« M J) -11-where ,h. * h. \ k / \ k V 3 j+1 n \ n h 4h. . n rrv-"^ . . " 5 - . . ^ i - ^ 1 : : ^ - ^ n ~ h* TTT J > 1 1 (h * ( h j • h j + l ) :•. h h +h ( n n n-1 n-1 1 q, u - q. ,u l h . J 3 J + 1 , 3-1 j - i y v j j , h .+h , j j+1 n n-f +• f . J+l J H -12-This difference scheme leads to matrix equations (n) (n) (n) where (n) u — D (Kn N , n n u = (U Q , u , n T and (») n n — = 0 (14. j 4 J - l ; U n« N) Now 2(h.+h. ,) 1 1 h.-h.., n l s J J+l „ , J J+l H > j a +- — 4- n and h J ^ - l 3kn J * 1 ^ ^ n* h j 3k n' + - ( S q n (Since 3 ^n * j - l 1 h / h • S i m i l a r l y n h j + l j + r * j+l n * •+• — q 3kn 6 -16-n Hence b n ( h . + h j + l ) q n h 3*n h J + l * „n k. n n > 0 since by assumption <( in Cor. 2.2 In the s p e c i a l case when p = 1, q = 0 i n (1.1), i . e . f o r the inhomogeneous heat equation y t = yxx - f ( * . * ) -17-subject to the i n i t i a l and boundary condition (1 . 2 ) there i s always a unique sol u t i o n to the matrix equations ( 2 . 2 ) , Proof. With p =. 1, q = 0 the c o e f f i c i e n t s ( 2 . 3 ) become h . 1 ( 2 . 5 ) a" 3k h j 2 ( h j + h j + 1 ) 1 h. . 1 d j = 8l(*n) -18-hi 1 \ 3kn h j / ( 2 ( h O + h J + l } 1 1 N 3k. - l n h - i 1 \ n 3^ ^ ( h j t h . ^ ) 1 (1< J « J - l ) - n c , 3 1 1 h. 1 J V 3kn h j 3kn -4-1 h j 1 h 1 3k n h j h j * i 3kn h. 3kn > 1 h J + h J - l 3k n -19-1.3 Consistency of the Difference Scheme. Theorem 3.1 The difference scheme (2.1) i s consistent to second order accuracy with the problem (1.1, 1.2). Proof. Define the l o c a l truncation error 7T n by 3 17° = V -where y i s the solution of the continuous problem (1.1, 1.2). Taylor expand to derive the f o l l o w i n g : ( 3 . i ) ffi - j»-Vfi .( ( ? y ) t ) ^ . un),)]'1 < « » W j " * k n 2 12 + higher order terms. -20-/ o ~. \ n n (,.2) y ^ - y . J + 1 n n«*. y i " y • J J - l |(h .fh * v J J * l y "( yxx) J - l + 2 ( y x x)3 +-J+l "(yxx)n J J+l h .+h . J J+l n j+l 3 •} J J+l 12 \ h .+h. . J J>1 ( y x X X X ) ^ f higher order terms Using the i d e n t i t i e s (3.1) and (3.2) and the f a c t that y s a t i s f i e s the d i f f e r e n t i a l equation ( 1 . 1 ) we see t h a t : <3'3> T n 1 j h .+h 3 3 h f h J j j+l j J+l ^xxxx-'j * " x x x x M k n 36 C(nv),..) ^ 2 _ _ ^ y W j - i . f C ( p y ) t t t ) ° 1 ; J j+l h . , i J C(py) 1 n ~ 2 ^ j _ J < p y > t t t ) j f l j f h j + i f higher order terms -21-If we l e t h ; max h and 0 £ ;ji2.J k = max k n l ^ n ^ N then (3*3) shows that 7J n. = 0 ( h 2 + k 2 ) 1.4 S t a b i l i t y of the Difference Scheme. In t h i s section i t w i l l be shown that the d i f -ference scheme defined by (2.1) applied to the problem (4.1) v. i s stable i n the sense of Von Neumann. We have to make the following assumptions: (4.2) The net H h(T) i s uniform. ( i . e . h = h ( l ^ j ' ^ J ) ; k n = k ( l £ n * N ) pyt = y x x + q y y(x.o) = g(x) y(o,t) = M*) y(i.t) = S i ( t ) -22-( 4 « 3 ) For each n (O^-n^N) the approximate sol u t i o n ( O ^ j ^ J ) is periodic with period 1. (A»A) The c o e f f i c i e n t functions p and q are constant with p^O and q^O. The basic idea of a Von Neumann sta b i l i t y -a n a lysis i s to write the f i n i t e difference s o l u t i o n as a Fourier series and then to show that none of the components can grow i n amplitude as time increases. Under the assumptions (4.2, 4«4) the d i f -ference equation f o r (4«l) i s ' ( 4 . 5 ) ( hp _1 _ hq 3k h 6 1 U j - l -r 4hp | 2 2hq 3k h 3 3 k h 6 U>1 ^ ^ I ^ U n - 1 ^ ^ _ £ + B ^ ] n - l + hp + l ^ h q 3 J v 3 k h 6 3k h 6 J \. 3k h By ( 4 * 3 ) we can write ' n. . y n iw2fT3i ( 4 . 6 ) ^ ( x ) ^ °. W 2 - OQ n-1 iw27Tx e w where { c 1 1 V and { c 1 1 ^ \ are the F U n(x) and ^ " ^ ( x ) r e s p e c t i v e l y . ourier c o e f f i c i e n t s of n Uj+1 - 2 3 -Using ( 4 . 6 ) i n ( 4 . 5 ) w e deduce that the Fourier coef-f i c i e n t s must s a t i s f y the following r e l a t i o n : ^ n A ? 4. 1 ^ h < l \ g - ^ ( ^ ^ h p ^ i , 2 h q \ /hp 1 ,n - l ~ , _ (IgE. - *g ) ^ ( 1 + ft) s i n * . h q ) ( ^ E _ ^ s i „ ^ / a ) ft) s i „ ^ / 2 _ h q) Upon studying t h i s l a s t expression and r e c a l l i n g that p ? 0 and q< 0 vs conclude that n c w n-1 c w £ 1 e ( 4 . 7 ) ^ f hp 1 hq\ - i l / 4hp , 2 2 h q \ / h p 11 hq ] i ^ Ik ~ h~ ~ £*"/ f 5 \ ~ 3 T h ~ ~ J " y 3k" ~ h" 6" where (• x 2 > wh Using the f a c t that e 3 + e 3 - 2 ( l - 2 s i n J / 2 ) we get e and hence s t a b i l i t y has been established. - 2 4 -1.5 Equivalence of the F i n i t e Difference Scheme and the C o l l o c a t i o n Scheme. The purpose of t h i s section i s to show that the c o l l o c a t i o n scheme defined by ( 1 . 3 , 1.4» 1 « 5 ) i s equivalent to the f i n i t e difference scheme ( 2 . 1 ) . Before e s t a b l i s h i n g t h i s f a c t i n Theorem 5»3 we s h a l l state a number of w e l l known properties of cubic splines that are needed i n the proof of Theorem 5«3« Since we have already shown that the difference scheme has a unique solution provided the conditions (2.4) hold the equivalence of the two schemes w i l l then imply that under the same conditions there i s a unique solution to the c o l l o c a t i o n equations. This consequence i s given i n C o r o l l a r y 5»4« The f i n i t e difference scheme only gives ap-proximate values of the s o l u t i o n at the mesh points. However, i n C o r o l l a r y 5 - 5 we s h a l l show how approximations to f i r s t and second derivatives at the mesh points as w e l l as approximations to the s o l u t i o n and it<$s f i r s t two derivatives at intermediate points can be obtained e a s i l y from the s o l u t i o n to the difference equations. - 2 5 -Lemma 5«1 Let ( 0 - ,1~ J; l^n<-N) be the unique solu t i o n of the f i n i t e difference equations (2.1) subject to the conditions (2.4). Then f o r each n ( l ^ n ^ N ) there i s a unique cubic spline S that interpolates the of the mesh points k. ( 0 - j — j ) and s a t i s f i e s 3 ( s ,n = 2 V P 0 S 0 ' P 0 S 0 / £ . n . x n-1 n- X ( Sxx)3 - 2 n n n-1 n-1 "\ J , n n n . . n-1 n - l _ n - l in-1 where ( 8 ^ ) ° = g ^ O ) . ( S j J : ^ ( 1 ) , .0 Q S u z g(0) and S i r g ( l ) . 0 u Proof. C l e a r l y the assertion follows i f we can show that for a r b i t r a r y data-points with J J J j - o 0 1 x < x, <.....<. x ^ 1 0 T - 2 6 -and for a r b i t r a r y numbers °^ a n d t h e r e exists a unique cubic spline S such that S(x ) = y ( 0 S j <: J) S 3(x ) - ° L~ XX J / But t h i s i s a standard i n t e r p o l a t i o n r e s u l t , that can e a s i l y be shown by making use of r e l a t i o n ( 5 « 1 ) of the following Lemma. Lemma 5«2 Let S be a cubic spline i n the i n t e r v a l ^ 0 ,1 J •vrith knots 0 = x Q < x^< x - 1 Then S s a t i s f i e s 1 fh- ( h - + i (5.1) " I - kc'j-l + ^ x x ' j ( S * * > j H 3 [ W W i j+i j J j - i h . , h . J+l J -27-( 5 . 2 ) ( 5 - 3 ) S. - S, j " — + - { 2< Sxx).i + (Sxx).i-l h. 6 • ( 1 * j £ J) ( S x ) j - 2 ( S x x > j +" ^ x x ^ l f ( 0 £ j ^ J - l ) Proof. These i d e n t i t i e s are easy to obtain, e.g. by using Taylor expansions f o r S noting that S - 0 f o r k ^ 4 dx See f o r example 9 pp. 31-32~J . Theorem 5 . 3 Let [ u n j ( 0 ^ 3 6 J 1 •=. n =^ N) be the unique so l u t i o n of the difference equations ( 2 . 1 ) subject to the conditions ( 2 . 4 ) . -28-(5«4) For each n (1- n^N) l e t S U be the unique cubic spline interpolant of j^U* J ( 0 * J * J . 1— n^ =. N) s a t i s f y i n g a I — / ^ - $. ( s ^ r - w - c n n ri-1 n-1 p S - p T S XX J I L n n n ^n . x n - l n-1 n-1 £n-l k - J 0 (5»5) Let S be the unique cubic spline interpolant of the i n i t i a l data g(x) at the points ^ x^ . ^ ( O ^ j ^ r j ) s a t i s f y i n g (S )° = g (0) xx 0 x x (s x x ) ° - g (1) x x J XX Then f o r each n ( 1 — n^-N), S n s a t i s f i e s the c o l l o c a t i o n equations (1.3» 1-4, 1.5). Conversely assume that S U i s a s o l u t i o n of the c o l l o c a t i o n equations (1.3» 1.4» 1.5). Then f o r each n ( 1 — n £ : N), S s a t i s f i e s the difference equations (2.1). - 2 9 -Proof. Let-^SnJ (0£.n£:N) be the cubic splines as defined i n 5 » 4 and 5 » 5 . (These are uniquely determined by Lemma 5»1)« Then S n s a t i s f i e s the difference equations (2.1) f o r (l^=n — N ) . Also Sn s a t i s f i e s the r e l a t i o n (5.1) f o r ( O ^ n ^ N ) . Hence fo r ( l ^ n ^ - K ) we have that 1 f h. 'P1} ,S n -p^V"1 J - l J - l j-1 j-1 3 h.+h. . 1 J J+l / J J J j j+i , n „n n-1 n-1 p 5 -p S j+l j + l j+l j41 W l n 1 h ( 5 . 6 ) 1 J *Sxx)il * 2 h j j+l - ( S x x ) j + l J , n _n n-1 n - l x , 0 / _n n , n-1 n - l N ( V i V i ^ d - i ^ +•2 ( q j s / q j s J ) h.+h I J j+l i n n n-1 n-1. c,n ,«n-l, j+l r.n n-h .th . , h .+h 0 j+l - 3 0 -Let n n n n-1 n-1 ^ s . P s -p s i r r B . j j j j , 'ri n n n Nn-1 ( n - 1 n-1 n-1 L S J { ( S ^ J j + q ^ j f f J K S J J +qd S j t f j k n 2 [ Then ( 5 « 6 ) becomes (5.7) J — t s ^ ^ T ^ + ^ i - 1 ^ 1 = 0 h .+h , h 4-h J J»l 0 j+l ( 1 £ i j ^ j - l , \6z N ) . and by hypothesis 5«4 we have that 2 V - 2! S n r 0 (1 n ^ N ) . Yfe can consider (5»7) as a system of l i n e a r equations i n the unknown j L s n { ^ for each ( 1 — n N) with (. j J J •= 0 the right-hand side being the zero vector. Since h h . , we must have that j J + l 4- : K 2 h.fh h +h J j+l J j+l L S n - 0 (0 £ j ^ J , l ^ n ^ H ) j - 3 1 -But t h i s implies that f o r each n ( 1 < n N) S s a t i s f i e s the c o l l o c a t i o n equation ( 1 . 5 ) . C l e a r l y f o r each n ( 1 ^ n ^ N) S n a l s o s a t i s f i e s the boundary conditions (1.4) and hence the f i r s t a ssertion has been proved. Conversely i f we assume that f o r each n ( 1 4 n ^ N) S n i s a solution of the c o l l o c a t i o n equations ( 1 . 3 . 1 . 4 . 1 . 5 ) then This r e l a t i o n , however, i s i d e n t i c a l to ( 5 « 6 ) and using ( 5 . 1 ) i t follows that each S n s a t i s f i e s the difference e-quation ( 2 . 1 ) . As an immediate consequence of Theorem 5 . 3 w e have the following: Ls = 0 (o 4 j 4 J. j U n U ) Therefore - 3 2 -Cor. 5«4 Under the conditions ( 2 . 4 ) of Theorem 2.1 there is a unique solution Sn (l 1«4» 1 « 5 ) with Sn =r Un ( 0 4 j < J; 1 4 n < N) j j where if; ( 0 ^ j ^ J; U n < N) is the unique solution of the difference equations (2.1). Cor. 5«5 Given the solution U1! ( 0 < j i J; 1* n 4N) of the finite difference equations (2.1) and setting Sn = U11 (0 4 j < Jj U n U ) we can explicitly compute (S )n and (S )n (0 < j U ; 1 S n < N) by x x j x J -33-( 5 . 8 ) (8 f xx j = 2 n n n > - Q.S ^ (8 ^ q-^-fS" 1 ] - xx J J J J (0 .< j 4 J; U < n 4 N) ( 5 - 9 ) n n S -S h J+l j J+l j+l ( 0 <. j 4 J - l ; K< n ( 5 . 1 0 ) ) n s n - s n • J J - l ~~J = 2(S ^ (S ) n T . v xx'J + v xx'J-1 where S i s sgain taken to be the cubic spline interpolant of the i n i t i a l data g(x) at the points I x . j ^ s a t i s f y i n g : L J J j = 0 and (s r = g (o) xx 0 xx (S )° = g ( 1 ) XX J XX Moreover i f x < x < x then j j - l -34-X j + l ~ x x ~ x j ( 5 . 1 1 ) (S Ax) = - ^ - r - (S ) n - + — ( 8 f xx h XX j h xx j + l j+l j+l (5.12) „ n „ ( ( ^ ) J U - < S » 2 h. J f l ( 5 . 1 3 ) n (x -x)S n+(x-x )S S , ^ J J J J + i h j + l (x r^ (x-x ) (S Ax) . (S + (S )* x xx' *>'v 4- xx;j 1 v xx'j+l Proof. ( 5 . 8 ) follow from the c o l l o c a t i o n condition ( 1 . 5 ) , ( 5 . 9 ) and ( 5 . 1 0 ) follow from Lemma ( 5 . 2 ) , and ( 5 . 1 1 ) , ( 5 . 1 2 ) and ( 5 . 1 3 ) are simple i n t e r p o l a t i o n formulas. For example ( 5 . 1 1 ) represents l i n e a r i n t e r p o l a t i o n , since the second d e r i v a t i v e of a cubic spline is a continuous piece w-ise l i n e a r function. -35-1 . 6 The Algorithm for Solving the C o l l o c a t i o n Equations. The algorithm for solving the problem (1.1) subject to the i n i t i a l amd boundary conditions (1.2) i s i n f a c t implied by Cor. 5«4 and Cor. 5«5» Below i t w i l l be siteafced more p r e c i s e l y , using the standard L U decomposition for t r i d i a g o n a l matrices (2.2) Given with x J 1 t r T N and netspacings h . i x . -: ( 1 ^ j ^ J ) k t - t n-1 (1 ^ n 6 M) n n the algorithm w i l l give approximations -36-n U = v, = (s ). n xx to y to (y ) n W. E. (S ) to (y ) n xx (0 — j — J; l ^ n ^ F i r s t we have to compute the cubic spline 0 , . interpolant S of the i n i t i a l data g(x) at the points f x "I s a t i s f y i n g the endpoint conditions } j J j = 0 (s )° - g (10) xx 0 x x and (S ) - g Y (1) XX ,' j — xx However, since we only need (S ) (0 j ^z. J) for xx j subsequent computations, these values can be computed e a s i l y by making use of r e l a t i o n (5«1) of Lemma 5*2. This r e l a t i o n leads to a matrix equation of the form: -37-(6.1) (2 ^ j ^ J-1) (1 & j ^ J-1) (1 £: j-=: J-2) and with 0 c . 0 h. h +h. , 0 b - 2 h ... 3 + l h.+h. , J j+l -38-d - 6 2 ) - g ( x 1 ) g ( x 1 ) - g ( x Q ) h 2 h l / i \ • h 2) Sxx(x(P 0 c. g(x.+1)-g(x.) gCx^-gCx.^) (h . f h. ,) (2 £ j 6 J-2) ^ J ^ W . ^ S ( X J - 1 ) _ 6 ( X J - 2 ) n £ j-i 6 x x v J ' (h * h ) V J-1 J1 J h tti J-1 J Applying the L U decomposition algorithm f o r t r i d i a g o n a l . 0 .rQ matrices, with intermediate quantities oC , ft and 0 j j w , we can now solve the system (6.1) by: j -39-w w 1 0 = (d ( 2 ^ j<£ J - l ) (S ) = E (x ) - g (1) XX j XX J XX (S ) - w xx J - l j _ i 0 0 , 0 0 (S ) - w - V- (S ) x x j j 6 j xx j+l (J-2 •1) -40-(S ) - S (* ) (0) xx 0 xx o xx 0 o We now i n i t i a l i z e U and W by j 0 Let U -n-1 (S )° XX j (o — j ~ J) n n n n Next compute the c o e f f i c i e n t s a ., b , c and d of 3 J j j the matrix equation ( 2 . 2 ) as defined by ( 2 . 3 ) , or by ( 2 . 5 ) i n case of the inhomogeneous neat equation. Again to solve ( 2 . 2 ) we use the L U decomposition n _n algorithm, with intermediate quantities , / and n j j w , which leads to the following computations: 0 -41-n n n n / n )r1-\lJ~ i .3 .i J * ' n 4/.n n n / . n j r . z c . . ( 2 ^ j-S.J-1) ( 2 - j - b J-2) w11 = dn -a n g (t ) 1 1 1 Ov n' wn - ( d n - a V ) I JL " (2£: j ^ r J-2) j J j j-1 ' J n n n n n . / , w - (d - c g (t ) - a w ) I oL. J - l J-1 J-1 1 n J-1 J-2 / J-1 The approximations to the solution and i t s f i r s t two derivatives are now obtained as f o l l o w s : U - w J-1 J-1 n n i_n n U ^ w - / U (J-2 ^ j ^ l ) •1 j j j+l u n ^ g (t ) 0 o n -42-p^-p^V" 1 2 JJ J J i j j - J - J ~ J J IT •n n j+l j J+l j j+i j+l (o ^ j — J) (0 ^ j ^ J-1 n J n 2 W J + S t a r t i n g at (*) we repeat the same computations with n increasing from 2 to K, -43-PART I I . Results of Numerical Computations. The c o l l o c a t i o n scheme was used to obtain ap-proximate solutions to four d i f f e r e n t i a l equations. For each problem solutions were computed with various mesh si z e s . The r e s u l t s are given below, as well as those Crank-Nicolson scheme [ 6 j to these problems. In these tables J denotes the number of mesh-i n t e r v a l s , both i n the x - d i r e c t i o n and i n the t - d i r e c t i o n . Each problem was solved on the unit square, so that f o r a uniform mesh the mesh-spacings i n the x- and t - directions are given by l / j . A non-uniform mesh was used i n Problem 4 only. Further by e, e and e we denote the maximum x xx error i n the approximate s o l u t i o n and i n the approximate f i r s t and second derivatives r e s p e c t i v e l y , the maximum being taken over the mesh points. The quantity between brackets i s the observed trate of convergence, obtained by computing the value of obtained by applying log log -44-Problem 1. The f i r s t problem was to solve the following homogeneous heat equation: y - y (0 6 x ^ 1 ; O i f i l ) z xx y(x,0) z e X t y(0,t) = e / ,\ t+.l y ( l , t ) = e x+t which has y(x,t) = e as i t s unique s o l u t i o n . The mesh was taken to be uniform, i . e . h = k - l / j . The r e s u l t s are given i n Table 1. We note that solution of the c o l l o c a t i o n scheme converges at the rate of 4.0 to the s o l u t i o n y(x,t) of the continuous problem at the mesh points. Such an accuracy cannot be expected i n general, however, The reason that i t happens here is that ap-parently f o r th i s p a r t i c u l a r problem the next two terms of the error expansion (3.4) cancel. Further we note that the error e for the Box-scheme i s approximately one h a l f of the error e f o r the Crank-Nicolson scheme. -45-Froblem 2. The second problem was to solve the following inhomogeneous heat equation: y - y 4 cos(xft) - s i n ( x t t ) (0 ^ x 1; O ^ t t _ xx y(x,0) - sin(x) y(0,t) - s i n ( t ) y ( l , t ) sin(t+l) which has y(x,t) = sin(x4-t) as i t s unique s o l u t i o n . As i n Problem 1, t h i s s o l u t i o n i s p a r t i c u l a r l y smooth and well-behaved so that good r e s u l t s may be expected. The mesh was taken to be uniform. The r e s u l t s are given i n Table 2, and indi c a t e that the cubic spline c o l l o c a t i o n scheme i s second order accurate i n i t s approx-imation to the solution as w e l l as i n i t s approximation to the f i r s t and second de r i v a t i v e of the s o l u t i o n . In t h i s problem the Crank-Nicolson scheme gives the best approximation to the s o l u t i o n . - 4 6 -TABLE 1 • Cubic Spline Scheme Box Scheme C. N. J e e x e XX e e X e 5 . 4 5 6 " 5 • 5 6 4 " 3 . 2 5 5 " 1 - 2 .173 . 1 6 2 " 1 . 3 4 0 " 2 10 ,287~ 6 . 6 7 7 " ^ . 7 8 3 ' 2 . 4 3 2 J - 2 . 4 0 4 . 8 6 l " 3 ( .4 .00) ( 3 . 0 6 ) (1.70) (2.00) (2.00) ( 1 . 9 8 ) 2 0 .181" 7 . 8 3 3 _ 5 - 2 . 2 0 3 - 3 . 1 0 9 * - 2 . 1 0 1 - 3 .217 J ( 4 . 0 0 ) ( 3 . 0 2 ) ( 1 . 9 5 ) (2.00) ( 2 . 0 0 ) ( 2 . 0 0 ) 4 0 - 8 . 113 . 1 0 3 " 5 - 3 .517 J - 4 .272 * • 2 5 3 " 3 - 4 • 543 ( 4 . 0 0 ) ( 3 . 0 2 ) ( 1 . 9 7 ) ( 2 . 0 0 ) ( 2 . 0 0 ) (2.00) -f -47-TABLB 2 i Cubic | 'j Spline Scheme Box Scheme C. N, J J e e X e XX e e X e 5 • 5 6 l " 3 „ -2 .257 . 6 6 3 * 2 - 3 .970 J -2 .777 - 3 . 3 3 0 J 10 . 1 4 2 - 3 (1.98) . 6 1 4 " 3 (2.07) . 1 6 6 " 2 ( 2 . 0 0 ) .213' 3 ( 2 . 1 9 ) .177~2 ( 2 . 1 3 ) .826" 4 (2.00) 20 •355~ A (2.00) .. - 3 .148 J (2 . 0 5 ) .417" 3 (2.00) . 5 2 9 " 4 (2.01) .428' 3 ( 2 . 0 5 ) .208"^ (2.00) - 4 8 -Froblem 3. The t h i r d problem was to solve 400(x-t+l) _ 2(400(x-t) ) 2 ( l f 2 0 0 ( x - t ) S l ) 2 ( l + 2 0 0 ( x - t ) 2 ) 3 ( 0 £ x £ l j 0 ^ - t ^ 1) y(*,o) - l y(o,t) z y(-i,t) = 1 - - . 3 2 lt200x 1 l+200t 2 1 l + 2 0 0 ( l - t ) 2 which has 1 y ^ x * ' 2 a s i t s unique solution. l+-200(x-t) -49-The s o l u t i o n of t h i s problem involves a narrow wave t r a v e l l i n g along the l i n e t = x. The mesh was again taken to be uniform, since what may be a good choice of mesh points i n the x - d i r e c t i o n at one t i m e - l e v e l , w i l l be a poor choice at the next t i m e - l e v e l . The numerical r e -s u l t s are given i n Table 3« For t h i s problem the cubic spline c o l l o c a t i o n scheme can be seen to give the best approximation to the solution and i t s f i r s t d e r i v a t i v e . The Box scheme shows a large error i n the ap-proximation to the f i r s t d e r i v a t i v e , when U i s small. In f a c t heavy o s c i l l a t i o n was observed i n the approximation to the f i r s t d e r i v a t i v e of the s o l u t i o n , even where t h i s d e r i v a t i v e should be small. It has been noted i n f_7 pp.25-26 that the Box scheme does lead to such o s c i l l a t i o n s i n c e r t a i n problems with boundary conditions as above. The rate of convergence does not approach 2 . 0 as c l e a r l y as i n the preceding problem, however, t h i s i s u s u a l l y the case when less well-behaved functions are used to t e s t a scheme. - 5 0 -TABLE 3 Cubic Spline Scheme Box Scheme C. N. J e e X e XX e e x e 10 . 4 6 4 + 1 . 2 4 9 + 2 . 6 3 2 + 2 • 5 3 3 + 1 6 6 4 8 + 2 . 5 2 4 + 1 20 . 2 2 0 ° ° [ 4 . 4 0 ) . 1 3 8 + 1 (4.17) .147+2 ( 2 . 1 0 ) . 3 5 2 ° ° ( 3 . 9 2 ) . 2 3 0 + 2 ( 1 . 4 9 ) 00 . 3 4 4 ( 3 . 9 3 ) 40 . 2 6 4 " 1 [ 3 . 0 6 ) . 9 0 3 ° ° ( 0 . 6 1 ) . 2 3 2 + 2 ( - 0 . 6 6 ) . 4 3 9 " 1 (3.oo) . 3 5 0 + 1 (2.72) • 3 2 5 " 1 ( 3 . 4 0 ) 80 .617"2 [ 2 . 1 0 ) . 2 2 0 ° ° ( 2 . 0 4 ) . 1 0 9 f 2 ( 1 . 0 9 ) -2 . 9 2 2 ( 2 . 2 5 ) .502°° (2.80) . 6 5 4 ~ 2 ( 2 . 3 1 ) - 5 1 -Problem 4. F i n a l l y we consider t h ^ following d i f f e r e n t i a l e-quation: 2 t x 2 0 380x 1 8 y t = ^xx o Q ( O ^ x ^ l ; 0 £ t 6 1 ) y(x,0) = x y(0,t) = 0 1 y ( l , t ) a 1+t 2 which has sol u t i o n 20 y ( * , t ) c -14t 2 Since the sol u t i o n involves a boundary layer near x - 1 i t i s reasonable to place more net points near x = 1 than else where. Table 4 shows the r e s u l t s that were obtained with a uniform mesh, and Table 5 gives the r e s u l t s - 5 2 -obtained with non-uniform mesh i n the x - d i r e c t i o n . The net spacings i n the time d i r e c t i o n were taken to be uniform i n both cases. In the case of non-uniform netspacing i n the x - d i r e c t i o n with J = 7 the net points were: x Q = 0 . 0 xl = 0.47 x 2 - 0.67 X o - 0.78 x^ = 0 . 8 5 x c = 0 . 9 1 5 x 6 -O . 9 6 = 1 . 0 0 With J - 12 the netpoints were chosen to be: x Q = 0 . 0 x ? = 0.87 x x = 0.47 x 2 = 0 . 6 1 = O .69 Xg = 0 . 9 3 4 x^ Q = 0 . 9 5 8 x,, = 0 . 9 8 0 = 0 . 9 0 x^ - 0 . 7 5 Xr = 0.80 ^1>.00 x 6 = 0.84 - 5 3 -With J = 28, 43 the net points were chosen i n a si m i l a r fashion. For t h i s problem the Box scheme gave the best r e s u l t s , and comparing Table 4 and Table 5 we note that considerable improvement i n accuracy i s obtained i f we place more net points i n the region where the solution changes most r a p i d l y . TABLE' 4 (Uniform mesh) Cubic Spline Scheme Box Scheme C. N. J e e x exx e e X e 7 * ^oo . 8 9 2 +•2 . 1 2 5 .423 + 2 . 2 6 9 ° ° . 9 0 8 ' 1 . 5 1 4 ° ° 12 00 . 2 3 7 ( 2 . 4 6 ) . 4 7 6 + 1 ( 1 . 7 9 ) .267 f 2 ( 0 . 8 5 ) • 9 9 5 " 1 ( 1 . 8 5 ) . 3 5 7 " 1 ( 1 . 7 3 ) . 1 9 2 ° ° ( 1 . 8 3 ) 228 .316* 1 ( 2 . 3 8 ) 00 . 9 0 0 ( 1 . 9 7 ) .875 + 1 ( 1 . 3 2 ) . 1 5 4 " 1 ( 2 . 2 0 ) .705°° (1.91) . 3 0 2 " 1 (2.13) 4 3 .128'1 ( 2 . 1 1 ) .370°° ,(2.07) • 4 3 7 + 1 (1.62) . 6 3 0 " 2 -(2.08) . 3 0 2 0 0 ( 1 . 9 8 ) . 1 2 6 " 1 ( 2 . 0 4 ) -54-TABLE 5 (Non-uniform mesh) Cubic Spline Scheme Box Scheme J e e X exx e e x 7 . 1 0 9 ° ° .185 + 1 . 1 1 9 + 2 .221"1 . 1 0 3 + 1 12 .258"1 . 4 9 0 ° ° .377 + 1 • 478"2 .269°° (2.66) (2.46) (2 . 1 3 ) (2.84) ( 2 . 4 9 ) 28 . 4 0 6 " 2 .738"1 .7 56 0 0 •517"3 • 4 3 8 " 1 (2.18) (2 . 2 4 ) ( 1 . 9 0 ) (2.62) (2 . 1 4 ) A3 .163'2 .2J71"1 . 3 3 0 ° ° .168"3 .171"1 (2.13) (2.33) ( 1 . 9 3 ) (2 . 6 2 ) ( 2 . 1 9 ) -55-PART I I I . A p p l i c a t i o n to the Stefan Problem. I I I . l Introduction. In t h i s section we s h a l l i n d i c a t e how our c o l l o c a t i o n method with cubic splines may be used to obtain an approximate solution to a moving boundary value problem, usually known as the -'Stefan problem. The Stefan problem a r i s e s i n ce r t a i n physical processes that are governed by a parabolic p a r t i a l d i f f e r e n t i a l equation, together with a moving boundary condition. As an example we mention the melting of i c e . We s h a l l consider a simple type of moving boundary problem, namely that of f i n d i n g s ( t ) (the moving boundary) and y(x,t) f o r which (1.1) A) y t - y x x + f(x,£) ( O ^ x ^ s ( t ) , O ^ t ^ T ) B) y(x,0) - g(x) C) yjO,t) - 0 fs(0) = 1 s(T) ? 0 D) y(s(t),t) = o E) ^ y x ( s ( t ) , t ) - i ( t ) = H(t) -/,\ — ds - 5 6 -(Condition. (1.1 B) i s the moving boundary c o n d i t i o n ) . We require s ( t ) to be s t r i c t l y decreasing and j s(t) l ^ o ^ f o r (0 ^ t •^.T). For a discussion on the existence and uniqueness of sol u t i o n to (1.1) see Friedman [S J or Cannon and H i l l £2]. Numerical so l u t i o n of the Stefan problem has been treated by for example J. Douglas, J r . and T. M. G a l l i e , J r . C Uj, who ma'de use of an i m p l i c i t (forward) f i n i t e d i f -ference scheme. Their a n a l y s i s required the f l u x a t the f i x e d boundary to be constant. Also we mention G. H. Meyer [ 12J , who considered a more general moving boundary problem. By d i s c r e t i z i n g the time-variable he reduced the Stefan problem to a sequence of free boundary value problems of ordinary d i f f e r e n t i a l equations, which are solved by conversion to i n i t i a l value problems. -57-III.2 Description of the Scheme. F i r s t we s h a l l replace the moving boundary condition (1.1.E) by an equivalent condition that does not involve the de r i v a t i v e s of the moving boundary. Since the der i v a t i v e dy/dt i s zero i n the d i r e c t i o n of the boundary s ( t ) i t follows that (2.1) s ( t ) y x ( s ( t ) , t ) z - y t ( s ( t ) , t ) r - y T X ( s ( t ) , t ) - f ( s ( t ) , t ) Using (2.1)tthe moving boundary condition (l.l.B) can now be written as (2.2) |H(t)-^ 2y x(s(t),t ) j y x ( s ( t ) , t ) . y x x ( s ( t ) , t ) + f ( s ( t ) , t ) The moving boundary condition written i n form (2.2) w i l l be used i n l a t e r discussion. Let a p a r t i t i o n 0 = xn^- x,<< T - 1 U 1 u of the unit i n t e r v a l Co,lJ be given. The ttimel&evels denoted by 0 ~tQ<^t^<£ C t n with N — J-2 are va.?yigg. Mesh spacings are given by and *n = *n * V l ( 1 ~ n - N ^ J-2) -58-At t i m e - l e v e l t the meshooints are •> ^ . , , n - I J I i = o so t h a t whan proceeding from one t i m e - l e v e l to the next one, the rightmost meshpoint w i l l be dropped. At t i m e - l e v e l t t h i s rightmost meshpoint xT i s to n . J-n be the approximation to the moving boundary. (See f i g . 2.1). t -t. J-2. •it to 1 l \ 1 1 X o so X/ -59-As i n Part I we require the approximation to the solution y(x,t) of (1.1) at each time t to be a cubic spline S n(x) s a t i s f y i n g : The boundary conditions: ( 2 . 3.A) (S X)Q = 0 ( 2 . 3.B) S j _ n = 0 and the c o l l o c a t i o n conditions : 1 / ( O ^ j ^ J - n ) The timesteps Ik ( ^ 2 are unknown before hand; r I n J n - 1 they w i l l be determined by req u i r i n g that S (x) s a t i s f y the moving boundarv condition ( 2 . 2 ) . That i s S n(x) must s a t i s f y (2.A)(H(t„) - J-*^)"^ I ( \ ) " . N - ( S ^ - 3, 5 G ( V S n ) = 0 n—1 By the r e s u l t s of part I we can, f o r given S (x) and k , e a s i l y obtain a solution to the c o l l o c a t i o n n* J scheme ( 2 . 3 ) by solving and equivalent f i n i t e difference scheme. This f i n i t e d ifference scheme again leads to matrix equations that involve only t r i d i a g o n a l matrices, and that are therefore r e a d i l y solved by the standard L_U decomposition algorithm^,pp. 35 - 42. J . -60-Appftoximations to the f i r s t two derivatives of the solution y(x,t) can be obtained by e x p l i c i t com-putation from the sol u t i o n of the f i n i t e difference equations. These approximations are needed to com-pute the' value of G(k n, S n) as defined by ( 2 . 4 ) . The basic idea of the scheme i s to f i n d k n for which the sol u t i o n 8 n ( x ) of the c o l l o c a t i o n e-quations ( 2 . 3 ) s a t i s f i e s ( 2 . 4 ) . The only place where n o n - l i n e a r i t y enters the scheme i s i n r e -l a t i o n ( 2 . 4 ) . In Part I we did not consider boundary conditions involving a d e r i v a t i v e . We s h a l l therefore f i r s t i n -dicate how to replace the boundary condition ( 2 . 3 .A) by a f i n i t e difference equation. By r e l a t i o n ( 5 « 3 ) of Part I we can write ( 2 . 3 .A) as: and the c o l l o c a t i o n condition ( 2 . 3.C) can be written as: ( 2 . 5 ) Sj-sg 6 I J n-1 -61-Using ( 2 . 6 ) i n ( 2 . 5 ) we get (2.7) B h ( S Q , S n ) 3 (-lAx " 2 V 3 k n ) S 0 + ( 1 / h 1 ~ h l f U *~ \ n-1 ^, , n - l , 0 ~ n , 0«n-l . n-1 , h J j+l For the f i n i t e difference scheme ( 2 . 8 ) to be equivalent to the c o l l o c a t i o n scheme ( 2 . 3 ) i t i s necessary, as i n Part I, to define S (x) as the cubic spline i n t e r -polant of the i n i t i a l data g(x) at the meshpoints |x. o p) j = 0 That i s S should s a t i s f y S° ~- g ( x j (0 3 3^0 For two d i s t i n c t i n i t i a l approximations k ^ ^ and k ^ ^ n n to the si z e of the next timestep k we set n (0) (0) n n-1 n arid - * , • k f 1 ' n ~ n-1 n jRHd solve the f i n i t e difference equations (2.8), (only t r i d i a g o n a l matrices are involved i n the corresponding matrix equations), to obtain two sets of sol u t i o n , name -64-( i ) s n ) J " " ( 1 * 0 . 1 ) From these we can e x p l i c i t l y compute two i n i t i a l ap-proximations to the second derivative of the solution y(x,t) at the location of the approximate boundary x _ , name l y J-n n S i m i l a r l y two i n i t i a l approximations to the second d e r i v a t i v e at the meshpoint x > that i s nearest to the approximate J-n-1 boundary x can be computed by J-n . , » n ^^ 2 , n-1 . , , n - l , (*) n-1 (2.10) (S ) - — (S -S )- S ) - f ( x T . , t )-f V ' K x x ; j - n - l " ^ K J-n-1 j-n-1 xx J-n-1 V J-n-1' n ' J ^ n _ 1 k n ( i 0,1) Relations ( 2 . 9 ) and (2.10) follow from the equivalence of the c o l l o c a t i o n scheme ( 2 . 3 ) and the f i n i t e difference scheme ( 2.8). In f a c t these r e l a t i o n s are a s l i g h t l y rewritten form of the c o l l o c a t i o n equations ( 2 . 3 £ C ) . The equivalence of the two schemes w i l l not be shown here. -65-The proof requires only a minor modification of the proof following theorem 5«3 of Section I. By making use of r e l a t i o n (5»2) in Lemma 5.2 of Section I, we can a l s o get two i n i t i a l ap-proximations to the derivative of y(x,t) at the approximate boundary l o c a t i o n x , namely J-n _ v x ' (S~ - S" ) h ( i ) ( i ) (2.11) ( S i " ^ + ± 1 < 2(S )n_ +(8 )* \ v v x'J-n ^ ~ ^ v xx'J-n v xx'J-n-1 .T—n v. h T 6 J-n ( i = 0,1) The ultimate goal i s to f i n d a zero of the function 0(1^,S n) as defined by (2.4). Since S n depends on the choice of the timestep k^ (that i s unknown beforehand) we can consider G as a function of k only. If we f i n d k , with cor-n n responding approximation S n(x) to the s o l u t i o n y(x,t) of the continuous problem (1.1) at tims-l e v e l + - n - 1 + k n , such that G(k n,S n) = 0 then S n s a t i s f i e s the moving boundary condition (2.2) -66-exactly. We donnot know as yet under what conditions a root k of G e x i s t s , however, numerical computa-n tions have indicated that such conditions w i l l not be too r e s t r i c t i v e on the solutions y(x,t) and s(t) of the defining problem (1.1), apart from the f a c t that s(t) i s assumed to be s t r i c t l y decreasing with | s ( t ) | E / | f / < n(°) n^ 1) Having obtained S and S, as above (0) (0) f 1 ) we now compute G(k n , S n v ') and » S n ) (2) where a f t e r a new approximation k to k can be 1 r n n computed by for example the secant method 10, pp. 99 - 1 0 2 J . That i s (2) (1) (1) n ( 1 ) ( k n (D. k(0) n 1^ = k -G(k ,S ) (1) n ( l ) (0) (0) n G(k n ,S n )-G(k n ,S n ) \ S i - 0 Proceeding i t e r a t i v e l y we get a sequence •(<) v ' I (m) n'"' I The i t e r a t i o n is discontinued when j G(k Q ,S n j)j I ^—- ^ f n('Mm and corresponding < S r . ^ . where £ > 0 i s a pre-assigned small number. We then set -67-k _ k n — n (m) and S = S J j (0 < j ^ J - n ) The f i n a l approximations to the f i r s t two derivatives at the meshpoints / x . 7 ^ n [ 3 j j 0 time l e v e l t n are again determined by at ( 2 . 1 2 ) ( S )" ~- Us*-**'1) - ( S J " - 1 - f n - f " " 1 (0 - j - J-n -2) x x J k n 3 3 xx j j j >n Sn -S n h ( 2 . 1 3 ) (S )" - J+1 J J + 1 x 3 — V i 6 1 2 < S x x ) > ( S ^ , x x 3 xx j+l ( O ^ j ^ J - n - 1 ) n J m ) n n ( m ) n n ( m ) ( 2.H) ( S _ ) T , - ( S ) n : (S f - (S ) ; (S ) . - ( S ) v x x J - n - l - v x x ' j - n - l v xx'j-n v xx J-n v x'J-n v x ' j - n Having obtained the approximations S and x T to y(x,t ) and $ ( t ) we can now proceed to the j —n n next time level. III.3 Numerical Results. The f i r s t problem i s to solve y t = y x x (° ^ $ ( t ) ; o ±t<~ y ( x , 0 ) - 1-x2 y x(o,t) = o y ( s ( t ) , t ) = o y ^ s ( t ) , t ) - S ( t ) - At-i which has 2 y(x,t) - l - 2 t - x and s ( t ) z. l / l - 2 t as i t s unique s o l u t i o n . The mesh in x - d i r e c t i o n was taken uniform. For t h i s problem the numerical scheme yi e l d s the exact s o l u t i o n . That i s S n , (S ) n , (S ) n and x agree with J x J x x j J-n y(x , t ), y (x , t ), y (x , t ) and s ( t ) res-j n X j n xx J n n pectively almost up to machine accuracy. ( 0 ^ j ^ J-n, J7-2, l ^ n ^ J - 2 ) . Also i f we l e t S n(x) be the cubic spline interpolant of J .^ n _ s a t i s f y i n g ( 2 . 6 ) at x and x T then j 3 \ 3 - 0 0 J-n -6 9-w i l l be exactly equal to y(x,t ) and x - s(t ). n J-n n (This r e s u l t cannot of course, be expected i n general. In 2 t h i s problem the simple solution i s y(x,t) =• l-2t-x , which for f i x e d t i s a polynomial of degree two i n x, and for f i x e d x a l i n e a r function i n t ) . The second problem i s (-x 2-2t-l) y - y + —• — * " x x ( O ^ x ^ s ( t ) , 0 ^ t ^ | ) y(x,o) y x(o,t) y ( s ( t ) , t ) y x ( s ( t i ) , t ) - s ( t ) x 2 - i o 0 -2 - 4t t IT l-2t ' which has y(x,t) = \f l - 2 t ' ( x 2 4 2 t - l ) and s(t) - if l - 2 t ' as i t s unique sol u t i o n . We note that ds | —> oO as t dt The mesh in x - d i r e c t i o n was taken uniform i . e . h . = x . S x - h = -J 0 j-l J (1 — j — J) -70-An approximate solution was obtained with various values fo r J, up to time l e v e l t with N = J-2. The r e s u l t s N are given i n Table 1 with e = max I If 1 - y(x , t ) / 6 x - - ^ l ^ - ^ j ' V l exx = * « I - ^xy^l e - max x t - s(t ) I 3 J J-n v n / where the maximum i s taken over 1 ^-n J) with N( J) the largest integer such that t ^ ^ . 4 6 . For e (the maximum error i n the approximation to the moving boundary) the maximum i s taken over 0 — n — N(J). The quantity between brackets i s the observed rate of convergence. -71-Table 1 e e e e„ J X XX s as"2 . 6 7 ~ 2 • 5 9 ~ 2 . 10 .47-3 .IB" 2 . 1 1 " 1 . 20 (1-95) ( 2 . 0 ) ( 2 . 0 ) ( 1 . 9 0 ) .12-3 . 3 9 ~ 3 . 2 9 - 2 . 4 2 " 3 40 (2.0) ( 2 . 0 ) ( 1 . 9 5 ) ( 1 . 9 5 ) At each time t ^ we need two i n i t i a l approxim-ations to k . A good approximation can be derived from ( l . l . B ) , namely by computing r - + 2 t » j £ * - h< vi> ( y i s the approximation to s ( t n _ ^ ) ) a n d then setting ( 0 ) n - - h J-n 41 r ( 1 ) (o) k Sail be taken as a f r a c t i o n of k , for example ( 1 ) ( 0 ) take k n = 0 . 9 5 k n . In the problem above four to -72-f i v e a d d i t i o n a l i t e r a t i o n s per timestep were s u f f i c i e n t . (The value of £r was taken to be 1 0 " ^ . ) For given an approximate s o l u t i o n was com-puted for K - J - 2 time l e v e l s . For larger values of J the c r i t i c a l time t - |-, where the slope of the moving, boundary becomes i n f i n i t e , i s approached more c l o s e l y . Nevertheless a --ood accuracy was attained at t i n e t . J - 2 ' a.s can be seen in the table below. Table 2 t J - 2 * J - 2 e *3-2 * J - 2 exx t J - 2 e . 3 J . 4 3 1 5 4 4 as"2 . I I " 1 . 8 6 " 1 . 7 9 * 10 . 4 9 5 3 9 3 . 3 8 " 2 20 . 4 9 3 3 5 0 . i r 2 . 3 2 " " 1 _2 . 2 0 * 4 0 - 7 4 -- p • 73 de&s no'l exist] Bibliography [l J J. H. Ahlberg, E. N. Nilson, and J. L. Walsh; The theory of splines and t h e i r application; Academic Press; N. Y.; 1967. [_2 J J. R. Cannon and C. D. H i l l ; Existence, uniqueness, s t a b i l i t y and monotone dependence i n a Stefan problem forrthe heat equation; J. Math. Meoh.; 1 - 2 0 ; 1 9 6 7 . (_3 J C. de Boor and B. Swartz; Collocation at Gaussian points. (To appear i n SIAM J. Num. Anal.) £4 I J. Douglas J r . and T. M. Gal l i e J r . ; On the numerical integration of a parabolic d i f -f e r e n t i a l equation subject to a moving boundary condition; Duke Math. J. 2 2 ; 557 - 571; 1 9 5 5 . ["5 J J. Douglas J r . and T. Dupont; A f i n i t e element collocation method for non linear parabolic e-quations. (Manuscript). [& ] J. Douglas J r . ; Survey of numerical methods for parabolic d i f f e r e n t i a l equations; Advances i n Computers; Vol. 2 ; Academic Press; N. Y.; 1 9 6 1 . (7 ] K. W. Fong; Numerical solution of parabolic e-quations by the Box scheme; Thesis; C a l i f o r n i a Institute of Technology; 1 9 7 3 . - 7 5 -A. Friedman; P a r t i a l d i f f e r e n t i a l equations of parabolic type; Prentice H a l l ; I 9 6 4 . T. N. E. G r e v i l l e ; edJ;Theory and applications of spline functions; Academic Press; I 9 6 9 . E. Isaacson and H. B. K e l l e r ; Analysis of numerical methods; J. Wiley and Sons; N. Y.; 1 9 6 6 . H. B. K e l l e r ; A new differ e n c e scheme for parabolic problems; Numerical solutions of p a r t i a l d i f f e r e n t i a l equations; V o l . 2 ; J. Bramble, edl; Academic Press; N. Y.; 1971; pp. 327 - 3 5 0 . G. H. Meyer; On a free interface problem f o r l i n e a r ordinary d i f f e r e n t i a l equations and the one-phase Stefan problem; Numerische Mathematik 16; 2 4 8 - 267; 1 9 7 0 . R. D. Russell and L. F. ^hampine; A c o l l o c a t i o n method f o r boundary value problems; Numerische Mathematik; 1 - 28; 1 9 7 2 . M. H. Schultz; Spline Analysis; chapter 4 ; Prentice H a l l ; 1 9 7 3 . ~~