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<t 4T) -5-We s h a l l make the following assumptions: p,q and f are continuous functions of x and t , and p i s three times .continuously d i f f e r e n t i a b l e with respect to t i n R(T) - [ 0 , l ] x [ 0 , T J . The equation (1.1) must be parabolic i n the sense of petrowski, i . e . P ( * » t ) ^ p ^ O for (x,t)£ R ( T ) F i n a l l y we assume that there i s a unique solution to the d i f f e r e n t i a l equation (1.1) subject to the i n i t i a l and boundary conditions (1.2)., and that t h i s solution i s three times continuously d i f f e r e n t i a b l e with respect to t and four times continuously d i f f e r e n t i a b l e with respect to x« On the rectanble R(T) nlace a net R^T) such that x 0 » 0 X j = 1 t 0 i 0 t N ^ T with netspacings We s h a l l use the following notation: -6-For functions defined on H (T) h ^ =?(x.,t ) j J n For functions defined on H(T) y > 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) 1 „n , n S b i n °1 0 n a i2 n b2 n °2 o n n n a^ T i b c • J - l u j _ l J - l <n> , n n u = (U Q , u , n T and (») n n — = <V V n T • V with n _ h.P . . 1 (2.3) = 3k h . n J J - l /h . * h . . j 3k h . h . , J J+l - 1 3 -h j + l q j * l J = 3k n h J * l d = 6 ( t ) 0 0 n y d J = ^ ( t n ) and .h.-f' 1 1 v ' - i y - 1 ~ 3 k — + i r + — 5 — J'" 1 2(h +h ) ^ n - 1 1 1 n - 1 J j+l j \ 3k n h h U n - 1 n-1 / h "-P. i 1 + n -1 \3k n j+l U n -1 j+l + i IK ( f n , + f n _ 1 ) +2 (h. * h. j ( f n • 6 { j v j - l j - l ' v J j * l n j „n-l. , , n n -1 f ) *h . , ( f + f () j J+l j+l j + l ' - U -Theorem 2.1 For each n ( l ^ n ^ N ) there i s a unique solution to the matrix equation (2.2) provided that ^4 n* and n _n where n* J n j-q = max j | and n n p + = min (p ) j J Proof. It i s s u f f i c i e n t to show that f o r each n (I4 n^; N) the matrix A^ n^ of (2.2) i s diagonally dominant, provided (2.4) holds. S p e c i f i c a l l y we want to show that - 1 5 -n 1 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 <Sxx) j+ J+l 6 h +h. J J+l h.+h. , J j+l ( S x x ) j + l ^ n l n l h j + 1 — ( sxx ) r.i^ 2 ( ^ r 1 ^ — n-ll h .+h J j+l h >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<n^N) to the collocation equations (1.3> 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 <S*>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 , <fn .n-1 , 6 I k k n n The f i n i t e difference scheme that i s equivalent to the c o l l o c a t i o n scheme (2.3) i s now given by (2.8.A) I^S* z 0 ( 1 - ^ j ^ r J - n - l ) (2.8.B) Bh(s£, S*) z 0 (2.8.C) S n - 0 J-n , n n. , where B (S ,S ) i s defined by (2.7) and h 0 1 - 6 2 -L h s j - - < h . J ( n n-1 's. -s. J - l J - l h .fh \ J j+l \ n S -S \ +• 2 J J j+l ' j + l j+l V k / h .4h \ n ' J J+l n / ^s n -sn sn-sn s11"1^11-1 s""1^31"1 j+l j j j ; j l l j+l j j j - l - < 2 h . J J+l 4-* ( Y H j + i } j+l h >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 <L j ^ J) J J with endpoint conditions xx 0 0 xx' and (S ) = g (1) xx J x x - 6 3 -P r a c t i c a l computations have shown that the accuracy of the scheme i s maintained and i n c e r t a i n cases even im-proved i f instead we define (S )^ by xx j (S ) % g (x ) ( O i j i J ) ] x x J xx 3 J The procedure for obtaining an approximate solution to the problem (1.1) i s now as follows. Sup-pose that the approximations to y(x,t) and y (x,t) at xx time t ^ are given by ' n - l ) J-n+1 J f, .n-1 ) J-n+1 H i : : : 1 - «-> 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 ) <X> 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 .
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Numerical solution of linear second order parabolic...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Numerical solution of linear second order parabolic partial differential equations by the methods of… Doedel, Eusebius Jacobus 1973
pdf
Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
Page Metadata
Item Metadata
Title | Numerical solution of linear second order parabolic partial differential equations by the methods of collacation with cubic splines |
Creator |
Doedel, Eusebius Jacobus |
Publisher | University of British Columbia |
Date Issued | 1973 |
Description | Collocation with cubic splines is used as a method for solving Linear second order parabolic partial differential equations. The collocation method is shown to be equivalent to a finite difference method that is consistent with the differential equation and stable in the sense of Von Neumann. Results of numerical computations are given, as well as an application of the method to a moving boundary problem for the heat equation. |
Subject |
Differential equations, Parabolic Differential equations, Linear -- Numerical solution Spline theory |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2011-03-15 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0080130 |
URI | http://hdl.handle.net/2429/32482 |
Degree |
Master of Science - MSc |
Program |
Mathematics |
Affiliation |
Science, Faculty of Mathematics, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-UBC_1973_A6_7 D63_3.pdf [ 2.72MB ]
- Metadata
- JSON: 831-1.0080130.json
- JSON-LD: 831-1.0080130-ld.json
- RDF/XML (Pretty): 831-1.0080130-rdf.xml
- RDF/JSON: 831-1.0080130-rdf.json
- Turtle: 831-1.0080130-turtle.txt
- N-Triples: 831-1.0080130-rdf-ntriples.txt
- Original Record: 831-1.0080130-source.json
- Full Text
- 831-1.0080130-fulltext.txt
- Citation
- 831-1.0080130.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
data-media="{[{embed.selectedMedia}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0080130/manifest