Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Numerical algorithms for the solution of a single phase one-dimensional Stefan problem Milinazzo, Fausto 1974-01-22

You don't seem to have a PDF reader installed, try download the pdf

Item Metadata


831-UBC_1974_A1 M53_8.pdf [ 5.63MB ]
JSON: 831-1.0099968.json
JSON-LD: 831-1.0099968-ld.json
RDF/XML (Pretty): 831-1.0099968-rdf.xml
RDF/JSON: 831-1.0099968-rdf.json
Turtle: 831-1.0099968-turtle.txt
N-Triples: 831-1.0099968-rdf-ntriples.txt
Original Record: 831-1.0099968-source.json
Full Text

Full Text

NUMERICAL ALGORITHMS FOR THE SOLUTION OF A SINGLE PHASE ONE-DIMENSIONAL STEFAN PROBLEM by Fausto Milinazzo B.Sc, University of British Columbia, 1970 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in the Institute of Applied Mathematics rtrvd Statistics We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH April, 1974 COLUMBIA In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the Head of my Department nr by his representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of Aymlipr) I" a -h>iPTr);T+. -j n« rmrl fit.gtistins The University of British Columbia Vancouver 8, Canada Date June I>. ~i 97/ ABSTRACT . A one-dimensional, single phase Stefan Problem is considered. This problem is shown to have a unique solution which depends continuously on the boundary data. In addition two algorithms are formulated for its approximate numerical solution. ' . The first algorithm (the Similarity Algorithm), which is based on Similarity, is shown to converge with order of convergence between one half and one. Moreover, numerical examples illustrating various aspects of this algorithm are presented. In particular, modifications to the algorithm which are suggested by the proof of convergence are shown to improve the numerical results significantly. Furthermore, a brief comparison is made between the algorithm and a well-known difference scheme. The second algorithm (a Collocation Scheme) results from an attempt to reduce the problem to a set of ordinary differential equations. It is observed that this set of ordinary differential equations is stiff. Moreover, numerical examples indicate that this is a high order scheme capable of achieving very accurate approximations. It is observed that th (ii) apparent stiffness of the system of ordinary differential equations renders this second algorithm relatively inefficient. TABLE- OF CONTENTS ABSTRACT TABLE OF CONTENTS LIST OF TABLES LIST OF FIGURES ACKNOWLEDGMENTS < . INTRODUCTION CHAPTER I The Problem Existence and Uniqueness of the Solution The Numerical Stefan Problem Non-dimensionalizing EXISTENCE AND UNIQUENESS The Integral Equations Existence and Uniqueness of the Solution for Small Time CHAPTER II THE SIMILARITY METHOD Lie Group of Transformations Invariance (i) (iii) (vi) (vii) (viii) 1 2 4 5 7 11 11 18 Existence and Uniqueness of the Solution for All t 6 OiTl 29 32 33 35 The Most General m-Parameter Lie Group of Transformations Leaving Invariant (2.1) 37 Reducing the Number of Variables 42 The,Classical Group of the Heat Equation 44 (iii) (iv) CHAPTER III A USEFUL SIMILARITY SOLUTION OF THE HEAT EQUATION FOR.THE STEFAN PROBLEM 46 CHAPTER IV THE SIMILARITY ALGORITHM 60 CHAPTER V THE CONVERGENCE OF THE SIMILARITY ALGORITHM 65 Continuous Dependence for Small Time 65 Continuous Dependence for All Time 80 Convergence of the Similarity Algorithm 82 Order of Convergence 83 CHAPTER VI THE SIMILARITY ALGORITHM - NUMERICAL RESULTS 86 : Numerical Examples 8Optimization of the Similarity Algorithm 95 Order of Convergence of the Similarity - Algorithm 98 The Small Time Versus the Large Time Representation of the Similarity Solution 102 Comparison of the Similarity Algorithm with Lotkin's Difference Scheme 112 CHAPTER VII A COLLOCATION SCHEME 115 The Lagrangian Equations for Heat Conduction 115 A Galerkin Scheme 127 A Collocation Scheme 135 Numerical Results 142 CHAPTER VIII CONCLUSIONS 150 BIBLIOGRAPHY 152 APPENDIX A . APPENDIX B APPENDIX C APPENDIX D APPENDIX E APPENDIX F (v) 156 162 167 171 174 179 LIST OF TABLES TABLE 6.0 Errors Versus Time Increment (Boundary Data (6.0)) TABLE 6.1 Error Versus Time Increment Ai. Using Modifications I and II TABLE 6.2 Observed Order of Convergence TABLE 6.3 Approximate Operation Count TABLE 6.4 Similarity Algorithm Versus Lotkin's Difference Scheme TABLE 7.0 Errors in UU,t),W^'1t\Mlw(«1t)and JCtV (Exact Solution (7.25)) TABLE 7.1 Errors in v*w,t>, u. and Jtt) (Exact Solution (7.26)) 88 97 101 109 1.13 146 147 (vi) LIST OF FIGURES FIGURE 0.1 FIGURE 1.0 FIGURE 4.0 FIGURE 6.0 FIGURE! 6.1 FIGURE 6.2 FIGURE 6.3 FIGURE 6.4 A Melting Slab 3 The (weak) Maximum Principle on a Non-rectangular Domain 13 The Similarity Algorithm 62 Approximate Temperature Distributions at T =.40 for the Boundary Data (6.0) 90 Approximations to the Position of the Boundary >5C-i)up to T=.4 for the Boundary 91 Data (6.0) Comparing H(10 and Vi("t) for the Boundary Data (6.0) 92 The Approximate Temperature Distribution for "t Between 0 and .1 for the Boundary Data (6.1) 94 Comparing the Given Heat Flux with that Generated by the Similarity Algorithm Using Modification II 97 (vii) ACKNOWLEDGMENTS I am very grateful to my supervisor, Dr. George Bluman, for having suggested the topic and for his. patient, encouraging guidance during the course of the work. I am grateful to Dr. James Varah for many hours of helpful and useful discussions. I thank the National Research Council of Canada and the University of British Columbia Mathematics Department for their financial assistance. I also express my appreciation to my wife, Beverly, for carefully typing the manuscript. Finally, I thank my friends in the Institute of Applied Mathematics and Statistics, particularly Hart Katz and Eusebius Doedel; for each one has contributed significantly to whatever I have achieved. (viii) INTRODUCTION This Thesis is concerned with the diffusion of heat through a medium which is experiencing a change of phase. Characteristically such problems involve a "moving" surface made up of points at phase changes to another. If the position of the surface is given as a function of time, the problem, known as the Inverse Stefan Problem, reduces essentially to one of solving a parabolic differential equation with associated boundary conditions on an irregular domain. Evidently if the differential system is linear then so is the Inverse Stefan Problem. However, when the position of this surface is not given a priori, the problem, referred to as a Direct Stefan or Free -Boundary Problem,, becomes one of finding simultaneously the temperature distribution of the medium and the position of the "moving" surface. As can be seen readily the Direct Stefan Problem is non-linear. Although Free Boundary Problems date back to a work of « G. Lame and B.P. Clapeyron published in 1831 and to several - 1 -- 2 -papers of J. Stefan which appeared in 1889, not until the nineteen thirties did work on such problems begin in earnest. During the past twenty years a considerable amount has been published documenting the analytic properties of one-dimensional Stefan Problems. In addition, a number of schemes have been developed for their numerical solution. The Problem We consider a particular one-dimensional single phase Stefan Problem. More precisely, we wish to describe the melting of a homogeneous slab, which initially occupies the space between ^= 9 and ^* £Q , and whose initial temperature distribution is 3© ^) • We assume that the temperature distribution, ^\ytZ) > obeys the heat equation interior to the slab for T>o . Furthermore, we assume the slab to be insulated at o , while at > the position of the "moving" boundary at time 2" > a heat flux H0(f) causes an isothermal phase change. By having the melt removed immediately upon formation, we restrict our attention to the solid phase only. The following equations govern the temperature distribution in such a slab. - 3 (1) J-Crf Cr;, r) J 7-^ JyiOyZ) - O recoup), r t to,TF), W»(r> (0.1) r« Co.Tp), (0.1a) (0.1b) (0.1c) (O.ld) (O.le) Fig. 0.1 A Melting Slab (1) Here XT iV.r), J>, (y,*(y,*) denote partial derivatives of JT(y;r) and J'ti d_ ['jS'tJ:)] - 4 -Here, - the melting temperature, c " tne specific heat, p - the density, H - the conductivity, and .\» - the latent heat of fusion are characteristic of the material and are taken to be constant. Moreover, we restrict our attention to the case X * 3-**. , H„(T)fc6. Although the general problem (0.1,a,b,c,d,e) can be dealt with numerically, throughout the analysis we will assume that the heat flux, V-i6(r) » i-s sufficient to maintain melting, i.e. .J'Gs'Or), ~) » tne temperature of the slab at the melting boundary is never allowed to fall below vTwv ~ the melting temperature. Whether for arbitrary J*0 and M6 this can be guaranteed a priori will be discussed briefly at a later time. Under the above assumption, condition (O.le) becomes the Stefan or Free Boundary condition and determines the allocation of available energy to the diffusion and the melting processes. The complete solution of (0.1,a,b,c,d,e) is then the pair of functions (T^Z), »SHr>) , Existence and Uniqueness of the Solution The existence and uniqueness of solutions to Stefan Problems are established by one of several methods. Usually the solution is expressed in terms of a set of coupled Volterra - 5 - • . Integral Equations, then the proof proceeds by either using the Maximum Principle or a fixed point argument. Cannon and Denson Hill £ 7 ]J use the Strong Maximum Principle together with a retarded argument approach to establish the existence and uniqueness of the solution of the problem they consider. Priedman £ 17 3» refining the work of Rubinstein |~ 29 ]],' treats the same problem using a fixed point argument. ' Using methods as outlined in Friedman f 17 1 we establish the existence, uniqueness (Chapter I) and continuous dependence (Chapter V) on the boundary data y H» (r^ J"0 (*)| of the solution to the system of equations (0.1,a,b,c,d,e). The Convergence of the Similarity Algorithm (see Chapter IV) then follows from the continuous dependence of (0 .1 ,a,b ,c ,d,e) on Vl0l^» The Numerical Stefan Problem Usually numerical schemes dealing with Free Boundary Problems are particular to the boundary conditions being considered. For instance, the finite difference scheme developed by Douglas and Gallie £ 11 ] uses two boundary conditions to establish an iteration which at each step in time locates the position of the "moving" boundary. Similarly the continuous methods of Mason and Farkas £ 24 J rely on the appearance of ^{X) twice in the system of equations so that again an iteration to the solution can be - 6 -established. Several authors, following the lead of Landau [ 22 ], make the transformation Xs ¥/jS'ltr) , then construct approximating schemes for the resulting system on the fixed space interval £o,lJ. For example Lotkin £ 23 ~\ uses this transformation to obtain a finite difference approximation for (0.1 ,a ,b ,c ,d , e) . One alternative to difference schemes on the fixed space interval [o,lJ has been to reduce the Stefan Problem to a countable set of ordinary differential equations. This was first achieved by Melamed £ 25 J by expressing the temperature distribution, "J(X#<r) as an appropriate Fourier Series with time dependent coefficients. The system (0.1, a ,b,,c ,d , e) then yields a set of ordinary differential equations for the Fourier Coefficients and the position of the boundary. We propose several schemes. The first scheme, referred to as the Similarity Algorithm (Chapter IV) is based on an exact solution of the Inverse Stefan Problem obtained through the Similarity Method (Chapter II). That is, solutions of the Inverse Stefan Problem are pieced together in such a way as to given an approximate solution of the Direct Stefan Problem. In Chapter VI we give numerical examples illustrating the Similarity Algorithm. ~ - 7 -The second and third schemes (Chapter VII) are closely related, and arise from attempts to reduce the Direct Stefan Problem (0.1,a,b,c,d,e) to a countable system of ordinary differential equations. Hence they can be looked upon as extensions of the method of Melamed. However, instead of taking as a basis, functions which are global on £o,l] (such as the Trigonometric functions) we adopt a finite element approach. That is, we approximate vT(v'T^ by a finite linear combination (with time dependent coefficients) of functions which have support in a subinterval (finite element) of £o,l]. A system of differential equations for the coefficients and rt'{Z) can be obtained in several ways by using equations (0.1,e). We will derive two systems of equations. The first is known as a continuous Galerkin system while the second is referred to as a Collocation system. Some numerical results are given. Non-dimensionalizing Before proceeding further we non-dimensionalize the system of equations (0.1,a,b,c,d,e) by introducing the following - 8 -variables: (2) X= >Vet , where " a " is a characteristic length T pea* T> ucx.tvi CTo?r)- TvAi/Jv^ , (0.2) •T - K T Substituting the variables (0.2) into (0.1,a,b,c,d,e) ve obtain v The characteristic length "a" can be taken to be the initial length of the slab - 9 <*,*V= MfcU,*) O < X < 4 Ct) < JC«Js b I « C o,T^ J J (T)- b„>o (0.3) i<t>* o U (ojlj (0.3a) UUID^J-O t6(o,T), (0.3b) MRt«>*)-© .' tt (OjT), (0.3c) Utx.o)-. X t C«, (0.3d) it!) s «* |M«UU^-0- h(t\J t«Co,TK (0.3e) We note that (0.3e) with j(o); b can be written equivalently as (0.3f) We now seek the solution ( -u. (.*, i), j ct>) to the system of equations (0.3,a,b,c,d,e) . More precisely, we take V.U') e C'Co.b] with v.lOi© on [0,b], A0(»)= U0 ( b)r o and \-»Ct) continuous but for a finite number of jump .discontinuities on [o,l]J. Then we look for a solution (*u. j ct)). of (0.3,a,b,c,d,e) satisfying the conditions: (a) Mxxt*,t), M lw,t) 6 C t>,.scu), t <s <O,T); (b) v.(.x,t^ <S. cro^ct)! , t C Co.T); (c) U„ C*,t> € CCo.iti^ te Co,T); We begin in Chapter I by showing that the system equations(0.3,a,b,c,d,e) has a unique solution. CHAPTER I EXISTENCE AND UNIQUENESS In this chapter we show that the system of equations (0.3 ,a ,b ,c ,d ,e) has a unique solution (14,5) for all , <£) s j U,t) i o< x < 4<tl , o < t < T To accomplish this we follow the lead of Friedman [17 J, by constructing an equivalent system of Coupled Volterra Integral Equations and showing that there exists a C>© such that for all "K< 9- , these integral equations have a unique solution. We then show that this procedure can be repeated to yield existence and uniqueness of the solution of the system (0.3,a,b,c,d,e) for the interval of time (o>T). The Integral Equations .. Before constructing the integral equations, we state several useful lemmas. The first, due to Friedman [ 17 ], is a working lemma used extensively throughout the construction of the integral equations; while the other two establish properties of -SCt.^ (the position of the free boundary) and vt>x (..St-U^-fc^ (the - 11 -- 12 -amount''of. heat allocated to the diffusion process) respectively. Defining K ix,-t;^tt) to be the usual source solution of the heat equation, that is, - Cx-g)* we have the following lemma. Lemma 1.1. Let <fit^i -set) be continuous functions on the interval .In addition, let jet) satisfy the Lipschitz condition for some constant M . Then for all "k £ i0,0"! we have sctw© Ox ^ Proof. The proof, given in Appendix A, consists of showing - 13 -r1 In order to establish the next two lemmas we use the following auxiliary propositions whose proofs are given in Appendix B. Proposition 1.1 (The (weak) Maximum Principle). Suppose "UAx^) satisfies with <St-t) a positive continuous function and -uCx,-t> e C ) , 'M„U,t>, vi^ <.«,*> £ C (•& V <8X ) where .£) is the closure of ^ and Fig. 1.0 The (weak) Maximum Principle on  a Non-rectangular Domain - 14 -then v.U,t) assumes its maximum and minimum values on the data boundary - <Sy • Proposition 1.2 (The necessary condition for melting). If (-u,s) is a solution of the system (0.3,a,b,c ,d ,e) then Vj.uui.t^ jo for all t e CO,T3 . We are now ready to establish the following properties. Lemma 1.2. (0.3,a,b,c,d,e) then .SCiV satisfies the Lipschitz condition If SC*.) is a solution of the system of equations (1) for all t., ,tj. €. O.T] . Proof. Consider condition (0.3f) « o at t, and i i Then Notation: lt*il,-where o o Since A' - - H, (Jlt.+j),!,^) and /Wo) = o , Proposition 1.2 implies that ^^^5 > <3*°* Tnus t, and hence ••Uct,..>-s(i,)| > «<*«miT » t(-t,l . Lemma 1.3. If ("U,5) is a solution of the system of equations (0.3,a,b,c,d ,e) then 'M* C*,tV satisfies for all (x,t} C c5b , provided vt^t^-O satisfies the hypotheses of Proposition 1.1. Proof: trivial. Having established these preliminary results we can construct the integral equations which will ultimately allow us to establish existence and uniqueness of the solution to the system of equations (0 »3 ,a ,b ,c ,d ,e) . To this end. we introduce the Green's functions for the half plane 6*Ci,t;5,T)t K(*,ti^t)+ KU,i;-5,:r), - 16 -and note that any solution v.C^,t) of the heat equation satisfies Green's Identity §LfG+£M - ^ 3 6+l - &L f(5*vl =6 (1.0) Sk 5fl- art J in the domain • O8£ = JC?,^- o<^,<s(rv> o<f<T<t-f^ where £ >o. Integrating (1.0) over and using conditions (0.3b,c) and 6^ t*, i; o,2*V~ O we obtain r° O- ^ vt(^,£) G + CKlt;^<r)c|<F .sen + \ vUE;,t-0 G + (*,t;^, W^g; (1.1) In Appendix C we show that as f->o (1.1) becomes (1.2) - 17 -Differentiating (1.2) with respect to x and applying Lemma 1.1 as x-»J<-k)-o we find that -o l-O 2 -un c 5 tti satisfies •--vm* * .}*0<3>G£Utti,t;e,o)d£ (1.3) t ^ vcr> G* W(t>,ti3(*^r) dr. Since <S^.(*,i>?to)*- <»*<*•,t;£,o) , G~(x, t • o,o)s o and -u0t»=o, after making the appropriate substitution, we integrate by parts the first integral expression of (1.3) to obtain ° (1.4) it • © Moreover, as does Friedman £ 17 } we integrate the condition (0.3e) and find that .sm- b + \ (v(rl-h(r)) dt. (1.5) © We have that W-t) satisfies (1.4) where s(t) is given by (1.5), hence we refer to v(t) as the solution of (1.4), (1.5). Furthermore, we have the following equivalence between v(t)(the solution, of (1.4), (1.5)) and (v,s) (the solution of (0.3 ,a,b,c ,d ,e)) - 18 -Lemma 1.4 (The equivalence of the differential and integral sys terns). If v(t) is a solution of (1.4), where s(t> is given by (1.5), then (v,s) (v<(*,t) defined by (1.2) and jti) defined by (1.5)) forms a solution of (0.3,a,b,c,d,e). Conversely if (Mts) is a solution of (0.3,a,b,c,d,e), then \> (t> * H» Utt),t) is a solution of (1.4). Existence and Uniqueness of the Solution for Small Time From the equivalence of the system of differential equations (0.3,a,b,c ,d,e) and thesystem of coupled integral equations ((1.4), (1.5)), we see that showing that the former has a unique solution reduces to demonstrating that the latter has a unique solution. To establish that ((1.4), (1.5)) has a unique solution we make the following definitions. Definition 1.1: The proof is standard and is given in Appendix D. the set of bounded continuous functions on |0,tf" Definition 1.2: the closed M-sphere in Cr . - 19 -Definition 1.3. Define • H"* to be the transformation given, by (1.4), (1.5), that is, + 2 \vt*v G;csto,t;5«*i,t;H-t where It is easy ,to see that Moreover, we have the following theorem. Theorem 1.1. There exists a <r>o such that Hf1 : C^M -> CO-.M where Proof. - 20 -Suppose V« CO-,M then Uvll^ x< M and hence Let satisfy the following inequality 0" 5 "<v^ L~N-Since from (1.4) we conclude that Writing (1.6) where N*\\WI|T. Then for all 1 C- [o^] which in turn implies that v< *ft) t 1 b for all **<r- (1.7) \\T^1| * aniA + M^^^^M^ • . (1.8) where G, = \\K - 21 -we estimate &i and Ga in turn. Noting that J(t) is Lipschitz continuous we find that \ G,l i -?L? ,(^^N) -t"« . To estimate S, we use (1.7) to obtain (1.9) . A(2) 16,1 «: a «*&e(-^.i) « (1.10) we Now applying the inequality etftcfO,* — • JL to (1.10) have IGJ* J. / t» . (1.11) Combining (1.9) and (1.11) we see that (1.8) becomes HTvllt f 2»*0rt'b + M (o('(M4W)+ J.) tVa. ty'/j b Hence the conclusion of the theorem follows if we insist that o-also satisfy the inequality C- J X ! . (1.12) (2) „ , t \ -** v ' We use the notation erucvM- — Vc At , - 22 -The following theorem shows that we can further restrict the size of <r so that T1 is a contraction mapping on ^C,M and hence allows us to conclude that (1.4), (1.5) has a unique solution in C for a small time. Theorem 1.2. There exists a C>o such that is a contraction mapping on Cg-|(^ for all t 6 Co, «rl,. Proof. Initially let c be such that <r i «r«0 % where ce satisfies (1.6) and (1.12). If ^ 1 it), v> Ct> C- CVtf^ let V< fc>, 5 Ct> satisfy (1.5) with -o'Ct) and -JCt) respectively and define Since t-fe) ^ v'tt^ <c we have f« for all • . From (1.5) we have the following inequalities: | sit) - j'(t\V * e «•* (1.13) (1.14) (1.15) - 23 -and as before V $ 5lkV, S»<*>-.$ 3 b. for all r,t « Z*,<*1 . (1.16) Now consider TvTv'= v, - V, where We can write V, - V,' + V," where V> 2 ^o(%)[wCS(^^o)-W(5'(t)^;^Jo)]a^ % o v;.,s-2^v.^)[KUrt),trf,o)-K(s,a)Jti-^o)] d|. o Applying the Mean Value Theorem and the inequality (1.16) to V,' we obtain IVI * £ (1.17) To estimate V,' we assume that .s'ct) > .S ("fc) and consider the possible cases: 24 Case I: o< b * AU>< V<t > i |b> Case II: o < h * * x< * | b , Case III: o < | * 5It) < «'<t) .< b . Considering v;s^^C5>[K_<SW,t3fIo>-KCS,<i)>t;f,o)]e»5 in Case I, in Case II, © Jt-fc) in Case III. and applying the Mean Value Theorem an appropriate number of times, we arrive at the estimates: Case I: IV.'W i. II • «t't "* Case II: W\\ J^1t"* ^ (1.18) Case III: I v'. I * ^ •<* t"' Combining (1.18) with (1.17) we see that W,l S tlfM^A *'V (1.19) where 25 To estimate we write V* = ^[\vl*> ^l±i+i£*l> K(s<t),t--s «->,*> dr o *Ct-r) Since Jit) is Lipschitz continuous we see that 17 "* (1.20) Applying the Mean Value Theorem to I (w> (i"*> J i we obtain - 26 -lw*J $ t"* . (1.21) ft"*-Writing Wj as _ (sttwwf w . _1_ \ -u'trvU'ttr- &'c^) e «<fc-*> fi we see that the last term can be estimated as follows \ (S'CO- -S'ttV)*- <SU>- 3 < 2.'(M+A/) S U'til -Hence taking a- to further satisfy 3 M (ot^ CM +//) o- J i (1.22) and using the inequalities I I ~ «~&| $ 3l^l. (l^lj i) and (1.15) we find that Wv,» i 3M (+* + t3'2 . (1.23) To complete the estimate of Vg , we write as - 27 -where L = 9 \ (•-Qtt)-i>'<*»(.-SC-fc.)-«-3(tn e dt The estimation of U, involves a straightforward application of (1.7), (1.10) and yields * 3£.t>, (1.24) Tr"a b To obtain an estimate for L2, the Mean Value Theorem for a function of two variables must be applied to the function C*-*^)e (a-any non-zero constant). A simple calculation then leads to the estimate VLJ < »*« (1.25) Using (1.25) and (1.24) we see that |V'| < £ (lVM*'+3) t''1 (1.26) Hence (1.20), (1.21), (1.23) and (1.26) imply that satisfies 28 -I Val < J. THM** +1 +• 3M(Kf(M4«))' + -^W+iv)] t * . (1.27) Now combining the estimates (1.19) and (1.27) we see that where >A is a constant dependent only on the data Taking O" to further satisfy A o-"» < \ , the conclusion of tire Theorem follows. Theorems 1.1 and 1.2 imply that for o->o(given in Theorem 1.2) (1.4), (1.5) has a unique solution for all t < <r in Co-,*"1 where M- a«"u,t|^ +1 . Note that o~ depends only on the data (1.29) . To complete the proof of uniqueness of the solution of (1.4), (1.5) we must show that any solution of (1.4), (1.5), irrespective of whether it belongs to - C,. M (where <r is the " <r" of Theorem 1.2), must coincide with the fixed point of H"1 in C <j-#!^ say -o , in their common interval of existence. (1.28) (1.29) (1.30) - 29 -If v(tV is another solution of (1.4), (1.5) on the interval C°,?l then we must show that v<-k)«v<i* on Co,5-] where J: •MW|»,»|) the common interval of existence. Note that when in Theorems 1.1 and 1.2: M is replaced by M*= -w\«-y ^ n-Gw«tfA^ we have that v,-C are both fixed points of T1 in C<J.',M' where in general a*<?. Hence we conclude that "v<^=v(t> on the interval Now if or, (<r, < 5) is such that -0<*l=O(i) on OJ O then it is clear from the integral equations (1.4), (1.5) that v <«•,)! ••*>< «-,) ,• Hence if \ V 15 <rj} S IT, , ^ ut <r,) , s < <M] are the temperature distribution and positions of the boundary at "tr <f, corresponding to v(U,tf(l) respectively then v< (5,0-,) s ^-C^.o-^ t >= si*-,}. Shifting the origin of time in Theorems 1.1 and 1.2 to we~'can again conclude that there exists an f>o such that vt^ivU^ on Co, 0-, + 4) , Since the only restriction on c( was that it satisfy c-, < 5- we conclude that •u<^ s vtt^ on L"«>»3 their common interval of existence. Existence and Uniqueness of the Solution for all 1 6 Q.T1 Let o-W) satisfy (1.6), (1.12), (1.22) and (1.30); .then there exists a unique solution of (1.4), (1.5) for t< Moving the origin of time to is ^ we can find a such that the - 30 -solution of (1.4), (1.5) exists and is unique for all t i o*1^ Continuing inductively we see that we can generate a sequence _ ** ^ o^*^._ such that (1.4), (1.5) has a unique solution for all 1 j. i *" . If we can show that there exists a 0"V such for each o-0) > 0-° (1.31) then we can conclude that for some A/ and hence (1.4), (1.5) has a unique solution for all t C- <°,T). However, this is immediate if we can find global upper bounds for For then o-° determined by the inequalities (1.6), (1.12), (1.22) and (1.30) with M replaced by and replaced by T* satisfies (1.31). Since "MKU,!) is continuous on ob we see that Lemma 1.3 is applicable and hence Therefore we have that (1.4), (1.5) and hence (0.3,a,b,c,d,e) - 31 -has a unique solution for all "t-* ( 0,TV- provided is bounded.on £O,T] and v0<x) is uniformly bounded on £o,bJ. In Chapter II we will outline the Similarity Method which will be used to derive the Similarity Solution (Chapter III) upon which is based the Similarity Algorithm (Chapter IV). CHAPTER II THE SIMILARITY METHOD The algorithm to be introduced in Chapter IV is based on particular solutions of the diffusion equation found by the Similarity Method. The following provides the theoretical basis as well as the procedure for constructing such solutions of differential equations. A common method of solving differential equations is to change variables in order to transform the equation to one whose solution is more easily obtainable. The transformations which give results are often those which exploit a symmetry of the original system. The Similarity Method provides a systematic recipe for finding such transformations using Lie (continuous) Groups. Sophus Lie showed that invariance of an ordinary differential equation under a one parameter continuous group of transformations leads directly to a reduction by one in the order of an ordinary differential equation. He showed how to find the "Lie" Group of transformations leaving invariant an ordinary - 32 -- 33 -differential equation^^ and found a subgroup of the full group of (2) the heat equation. However, it remained for authors of more recent years to show how to use continuous groups of transformations to reduce the number of variables, and hence find particular (3) solutions, of partial differential equations. The major contributions in this regard have come from Ovsjannikov £ 28 J, Matschat and Muller £2^]» and Bluman [ 2 J. More recently, Bluman £33»L"^3 has applied the Similarity Method to boundary value problems. (4) Lie Group of Transformations Central to the theory is the concept of a Lie Group of Transformations. Definition 2.1: -(-a Lie Group of Transformations). A one parameter family of transformations For a treatment of the Similarity Method applied to ordinary differential equations see Bluman and Cole Part I. (2) Lie did not see how to use invariance to construct particular solutions to partial differential equations. (3) For a thorough treatment of the Similarity Method as applicable to partial differential equations see Bluman and Cole [: 6 3 Part II. (4) Since we are interested only in a partial differential equation involving one dependent and two independent variables we restrict our attention to this case. - 34 -where and forms a Lie Group of Transformations with parameter £ if: (a) (Associative Property) there exists a function with for all a.bjCCQ such that, for x*x**jXeS satisfying (b) (Identity Element) there exists an €„£Q such that for all x € S ; (c) (Inverse Element) for every £ £ Q there exists an € Q such that We note that conditions (a), (b), (c) make the family a group of transformations, while the continuity conditions on •J'CXJO, make it a Lie Group of transformations. We - 35 -remark that by a suitable rep-arameterization, the identity element can be assumed to be zero. To apply the Similarity Method to a second order partial differential equation we consider the following Lie Group of transformations: X*= X<U,X,t;£) > (2.0) where -u. is the dependent variable and X, "t are the independent variables. Invariance. A partial differential equation G('U^>Mxt,Mtt,v<k,v<t,-w<>xl-fc)= o (2<1) together with the boundary conditions By(-u«,v4tlu,x,-b)so (2.1a) on the boundary curves Wvtx,t>= o Y= '.•••,p (2.1b) is invariant under (2*0) provided - 36 -(2.2) and (2.2a) on the boundary curves (2.2b) hold whenever (2.1,a,b) hold. That is, the governing differential equation, the boundary curves and the boundary conditions on these curves take the same form in both transformed and original variables. Since a partial differential equation seldom has a group rich enough to leave invariant boundary data such as. (2.1,a,b), we seek a group leaving invariant only the governing differential equation (2.1). What boundary conditions cannot be left invariant can frequently be satisfied by superposition (cf. Bluman and Cole [] 6 J Part II Chapter 11 ). In addition we can construct numerical solutions by "almost"satisfying certain boundary conditions'^ Further, useful particular solutions of (2.1) may be obtained by formulating boundary conditions in terms of the generated by the Similarity Algorithm is obtained by leaving invariant (0.3,a,b,c), by satisfying (0.3d) by superposition and by "almost" satisfying (0.3e). (5) The approximate solution of (0.3,a,b,c,d,e) invariants of the group leaving invariant (2.1). The Most General m-Parameter Lie Group of Transformations  Leaving Invariant (2.1) First we note that the transformations (2.0) on the variables induce transformations on the derivatives, which together with (2.0) constitute what are referred to as the Extended Transformations. These also form a Lie Group of transformations. Before proceeding further it is necessary to reformulate invariance in a more useful way. To this end, we introduce the infinitesimal transformations. Noting that TJ'CH^t-O, X(-*,wt\,i) , T (M,x,t;<r) 6 C °°(//?3,) we expand about f = o (the identity) to obtain (2.0) in infinitesimal form, where 38 »i^v,^«l>T(u)»,tiOl The transformations (2.3) induce transformations on the derivatives, i.e., and where Similarly we obtain where To obtain the second extensions we write - 40 -i.e. 7x* £5* \S*4v< 8K'/M< C>*» * Ve>~1 ,3*0.*/ * * dv. 3R <JW V Now (2.2) can be expanded about fro to yield - 41 -G IK-*-. K- e» *Vt»» • **V»**>**) = ° + f GC^xJ«,UMtlv»iB-tlv«,Mt,wyicl-t) . +• ot£») where we have introduced the first order differential operator It can be seen that invariance of (2.1) under (2.0) is equivalent to XG = 0 (2.4) whenever G = O , With this formulation of invariance we are prepared to find the most general m-parameter group leaving (2.1) invariant Substituting f,?>r, *?t , 9?xt ,r?it into (2.4) and using the relation G~0, we obtain the determining equations for by setting equal to zero the coefficients of the independent derivative terms (^KI )• We are left with a set of linear partial differential equations for - 42 -Reducing the Number of Variables To every Lie Group of transformations there corresponds a set of Canonical Coordinates, in which the group is a translation of one of the variables. Using these Canonical Coordinates it can be shown that if the translated variable is an independent variable then invariance of a partial differential equation under a one parameter Lie Group of transformations leads to a reduction by one in the number of independent variables provided the solution is unique (cf. Bluman and Cole [^6j Part II Chapter 3). number of independent variables by one leaves us with an ordinary differential equation. Suppose (2.1,a,b) is invariant under (2.0), whose infinitesimal transformations are given by (2.3). If w* toe*,*) is a solution of (2.1) then both v«s &ix*,t') and vv."*U(©.**'tiO are solutions of (2.2). Now assuming (2.2) has a unique solution then v = v^*' . Expanding V, w.v about £zo and gathering-terms in powers of <f we find that It should be noted that in this instance reducing the (2.5) (The Invariant Surface Condition) must be satisfied if viv* and - 43 -conversely. The general solution of (2.5) can be found by solving the characteristic equations c\x s . - d ® > (2.6) • g<G.«,±Y • (6) If %/t is independent of © , then the general solution of (2.6) takes the form © = €(*,±> \, where fs$(*|t) (the Similarity Variable) and are the two constants generated by-solving (2.6). Substituting © into (2.1) and using the relation $ V*,*) we obtain an ordinary differential equation for ^C^) Hence the number of variables has been reduced by one. The complete solution of (2.1) can be found by solving the ordinary differential equation for £t0. However, if a two parameter Lie Group of transformations leaves (2.1,a,b) invariant and the invariants (the similarity variables) (8) associated with the two parameters are functionally independent If 2*/T depends on © then the general solution of (2.6) is of the form O = «&<€>.*, t;i,V CO) and $-^(G>,.*,tK The boundary conditions (2.1 a,b) become boundary, conditions to be satisfied by ?^V>-(Q\ Two invariants are functionally independent provided their respective infinitesimal operators are linearly independent over the field of complex functions (cf. Bluman and Cole [6^ Part II Chapter 8). - 44 -then the solution of (2.1) can be found directly using the invariants without recourse to (2.1) (cf. Bluman and Cole [ 6 ] Part II Chapter 8). In general, if an m-parameter Lie Group of transformations leaves invariant a partial differential equation with accompanying boundary conditions, it is necessary that the associated invariants (similarity variables) be functionally independent before we are assured that the number of variables can be reduced by m. (9) The Classical Group of the Heat Equation Considering the invariance condition (2.4) implies whose solution yields the six parameter group: 7cwue[.*\j\|^|x*i] • ACt> j (2.7) (9) See Bluman and Cole \j>~\-Here «<, v^'V, 6, K, \ are the parameters of th.e group while • is an arbitrary solution of the heat equation. The group (2.7) in the (x,t) plane is a subgroup of (10) the eight parameter projective group. The parameters <*, K represent translations in the i.. and x directions respectively; "O represents a stretching or similitudinuous transformation; while ^ is associated with the Galilean transformation. To find the form of the global transformation associated with V we solve, the set of characteristic equations The resulting transformations are given by In the next chapter a subgroup of (2.7) will be used to construct the similarity solution central to the Similarity Algorithm. See Bluman and Cole [^6] Part I Chapter 7. CHAPTER III A USEFUL SIMILARITY SOLUTION OF THE HEAT EQUATION FOR THE STEFAN PROBLEM In this chapter we will use the ..Similarity Method, as does Bluman [ 4 ], to derive the solution to an Inverse Stefan Problem corresponding to the boundary melting at a constant velocity. We proceed as follows. Given set) , the system (0.3,a,b,c,d) reduces to the Inverse Stefan Problem: (3.0) i^Co.T)-, (3.0a) u(ito,t)-o 't«Co/r>, (3.0b) MK<»,±)^ o (1) tCCo.r), (3.0c) *A(*,O\- Vi„(V> xeC«iCl. (3.0d) ^ The methods of this chapter may be used to deal with the boundary conditions vWo,t>- PCt) or uK(o,t>-RC-O • - 46 -- 47 -We will show using the Similarity Method that for a member of.a two parameter family of curves, the system (3.0,a,b,c,d) has a closed form analytic solution. For convenience the group (2.5) together with a first extension Is given below, «t + 2M t + Vtl %<*>±\-- K+• St f»x + Yxt (3-1) where FCX,+,)=-Y[|\ J] - |* +- A dx V ^x/ -ax' * a* We consider (3.1) with ^cx,t)5o and note that if the boundaries x=o, x.-sit) are invariant under (3.1) and l"^6'^1 °, then C3.0,a,b,c) is left invariant by (3-1). The condition if''^"»o is satisfied provided &=o ; x-o is invariant under (3.1) if and only if x*"= o whenever XJO , i.e., S*K=o . The invariance ofx=^Ct) under (3.1) implies that satisfy the differential equation - 48 -^"(.SCt),*) «. irt) 3T(k) . (3.2) •When combined with .S<<rt = c (3.2) implies that Hence the three parameter subgroup of (3-1) leaving invariant (3.0,a,b,c) whenever Jt-kJ«- (c'uvt+yt1)^ (3-3) is given by: ?(<,t)t x(vf Vt) f (3.4) Using (3.4) the Similarity Solution of the system ( 3. 0, a ,b , c , d) corresponding to the most general "moving" boundary (3.3) can be constructed (cf. Bluman [ 4 ] ). However, for our purposes we only consider a subgroup of (3.3) to obtain the Similarity Solution of the system (3.0,a,b,c,d) corresponding to the "moving" boundary - 49 -Letting v--/3c > V-^i* (3.4) reduces to + +>^| ^ (3-5) where The infinitesimals (3.5) yield the set of characteristic equations From the first equality of (3.6) we obtain the Similarity Variable V= x/(cy3i) » where K = o <=> $ = O, Integrating the second equality of (3.6) along the similarity curves ^ = constant we obtain the solution surface (3.7) of (3.6). Here ^((ijjx) must satisfy a certain ordinary differential equation together with the boundary conditions -^--•jx)-- X^Oj/*) - o for all ^ (3.8) if v-CXjij^*) is to be a solution of (3.0,a,b,c) with eyai', To derive the differential equation satisfied by we write vet*,!) (the solution of (3.0,a,b ,c ,d) with iCt)-c-/3± ) as a superposition of functions, £i j^O, and substitute the resulting expression into (3.0). Introducing the variables into (3.7) we obtain v^tj^-.vU.t^)-. ^JF^pe *tr+/U>c/(*,P> <P*. • (3.9) The form of (3.9) suggests we take - 51 -(3.10) Substituting (3.10) into (3.0) and taking into account the boundary conditions (3.8), ve see that c^f(^;j>) must satisfy ^ -As*1 The solution of (3.11) is where we have introduced $C^)- Sc e ^ v. <«^> (3.11) (3.12) Substituting (3.12) into (3.10) and evaluating the Laplace - 52 -Transform by closing the contour in the left half plane we obtain /3-aL <*=> - u£t . (2) U(*,0»J^erteytiJ^T.e «7yiU co«(wH.a_) (/S) (3.13) where © If we evaluate the Laplace Transform by expanding for large p we obtain the small time representation of (3.13) UCx,-t>-- ^ ^(SpGCx/t^) dS-, (3.14). where n r x* fa G(M-^)--T=Ue- ' JE(-») V + <? (2) _ The expression %*. C*, ) of (3.7) can be substituted directly into (3.0). The result when combined with the boundary conditions (3.8) is a Regular Sturm-Liouville System for the eigenfunctions ^-?((il**M}. > '•N* st*-i[hr . The solution (3.13) is then obtained as a linear combination of the eigenfunctions - 53 -We notice that when y3-ot ( 3. 0 ,a,b ,c,d) reduces to the usual fixed boundary problem and (3.13) reduces to the well-known Fourier Series Solution. Clearly the series (3.13) together with all its derivatives converges uniformly on t^Cy^t] for any fixed "fc>0. If "U^txi € c'CojCl and "uc (O =, iue lc) ~° then we will show that (3.13) is uniformly convergent on at -t - c , and that v„U,i:o) converges to' vve CO . Since the functions jc»s(i«3„)f){ satisfy a regular Sturm-Liouville Problem on C©. »"V > for any function ^<JO C C* C°I»"} we can write T"(x)= ZZZzos^x)} fct*0 C«J(W^,1^ , (3.15) If the sum in (3.15) is uniformly convergent we also have V'OO (WAK) ^ Vc^i cosK^ «J.»j . (3.16) We claim the following. Lemma 3.1. Let ^Cx) 6 C* Co,."] with £cn=V'<oV-o then the sums in (3.15) and (3.1-6) converge uniformly on C°, il to tCx) and V(x) respectively. - 54 -Proof. Let J (x) be the following extension to C--?)2! of (l) 5co* *M * s C°-.»l, (2) §^ = $(->o x eC-a,*7, (3) i(«+i)s-i <x-i) xe C*l,.ll. ' Clearly ^(0 is piecewise C* on we claim A € C*C-2, il , Only the points x-o, 1 need be checked to ensure this result. X=o: £(0-)- cfv^. I <-> o L © - C j « ?(o+) =. - +"'(0+) = O ; c->o L «r J CfM)-l J c ->•© L Since $"00 is piecewise C* as well as C we have that -55 -00 «l 4 where the sura is uniformly convergent on C©,il (cf. Courant and Hilbert [8] Chapter 2) . That the sum in (3.15) converges uniformly on C°> il follows trivially since § <*> G C* C-*,*J (cf. Courant and Hilbert ^8] Chapter 2), and hence the lemma is proven. Applying Lemma 3.1 to vt6(.e*> we see that - * €Vc ZT t^K^) ) . « ' cci(w,«l^ (3.17) and e - 2^ S^Kf) \ e %,.(^ -,Cw^Ua (3.18) where both sums converge uniformly on C°itl. It should be noted that the results (3.17) and (3.18) are essential since the Fourier Series Expansions for v*,U,t> - 56 -and U„ Cx*i(t),-t) are used near i-o in the Similarity Algorithm. Further note that the hypotheses of Lemma 3.1 are satisfied at each step of the Similarity Algorithm after the initial step (see Chapter IV). The initial condition "HeCx) need not satisfy the hypotheses of Lemma 3.1, that is, n0 <x) may not satisfy (0.3b,c). If V„(x5o)^o the Similarity Algorithm is not affected, that is, the sum in (3.17) still converges uniformly . on [ o,C ], while the sum in (3.18) only converges uniformly on [S,c ] for any b>o . However if v,tt^o (the slab is hot prepared to melt) then we must modify the algorithm and use the usual Fourier Series solution of (3.0,a,c,d) and (we refer to this as the fixed boundary solution) with \$(i}=o until \*.(c)t) = o (the standard derivation of this solution is given in Appendix E). We see that there is some difficulty with this algorithm if at any time the heat flux is insufficient to maintain melting. Since we provide no mechanism for freezing, we must, at each time step, determine if the heat flux is sufficient to maintain melting, - 57 -i.e., determine the sign of S(t;),t;) - W\ (!«•)), *;6^~.(cf. Chapter IV) If A>o we utilize the fixed boundary solution until the heat flux is again sufficient to maintain melting. In what follows we will assume the heat flux to always be sufficient to maintain melting. Returning to the group (3.1) we note that setting Y5 ° aid proceeding as above we generate the particular solutions of (3.0,a,b,c,d) given by Sanders £3l] . Because of their relative simplicity, the Trigonometric functions lend themselves much more easily to numerical calculation than "do the Confluent Hypergeometric functions which are the basis of Sanders' solutions. Apart from the inherent difficulties in calculating with the Confluent Hypergeometric functions, there are also convergence questions which would place in doubt the utility of such a scheme. For much the same reasons, an algorithm, similar to that given in Chapter IV, based on the most general "moving boundary" would encounter difficulties from the onset, as here too the solution of (3.0,a,b ,c,d) is expressed as a sum of Confluent Hypergeometric functions. - 58 -/ O \ We remark that in 1939 Huber [20^ proposed essentially the algorithm of Chapter IV to approximate solutions of certain one-dimensional two-phase Stefan Problems. Huber's solution, however, is not based on a similarity solution. He eliminates the intial condition by introducing the usual source solution of the heat equation (J*C (*,i:,i;,r.) of .Ch.ap.ter ..I.)., then uses a set of Appell Transformations to transform X-JC-t)=. c.-y3t to the fixed boundary y- X , while leaving invariant the heat equation. The solution is then given as a sum of a source term plus a complicated convolution integral. Huber's solution is. too unwieldy for numerical purposes. Recently Rubinstein [30j has suggested that Huber's method can be significantly simplified by using a Green's function on the domain x e (0, cys-fc)^ t <c (o.,i)j first derived by Soloviev [33]. By so doing the need for the complicated Appell transformations is eliminated and the solution can be given directly ( 3 ) Recently A. Fasano and M. Primcerio [l4] have demonstrated convergence of Huber's Method for a one-dimensional single phase Stefan Problem. (4) The representation suggested by Rubinstein is *u<*,tU %cC$)r/rrx,t-,5,.w •*•<*,t',§>}«l$ where .5'(*,f,^ ) • is the Green's function of Soloviev, and v,,^) is the initial condition. We remark that this is basically the representation given by (3.14). However, the Green's function GCx.-t-,^) in (3.14) is given in a more compact form than is the Green's function, K(x,t'>5,o) + ^(x,. - 59 -The calculations involved in using Ruber's Method, with Rubinstein's simplification are very similar to those involved in the Similarity Algorithm if we were to use the small time representation (3.14) of the Similarity solution. We will argue in Chapter VI that the large time representation (3.13) is more practical than its small time counterpart, (3.14), and hence that the Similarity Algorithm is more useful than Huber's Algorithm. ! CHAPTER IV THE SIMILARITY ALGORITHM Having derived the Similarity Solution (3.13) we are now ready to outline the Similarity Algorithm. Suppose that we are given the pair of functions C"VA,S) satisfying the system (0.3,a,b,c,d,e) and we wish to find an approximation (£,$) to (vision I[O,T]. We proceed by partitioning the time interval (not necessarily a uniform partition) and estimate >5("t) on £to,\, by Sit)'- Ce - /3„ ( t-to) , t0<t<t, where As we have seen (cf. Chapter III), for 5Ct) so defined Similarity Methods yield a closed form solution of (0.3,a,b,c,d). We denote - 60 -- 61 -this solution by V'tK, i't,) and remark that it is valid on the domain To extend our estimate to (^»;^al we define |<t>* <t-f >, t, < t 5-t4 where and /3, is calculated by substituting y.*(e,;-ai,) into (0.3e) to obtain /3,« *'< Wet,) - U„*t,)) . Now considering as the origin of time in (0 .3 , a ,b ,c ,d) and Vt°Cx, ^i|) as the initial condition UetO we can, as before generate a solution M.' (*) •i.-t,) of (0.3 ,a,b ,c ,d) valid on the domain Continuing inductively we obtain the approximate solution on the interval ( t^ii+j^ by defining where - 62 -We obtain "U* Ix, i-t^) by taking to be the origin of time in (0.3,a,b,c ,d) and the initial condition Uc Lx) to be w**'( x> • Again the Similarity Method yields "u^x, t-t{) on Fig. 4.0 • - The Similarity Algorithm Time - t The approximate solution C t S) is then taken to be We prove in Chapter V that as we refine the partition in such a manner that -w\ <x^ J4,t^ —> O the approximation C"M.j5 ) converges in the supremum norm to (\*ts) - 63 -That is, given £>o there exists a S>e such that w\ew» < & implies t and . A (1) ? In passing we remark that here we have considered the case where the flux condition (0.3e) is satisfied at the endpoints t-i of the subintervals C "t«^ . I*1 fact, we could have satisfied (0.3e) anywhere in the subinterval. However we see no particular advantage in doing so, while, as will become apparent, the choice of a point inside the subinterval complicates the actual numerical procedure. A. Note that a heat flux V* Ct) is induced by the approximation {y^tt) i.e., one can calculate the heat flux which produces the melting described by This heat flux satisfies for all "kj € 7T , and as will be shown in Chapter VI, can be used as a indicator of the errors ^ Notation: given 4~Cx,t,> > <^t*(t) and e!tt)^0 then «»i-.«-jt--.«iUWM:f4«ki\,*i*'*,-'S,KA,'lJ-- 64 -and ( SCt), 3 (tl) i.e., the quantity \ UU>-KCT>| 4* (4.0) can be computed and used to determine whether 77* should be refined further. Although in principle (4.0) can be computed, it requires a good deal of labour, and instead we define ECi) = «<* \ (van- wunar c and calculate C1 /Ci> 1 0 ^ o o J (a relatively inexpensive calculation). We see that £Ct) gives an indication of how closely (0.3f) is satisfied by CviJ ) Moreover, if £(t) is small for all Co,T] we expect (vlji) to be a good approximation to Cvt.s) . We remark that tu,s) provides an analytic approximation to the solution < v\,s) °f (0.3 ,a ,b ,c ,d ,e) . That is j satisfies (0.3,a,b,c,d) exactly while it "almost" satisfies (0.3e). We will use this property in Chapter V to demonstrate the convergence of the Similarity Algorithm. CHAPTER V THE CONVERGENCE OF THE SIMILARITY ALGORITHM Since the pair of functions generated by the Similarity Algorithm is an exact solution of (0 .3 ,a,b ,c ,d, e) A with V^tt) replaced by the interpolate, V\(t) , it follows that the convergence of the Similarity Algorithm is equivalent to the continuous dependence of the solution of (0.3,a,b,c,d,e) (1) A on the boundary data V\,(i) , if we can show that Vti) tends uniformly to V\(t) on £0;T*3 as we refine the partition (cf. Chapter IV) such that v^»^ 1 At;|->o, Continuous Dependence for Small Time r We proceed by demonstrating that the system of equations (0.3,a,b,c,d ,e) is continuously dependent on the boundary data j b, Vlcl.x), Mtlj for a small time o->o. More precisely, given CM,S) and satisfying (0.3 , a,b , c ,d, e) with H(t) Co, b,] j (5.0) K,U)! THx"* on In fact we will show that <V,A> depends continuously on the boundary data | b, "ucu>, W t-uj . - 65 -and W(-M- R(t) (5.1) V.OO* on Co,bal (2) respectively, where W, > bj and j (T> - b, -r(T), we have the following Theorem. Theorem 5.1. If U,s), satisfy (0 .3 ,a ,b ,c ,d ,e) with (5.0), (5.1) respectively, then there exists a c->o such that the following inequalities hold: r^-sii^j: d>)v>;-fcai + o-Vuv- 1^ A3 c-|lH-R«t (5.2) where cl -r^l^. ^ Tli) , s tt^ and the constants depend only on data £ T; H<«,R't^u^^^oo, b^b^dfo),^ • (5.4) That is, the system of equations (0.3,a,b,c,d,e) is continuously (2) The functions H(t^, R(t^ are taken large enough to sustain melting. - 67 -dependent on the boundary data |b)M0Cx^ h (Uj . In other words, the,system (0.3,a,b,c,d ,e) is ,s.t able .with respect to variations in the boundary data, . Proof. With the definitions we note that the following equations hold (cf. equations (1.4), (1.5)) (5.5) X (5.5a) and o (5.6) c .t (5.6a) o - 68 -In addition, subtracting (5.6a) from (5.5a) we find that « lb,-b4J + v->i¥--t + K'-rw-'RU^t . (5.7) To obtain the inequality (5.2) we first subtract (5.6) from (5.5).and write the resulting expression as where and Since G^-Ot,t; r\-- <5*<*,t; £ and NKb^- 6"<*,t;oj0^ O, the expression for I can be rewritten in the form where r * © 69 o It is easy to see that V/f satisfies w,\ * "y-^'ib, . (5.8) To estimate Vj we write where o Proceeding, as in Theorem 1.2 with v^' we obtain |v;i v< A «i\ ,. (5.9) while applying the Mean Value Theorem to v^" and evaluating the resulting integral we can show that t V4" | N< JL liil^ Hs-*H, . (5.10) - 70 -From (5.9) and (5.10) we conclude that satisfies I V.I $ 1 lla- . (5.11) Finally we estimate by writing where V^, -z) *>(?>[ <ii!!if> >] AS,. Substituting v* (Sl*)-%)/ztv* into Vj and using the inequality Iwle $1 we obtain \VJl * ftJPflv,, (5.12) T*» t. Similarly we find that V" satisfies Tt"t t The estimates (5.12) and (5.13) imply that satisfies - 71 -7T '* t Combining (5.8), (5.11) and (5.14) we find that \1\< ±. W"LH*->H. where We now estimate 7" by writing wt * 2 J ) | vet) isa^sm KU<t),t-jStr>tti 3(t-t) J where 72 The term WJ can be expressed as the sum W, = VA/,' + Wj" + W(MI 2 V -tSiS (»(t)»»(«> J ,< (.set,, t • jew, r) Since .5 <t.V is Lipschitz continuous we see that iw/ix< »HJ!T II^.^II t"». (5.i6) Then noting that and - 73 -we find that \ Vct1-5<f)X _ .(t(t>-r(t)H < otitfty-^H. + itH-nirl (5 17) and hence that w," satisfies Iw,"!* Q('u^*»i; [l\v-xnt+ l>H-^nT] t"*. (5.18) Finally to estimate VA/,'" we first apply the Mean Value Theorem for a function of two variables (.x,^) to K (*, t; ^f, t) and use the fact that V-CtV is Lipschitz continuous. To the resulting expression — * J Ct-*> U-ti I (where X is between -cctN and 4<t^ and Cj is between +U) and Hz) ) we apply the inequalities |v| e * I and (5.17) to obtain lw,'"| 5 ll/^i^nftll^. ^Ilv-,*.!^ 4 |IH-RIT|t. (5.19) - 74 -Together (5.16), (5.18) and (5.19) imply that Ti"t ( Now-we write VA/^ as where (5.20) - K(f<t\,t;-*<*),7U dr Since b, >*(t^ s(Vt > b* we obtain (5.21) and - 75 -|W;w ajMI4U3-**. \ 1 ±*> A7 (5.22) To estimate w^'" we apply the Mean Value Theorem for a function two variables (*,£) to K(*,t;-f, <r) and obtain .A. A (where if is between 5 ft* and r(t> and «^ is between sm and •+(*) ) . Since * b« $ 2+ a t>, we have tw^'l* CV^'** dr. (5.23) o Combining (5.21), (5.22) and (5.23) we see that WA satisfies **** 1 (t-t>*4 J Making the substituting v = W/(t-*)v* it is easY to see that - 76 -and hence using the inequalities «t^e _L J. and ™ 3 <2 J ~ we obtain r (5*24) fit — <''4><V)! Thus -X satisfies Combined with (5.15), (5.25) implies k ill (5.25) (5.26) (5.27) - 77 -where C, is a constant dependent only on the data (5.4). Now using (5.26) and (5.7) we obtain n v s< Cz ^ 1 I «,,. ^ , + ^ II V- v«bj + t'SlVviAi^ + tv* IIH-RHT^t where is a constant dependent only on the data (5.4). Now take «j?c(<v) and let c-> o satisfy CJ»-,/*s I-y . For i«(o,») we see from (5.27) that and hence using (5.7) we have fc*-*^ j A, lV>,-b2l +^,*,'4|^-vf|l|,1 1 Aj»'l,H'fillT. That is, (5.2) holds. To obtain (5.3) we note that eu,U"> w(«,i) satisfies the heat equation in the region j C'f.tV: o<-t<-r| . Hence by the Maximum Principle (Proposition 1.1) we see that the - 78 -maximum and minimum values of 9Cx,i.) occur at o , x - o or x - d(t) . There are three possible cases. Case I: )e(*,i)| attains its maximum value at 1=0, here II6C -,a-> II < || Vii. , (5.29) Case II: le(x, i>| takes on it maximum value at * d Ci) (say for fixed "fc eUt)- •*•(+.) ), here where ^€ Cx, + tt^ and hence xe(©,5(t.) t€<o,T) Using Lemma 1.3 we obtain net',*">tld(0 5 ™**{BH!V,M*rtk,| lU'^V • (5.30) Case III: . takes on its maximum value at -X=o , here we use the integral equation (1.2) to obtain where - 79 -o It is easy to see that and hence that - 80 -Using (5.29), (5.30), (5.31), (5.28) and (5.2) it is easy to see that Continuous Dependence for All Time ^ £°> T3 If in Theorem 5.1 we replace ftVllu by rtv(-,tH!3(t, ft^rl^ by u W(., t)H . ... <*>H. by ft-u, (-.t^,,., H VII^ by « ^.(-jtiH^, then we can find a new t\ say and hence a new C" say CT such that (5.2) and (5.3) hold at any time y i, C C°. T- <rc] with S>(^)= w(^,) , - 81 -That is |j(t)-v<t)| (5.32) V + o-„ HH-ftllT } d(i1+ cr.) < G, ISCV>-V<i,)\ +• T^llViC-.t.l-vA/te.lllj^ (5.33) for new constants ^ At , A^ , Aj , Q, , B, j 84^ . Now (5.32) and (5.33) imply that (0. 3 , a ,b ,c , d , e) is continuously dependent on the boundary data | b, "M6<*I( V>(t)| for all time *t € CoiTJ^ since we can take A> such that '//j < ^ < y^N^ and apply (5.32) and (5.33) successively at the times ^ i* o, <r0, a <r0l.. •, Vi(.*/-/)J to obtain M-+»IT \ / Ai'- A, <r/» A3o; \ / li>,-ij \ IIH- Rll, - 82 -Hence (0.3 ,a,b ,c ,d ,e) is continuously dependent on the initial data \b, M0oo,Mtt|for all i» <c O ,T3 . We remark that for numerical considerations the size of the constants Qt } i* •» *j 3is of considerable importance. That they appear to be large, we feel is a failing of the. method of proof and not characteristic.of the actual system. Convergence of the Similarity Algorithm If (vi,4) (cf. Chapter IV) is the solution of (0.3,a,b,c,d,e) generated by the Similarity Algorithm with V\tt> (the induced heat flux) and (vi.,s) is the exact solution of (0.3,a,b,c ,d,e) we see that fts-slb > \UV-l-^U- u(•,1 • t-« E»'.*rJ ( •*wu^l.jJ4>i (where T=-\"TJT^ and T satisfies J(.T)-b^ ) depend <v A only on 1\*A»H\\ — # Hence we must show that V\ tends to as vv^ft^t <^-t^ _>e (cf. Chapter IV) in order to prove that the Similarity Algorithm converges. Suppose "tj is a point of the partition 7T then we will show To accomplish this we return to the expression (3.14) for the Similarity Solution. Since we want *u.*K ( S (i^t^ i) for small "t we differentiate (3.14) and evaluate "M. * ( s (-v -i), t ) - 83 -asymptotically for small "t (the detailed calculation is given in Appendix F), To first order we obtain "II 4 If K tx) is continuous from the right at all points of the partition 7t we have hCti+t>*htt^ 4 o(i) as t->o and hence for i <f ( ©, ,) ©<* ( Klii + t) * h + Thus V\(t} tends to on £O^T] as max ^1- —> o and hence we have shown that the Similarity Algorithm coverges in the sense of Chapter IV. Order of Convergence Having shown that the Similarity Algorithm converges we turn our attention to the rate at which it converges, i.e. its order of convergence. - 84 -If (u,s), (u,S) are as above, then we say that Cltx, tW> and j ttl -> Sit) with order of convergence and p: respectively provided and respectively as (<vtt) —> d . To establish the order of convergence of the Similarity Algorithm we assume that W (t) satisfies WUi+t)* K(ti)-+O.C.t.) (5.34) for all points "fc.^ of the partition If , then we have that «W-UlT < 1 ^^^M^Cc;j<iti)--/J.<'(c^t<.U(^t; + ()'/i + 0(*t Since V.Cx.'t) satisfies the heat equation at x=.S(i) the expression tL P vi. ( s tiXtll - o implies that ett J and hence that - 85 -Thus we have Hence if V\ Ct) satisfies (5.34), we have that the order of convergence of the Similarity Algorithm is one half. Since in practise is small we assert that the effective order of convergence of the Similarity Algorithm is between one half and one. CHAPTER VI THE SIMILARITY ALGORITHM NUMERICAL RESULTS In this chapter the results of our numerical experiments with the Similarity Algorithm are given. We present several examples illustrating the properties of the algorithm, including its order of convergence, and suggest two ways of increasing its accuracy. In addition, we attempt to justify the use of the large time representation (3.1.3) rather than the small time representation (3.14) in the Similarity Algorithm. We conclude the chapter by comparing the Similarity Algorithm with Lotkin's Difference Scheme. Numerical Examples By presenting the following numerical examples, we attempt to bring to light the advantages as well as the dis advantages of using the Similarity Algorithm. We first consider (0.3,a,b,c,d,e) with - 86 -- 87 -3 / (6.0) For the data (6.0) Sanders [ 31 ] has given the exact solution m«(t)' JtJ-"(i-a) on o < * •< .sc-n * 0 - * t)v*. Comparisons of the exact solution with the approximating solutions are- summarized by Table 6.0 and Figures 6.0, 6.1, 6.2. Here six terms of the Similarity Solution (3.13) and equal time increments ( AI^ - at for all i ) are used. In each case the approximation is used to 80% of the total melting time, i.e. T=.4. In what follows we use the notation e^(Ty = IU-J«T, If more than two or three terms of the series (3.13) are used, we have found Filon's Rule for integrating S* * coi eU ( K a real number) (cf. Filon £15 ], Davis and Rabinowitz [] 9 ] page 62) to be the most efficient method of generating the coefficients ^n(/3) (cf. Chapter III). - 88 -where ^ u. > £ C*,*), J C-O, .r Ci)^ £ Ct) j are defined as in Chapter IV. Table 6.0 Errors Versus Time Increment (Boundary Data (6.0)) % % •<wt gw( .4) Error (.4) Error ew (.4) .200 .87(-l)(2) 11 .89(-l) 20 .96(^1) .05 .25(-l) 3 .27(-l) 6 .30(-l) .01 .8 (-2) 1 .8 (-2) 2 .9 (-2) .001 .6 (-2) .75 .3 (-2) .6 .4 (-2) From Table 6.0 and our 'numerical examples it seems that an accuracy of one to five percent .is easily obtained.. However, higher accuracy is difficult to achieve. For instance, we see that with -6.t =.01 the algorithm requires 40 time steps and leads to errors ^(.4), 6S(.4) which.are smaller than 1% and 27„ respectively. However, for accuracy better than ?^(.4) = .006 (,757o) and es(.4> = .003 (.6%) more than 400 time steps are necessary. The reliability of the last column of Table 6.0, ©^CT) t as an indicator of the errors CT) and 6a CT) is difficult to assess. However, the solution generated by the (2) v\ Here we introduce the notation ek.C>M= a. x \o a € (0,1), vi an integer, - 89 -Similarity Algorithm satisfies (0.3,a,b,c ,d) exactly. Hence, for the Similarity Algorithm (,TV-> O is a necessary and sufficient condition for the convergence of the algorithm. In addition, a rough calculation shows that if i ' ' > i ' then vv% »-y /i I- —> <j 0 ' A where V\ CIO is the given heat flux and \-»(t) is the heat flux generated by the algorithm. Hence we consider ^j) to be a "rough" indicator of the errors S^CT) a«d es (.T) • Moreover, ve take the order of convergence of Q. (y)-> o * A to be an estimate of the orders of convergence of ©^CTi-^o and C -> O . We remark that for any scheme leading to an approximate solution of (0.3,a,b,c,d,e), the quantity &^ (T) can be calculated. However, in these cases ( T) -> u would not in general be equivalent to convergence of the corresponding scheme. For instance, in the cases of finite difference schemes and the Collocation scheme (cf. Chapter VII) C?w KT) is only an indicator of the truncation errors of the schemes. - 90 -R U(X.T) VS X TIME=G.40 - 91 -Fig. 6.1 Approximations to the Position of the Boundary Stt) up to T=.4 for the  Boundary Data (6.0) Fig. 6.2 - 92 -Comparing ^> (V> and Mi.) for the Boundary Data (6.0) HEAT FLUX VS TIME -i— -r-II - 93 -Figures 6.0., .6.1 and 6.2' show that the accuracy of the Similarity Algorithm depends on how closely the generated heat flux, W("t) , approximates the given heat flux, V% tt} . This illustrates the proof of convergence. For the Similarity Algorithm to be practical we should have to use at most six .to eight" terms ,o,f~ the (3) (3.13) during most of the calculation. Experimentally, fo^ smooth initial temperature distributions, the one given in (6.0), we find that six to eight terms is more than adequate for results similar to those given in Table 6.0. However, more terms- of the series in (3.13) are necessary when the initial temperature distribution is rich in the higher frequencies. The number of required terms is governed by how closely (3.13) evaluated at "t50 reproduces the initial condition. Although initially a relatively large number of terms may be required, the following example (see Fig. 6.3) shows that during a relatively short initial period of time (short compared to the total melting time) the higher frequencies are largely attenuated. This is a consequence of the dissipative character of the.heat equation. Hence only the first few terms (3) A bound on the error made in truncating the series in (3.13) will be given later in this chapter. - 94 -of (3.13) need be retained for most of the calculation. As an example, we consider (0.3,a,b,c,d,e) with ~u0u\ - u~>) e (t) - '°/3 J> = 4 and use the Similarity Algorithm with equal time increments •oA =.001 and fifteen terms of the series in (3.13) to obtain an approximate solution. Fig. 6.3 The Approximate Temperature Distr.ibu'tion for t. B'etw.een 0 .and .1 (6.1) for the Boundary Data (6.1) - 95 -Figure 6.3 shows the approximate temperature distribution at t =0, .001, .01 for the boundary data (6.1). It can be seen that at t. =.01 (about five percent of the total melting time) the high frequency components of the initial temperature distribution have been significantly damped. Hence by this time fewer than fifteen terms of (3.13) are needed in the calculation. Furthermore, our numerical experiments indicate that the initial temperature distribution is of importance only initially during the calculation. To a large extent the long term behaviour of the solution of (0.3,a,b,c,d,e) seems to be independent of the shape of the initial temperature distribution. Optimization of the Similarity Algorithm So far we have made no attempt to optimize the accuracy of the algorithm. This can be accomplished by choosing an optimal partition 7T and, or an optimal value at each time ~t± of IT . It is clear that any strategy aimed at optimizing these choices should be guided by a desire to have V\ (t 1 approximate ta<i} closely (the proof of convergence) or at least that j Kir),!? approximate ) MtWr well for o o all i. . Below we introduce two modifications to the Similarity - 96 -Algorithm as a step towards optimization. We first note that the Similarity Algorithm provides an exact solution when the boundary moves at a constant speed. Hence when the boundary moves at a slowly varying speed, i.e. |i*Ct)| small, the straight line approximation to the boundary should lead to better results than when \^Ctl\ is large. Thus we concentrate points of the partition II during periods of time when \ J C--k\\ is large, which, for the most part, corresponds to periods of time when V»(t) ±s changing most rapidly. As an example, for the data (6.0), the points 4.^ of 77* can be taken to satisfy I: j U(r).lT= C • ti for an appropriate constant C>0 . The next modification is motivated by the proof of convergence and is aimed at optimizing the choice of on the time interval £t^v"^<+\) f°r a given partition 7T The strategy is to add in the time interval r^o^< + i) a portion of the heat which was missing in the previous time interval E^N'-I, "tj ) • More precisely, we choose /3(-on t *i , tj-M^ to satisfy II: /3.-M«| MV) + *-(W(iiO-£<t;-))- <"fC;,*t;)j - 97 -for some >,«? (see Fig. 6.4). Fig. 6.4 Comparing the Given Heat Flux with that Generated by the Similarity  Algorithm Using Modification II Time - t We remark that both modifications I and II can be implemented at very little computational 'expense. To illustrate the utility of the above modifications, we consider (0.3,a,b ,c,d,e) with the data (6.0). With C= )K(Mc!r and V, - \ for all i we obtain the following results. Table 6.1 Error Versus Time Increment -<=>*t Using Modifications I and II (a) fit ( .4) • - Error in JCt) Unaltered Algorithm Modifications i II. I and II .05 . -27(-l) 6% T22(-l) 57. .13C-1) 37. .99(-2) 27. .01 •75(-2) 1.77. •65(-2) 1.5% •40(-2) .9% •36(-2) .8% .005 •58(-2) 1.3% .44<-2) 1.07. • .31(-2) .7% .29(-2) .77. - 98 -(b) e..(.4) - Error in v. l*,t) x Unaltered Algorithm .05 -25(-l) 31 .01 .84(-l) IS .005 .72(-l) .9% .22(-l) 31 .80(-2) 1% .70(-2) .9% Modifications II .13(-1) 2% .68(-2) .9% .64(-2) .87. I and II .ll(-l) 1% .66(-2) .8% .64(-2) .8% (c) gh(.4) - Error in HCA) A! .05 .01 .005 Unaltered Algorithm .30(-l) .86(-2) .58(-2) I •30(-l) -76(-2) .53 (-2) Modifications II •15(-1) .48(-2) .39(-2) I and II .14(-1) .44(-2) •36(-2) As is to be expected the modifications are most effective for the larger time increments. However, even for the shorter time steps the improvement in accuracy is significant. Modification II proves to be very useful, reducing ^-u. («4)> Q_5(.4) and e^(.4) from ten to fifty percent. We remark that Modification I increases the accuracy of the approximations although WCti of (6.0) is actually slowly varying for te [o.,.4] ( h(0)=5.33, h(.4)=7.36). Order of Convergence of the Similarity Algorithm In Chapter V we showed that the order of convergence - 99 -of the Similarity Algorithm is one half. That is, we showed that for all T less than the total melting time ©4 CT) " OU^** as wi e-y -a "t • —> i Pt where ^j,:^ : . This is important in that it. explains our observation that the Similarity Algorithm should be used only to obtain coarse accuracy. In Chapter V we observed that the coefficient multiplying the order one half term of the error expansions of ?U("T) and &S{T) is usually small, hence, the effective values of and are larger than one half. Hence the Similarity Algorithm should be significantly better than an order one half scheme. Here we give some numerical examples which support that -claim. For exact solutions we use those given by Sanders £ 31 In particular, for the data fact (6, - 100 -Sanders gives the solutions J 2 J 0-**t) (6.3) on o < A < ittr- V \-«A t) * . Here M(n; t;^) are the Confluent Hypergeometric Functions, and arid A are related by the condition that: \,a is the smalle.s„t positive root of the equation M (-X0 ^ -i;^ ) - o The data (6.0) corresponds to A =.5, A0=l. We also consider the data i with various initial temperature data •U, Ur- x'-l) (6.3a) Ue0O = (x-Oe . (6.3b) V0«,XN-- x- I •. (6.3c) To the data (6.2) and (6.3) we apply the Similarity Algorithm with equal time steps, A "t , varying from five to thirty percent of the total melting time. For the data (6.2) we are able to estimate the order of convergence of - 101 -CT) -50, However, for the data (6.3) we must settle for the order of convergence of since the corresponding exact solutions are not available. In each case the algorithm is used to approximately fifty percent of the total melting time. Table 6.2 provides a summary of the results. Table 6.2 Observed Order of Convergence, Data g,(T)->o <=>V4CT)->o gy, (T) ->o (6.2) A=.5 .8 .8 .8 (6.2) A = .85403 .8 .7 .7 (6.2) A=l. .7 .7 .6 (6.3,a) - - .8 (6.3,b) - - .8 (6.3,c) - - .9 - 102 -Table 6.2 supports our claim that the order of convergence of the Similarity Algorithm is between one half and one. Moreover, it provides some evidence that the order of convergence of ;T> -> o is closely related to the orders of convergence of f>jCT)->o and _> 0 The Small Time Versus the Large Time Representation of the  Similarity Solution In this section we give operation counts for one iteration (one time step, i.e. t-i to t^, ) of the Similarity Algorithm using the representations (3.13) (large time) and (3.14) (small time) respectively. Multiplications, divisions and additions are classified as equivalent operations, while exponentiations and square roots are taken to be equivalent to twenty and five operations respectively. One iteration of the Similarity Algorithm using the large time representation (3.13) involves approximately 40 + 50 m +3n+4mn operations, where n terms of the series in (3.13) are retained and the necessary quadratures are performed by means' of a 2 m - point Filon Integration Rule. On the other hand, one iteration of the Similarity Algorithm using the small time representation (3.14) requires 2 30 + 70 J + 90 JK + 40 J + 80 J K operations, where 2 + K terms of the series for GCv,t-^)(see (6 .5)) are retained and - 103 the integral in (6.5) is evaluated using a j node quadrature rule. To initiate the comparison we write (3.13) ( "u.'^ (3.13)) and (3.14) ( u*4 <—> (3.14)) as (6.4) s&n(P<y- ^i$C' & S,«t;) cos(w^«l. and (6.5) where T^S X/C,^, } i Ac(C;4L and the notation of Chapters III and IV has been used. Our aim is to compare the number of operations necessary to calculate u* (C^4( ^, -al±+l) to a given accuracy 104 using (6.4) and (6.5) respectively. However, the -term ? in both (6.4) and (6.5) motivates us to consider instead, the number of operations necessary to calculate ^•c»+>\> A"^t-n ^ e *** ^ to a prescribed accuracy, say s 0' ^ . In each case the error enters from two sources - the error made by truncating the series and quadrature, error made in evaluating the integrals. We first consider the large time representation (6.4) by writing Integrating the expression for -^>^(j^-^ twice by parts we obtain "U,(ft>l5 -L-l ^'^ (c"/J.C^V,(c,3iai))ce,(^)cl^ ^ 0 c^ Hence it is easy to see that 00 , , where 105 -Evaluating the above integral by parts we have \ R ^ S -TT*(2n-l>Ht (6.6) Moreover, if, TT, (2^,-0 ^ * is. lftfge, enp.ugh we,, obtain the. asymptotic estimate le^(oK< ^-\H L4 e_ . (6.6a) The expressions (6.6,a) give us a "rough" estimate of the number of terms, n , of the series (6.4) necessary to achieve a pres'crrbed -accuracy for a given We now focus our attention on the quadrature error arising from the calculation of ^^^^ j^'1!'"]*1 by a 2 m -point Filon Integration Rule. Suppose -Sv-(p') is the Filon approximation to y^^.(^-), then we can write where o<t^<i (cf. Davis and Rabinowitz [9 ] p. 64). If we take w large enough so that v-j^ / %^ .$ "V^ t then the total - 106 -error due to the quadratures, call it , can be seen to satisfy where Evaluating the above trigonometric sum we arrive at the estimate Since j^ixjjx for x >, o we can write (6 Hence we have Now turning our attention to the small time representation (6.5), we write 107 *3 <3>s e '^C<,4,SU*(C;4|j;*ti+1) 11 • ^ 4-. + e Since we have a series of positive monotonically decreasing terms we obtain iestoi.« JL'rir A« e-WR,V*«-where K. = 0<H <• » (6.8) To investigate the quadrature error involved in - 108 -evaluating the integral in (6.5), we note that the dominant term of G ( C-4l ^, At;+1 j c-^) is always Q'^&'&^I . That is to say, since all other terms of G ( C^, ^, , \ ^i^) are well-behaved for °* 3 > * ' > tne source term largely determines the required number of quadrature nodes, J Hence we choose J large enough to evaluate to a prescribed accuracy and assume this J to be representative of the number of nodes necessary to evaluate to the same accuracy. To compare the operation counts we note that \TT js. 1 . Hence we set j — = i and find n, m and K so that l^C^l) of (6.6), )c^| of (6.7) and I e j | of (6.8) are each less than 10*e' for given values of . d5 La and K, - 109 -For values of S; we take .001, .01, .1, 1. (a typical range over which varies during the course of a calculation). Since K, has very little influence on the magnitude of K ( X =0,1,2) we set V<j=l. Furthermore, to assess the effect of the magnitudes of L.( and on the operation count, for*, t-hse* larger time., representation-, we3 vary both and between 1 and 100. Table 6.3 provides a summary of the results. Table 6.3 Approximate Operation Count (a) *d -=4 .001 .01 .1 1. 0 1 0 1 1 2 ,(4) 25+ 25+ 15 15 5 5 # of Operations Small Time Solution n m 27,000 75,000 9 9 11,000 28,000 4 4,000 2 6,000 1 of Operations Large Time Solution 1,000 700 300 200 14 38 5 19 2 10 1 10 # of Operations Large Time . Solution L, : I,: 100 7,000 1,700 700 700 (4) Here we have used a Gaussian quadrature scheme (cf. Isaacson and Keller [2lJ p. 327) to evaluate the integral in (6.5). - 110 -(b) cl =6 r .001 .01 .1 K l l l 2 J 25+ 15 10 5 t of Operations Small Time Solu tion n 75,000 28,000 13,000 6,000 14 36 5 18 3 18 1. 9 # of Operations Large Time Solution 5,000 2,000 1,500 700 n m 17 114 6 52 3 52 2 29 # of Operations Large Time Solution 21,000 5,000 4,000 2,000 We remark that the number of nodes given in Table 6.3 for Filon's Integration Rule is larger than was used for any of our numerical experiments (usually m was taken between 10 and 6 20). Moreover, we note that 10 is normally the limit of accuracy which one would want, since we are wasting computing time if we try to make the truncation and quadrature errors significantly smaller than the inherent error in the actual approximation to (0.3,a,b,c,d,e) generated by the Similarity Algorithm. Although (3.14) is the small time representation of the Similarity Solution, Table 6.3 shows that it is numerically impractical to use if the integral is calculated by a - Ul -conventional quadrature rule. That is, because the dominant term of G ^U,y C;^) for small behaves like a delta function in about ^ , a conventional quadrature scheme requires a relatively large number of nodes, covering the whole of the interval [0,lj, to achieve the necessary accuracy. Moreover, we remark- that in estimating J we assumed u v C,-^, <&X;)<2 * a to be a constant. Hence the number of operations given in Table 6.3 for the small time representation could be an underestimate. Furthermore, since ^ (^>,")j are independent of <it- and ^. , the cosine, terms which enter Filon's Rule need be calculated only once during the entire calculation. At each step these constitute the weights of the Filon Quadrature Rule. This is in sharp contrast with the calculation of the integral appearing in the small time representation (3.14) where at each time step "i; } G ( C;+) ^t^,' C;vj\ must be calculated at all quadrature points in ^ and in ^ . That is, if . are the nodes of the quadrature scheme then at each time "t^' GCCu,^,^.-/^^) must be calculated. These results indicate that the large time represent ation (3.13) is better than the small time representation (3.14) for the numerical solution of (0.3,a,b,c,d,e) using the - 112 -Similarity Algorithm. Table 7.3 also provides us with an estimate for the number of required terms for -the large time representation (3.13). We can see that unless -u.*"' (C; ^ > o,*t; ) is exceptionally misbehaved (reflected in the values of L, and Li ) at most twenty terms of the series in (3.1,3) need b.e us.ed. For reasonably behaved functions three to ten terms are adequate. Our numerical experiments support these statements. Comparison of the Similarity Algorithm with Lotkin's  Difference Scheme We conclude this chapter by comparing the Similarity Algorithm with Lotkin's Difference Scheme. Lotkin [ 23 ] transforms to a fixed boundary by making the transformation ^ = X/.SCi) in (0.3 ,a ,b ,c , d , e) . Then he employs centered difference approximations (cf. Isaacson and Keller [ 21 J p. 445) for the spatial derivatives appearing in the resulting diffusion equation, together with backward > difference approximations (cf. Isaacson and Keller [ 21 ] P. 445). for both -u..C»,"t) and i(fc) in the transformed version of (0.3e). The resulting non-linear system of difference equations is solved iteratively. If a uniform mesh is taken in ^ , then the above scheme is second order accurate in space and first order - 113 -accurate in time. In comparing the schemes, we use the data (6.2) with A = .5, .85403, and 1. In each case the approximations are employed to approximately ninety percent of the total melting time. In the Similarity Algorithm three terms of the series in (3.13) are used and the errors given are »?j(T) and 6vi. CT) • For Lotkin's Scheme nine interior mesh points are used and the errors given are the maximum absolute errors at these mesh points. The results are summarized by Table 6.4. Table 6.4 Similarity Algorithm Versus  Lotkin's Difference Scheme Similarity Algorithm Lotkin's Difference Scheme T=.45 A=.5 T=.22 A=l: .01 .005 T=.26 A-.85403 .01 .005 .01 .005 .32C-1) .32C-1) -(5), Computer 'Error in Error in Computer g» CT) Time(Sec) uU.-t) JCt) Time(Sec) .62(-2) .69(-2) ,57(-2) .45(-2) .24(-l) .10(-1) ,24(-l) .79(-2) .12C-1) •91(-2) .14 .24 .07 .16 .06 .12 .70(-2) .97(-2) .08 •37(-2) .55(-2) .12 .45(-l) •24(-l) .89(-l) •58(-l) •16(-1) •96(-2) •22(-l) •13(-1) .06 .09 .05 .09 (5) All calculations were done on the IBM 370/168. - 114 -It can be seen for these examples, that Lotkin's Scheme and the Similarity Algorithm give comparable accuracy for approximately the same amount of computing time. We remark that if greater accuracy is required, then Lotkin's Scheme is the more efficient algorithm. CHAPTER VII A COLLOCATION SCHEME In this chapter we consider (0v3?,a ,b;,c ,d , e) from a variational point of view in order to develop algorithms for approximating its solution. Our ultimate aim is to achieve a finite element formulation of (0. 3,a,b,c ,d , e) . The Lagrangian Equations for Heat.Conduction To initiate a variational-formuTation of (0.3,a,b,c, d,e) we follow the lead of Biot £ .1 ] by defining §f-C*it), referred to as the heat displacement field, to be the time in tegral of the rate of heat flow across a unit cross sectional area of a given slab. With this definition the equation of heat conduction can be written as |M)*5 L[§C*(t>}v ~ Uy U,t). (7.0) In addition, the law of conservation of energy is expressed by the relation IxM: - «*HU,t). (7.1) - 116 -To obtain the Lagrangian equations for heat conduction as derived by Biot [ 1 ] we first let "uC*,"i\ and §(x,t) be the temperature distribution and heat displacement field respectively associated with (0.3,a,b,c ,d ,e) . Then we consider arbitrary variations of the heat displacement field which are consistent with the conservation of energy relation (7.1) and the boundary conditions (0.3,c,e), i.e. Su(«,u-- §K<«,t) and For -any interval (a,.b) along the slab, (7.0) implies that ^ b . (1) °- ) [^^x^.tH § U,t>] SfU.tWU. (7.3) a Upon integrating by parts the first term of (7.3) and using the constraining relations (7.2) we obtain (7.2) 4) x-a. where Since (7.3) must be satisfied for all time, the limits of integration a and b can be taken to be functions of time. In fact, we will take es:o)b--i(t) ((/u,s) satisfying (0.3,a,b,c,d,e)). - 117 -The variational principle (7.4) leads to a set of equations referred to by Biot ^ 1 ] as the Lagrangian equations for heat conduction, if we assume that $(*,i> can be expressed as a given function of x and %. and at most a countable set of independent parameters (generalized coordinates) ^ ?«^^|- » i.e. $(«,tu $<•$,,... 'hen for arbitrary variations ^^^'^^j. *n t^ie Parameters \ V^\- consistent with (7.2), &$(x,i), the it* variation of the heat displacement field is given by In addition we have the relation and hence *J . (7.6) Moreover, since V(«i,b^"fcV is also a given function of the parameters ^ Ct)| % we have - 118 -(7.7) Introducing (7.5), (7.6) and (7.7) into the variational principle (7.4) we obtain f °<M i?* F + ^ r? ^ °i (7 8) where we have introduced the dissipation function Since the parameters | Ctlj, can be varied independently, (7.8) implies that **« *l. ^ « ( } - the Lagrangian Equations for heat conduction. In order to use ( 7.9) we take (w,i) to be the solution of the system of equations (0.3,a,b,c,d,e) and let iv;t,cM. be a set of basis functions for L lT«, i"] (the set of square integrable functions on C°» ll ,) with the property \>; Co\~ -v-(I) = ° for all i*vtlt... , Furthermore we assume that 5(«|t) can be written as the sum - 119 ^;(5t,). (7.10) The equations (7.9) then become the determining equations for * That iS' We takS Ck:oi'3:S(tl; substitute Vfo.itt^t^ O <°> -s<*»;t) into (7.9) and note that 5 J t) = O, S2?te,t)=o to obtain S*<t) C,^(*j • ( A,- i(.tu(t) Q,) ^<t> = o (7.11) where o .1 « o f (0=..(,|,C«,..;;^(«,...)T. Moreover, if (-u.s) is to satisfy (0.3d) then |^<t) must satisfy an appropriate initial condition. If 5 (*<°)= °<* $0**) then we take ^o") to satisfy ^ [fo^'I^^^C^l^Ctl^x^ for all y.lJ?J... ; that is where - 120 -<2> and r1 Lf-V } $0Ub) ^-txUx.-Roughly speaking, ^i(o) is the projection of 58<x) onto ^-(Vjj) • Finally the Stefan Condition, (0.3e), becomes (3) where (7 iC©)=b Even if the initial value problem (7.11), (7.12), (7.13) has a solution it would be difficult to obtain because of the coupling of the derivative terms -i(t), »0 in (7.11) and (7.13)* Hence we seek ways of reformulating (7.11) and (7.13) in order to avoid this difficulty. (2) For a given temperature distribution, Mt.Xj"fc) the heat displacement field, <f(x,t) , is not uniquely determined.. Since we are interested in *u. O,-t,) and M^(v;-t) only we take § U,-tl- - 5 ^^."t'1 c^ • As noted by Biot £ 1] the determination of JCt) is not part of the variational procedure but merely another ordinary differential equation added to. the Lagrangian equations - 121 -If in the Stefan Condition (0.3e) we use ~ <s<*\J^ for ©('v^Utt), t> instead of -§Ucu,-t) then (7.13) becomes and the initial value problem to be solved becomes - 0 where (7 . JCo)-b Another way by which we can eliminate the difficulties of (7.11) and (7.13) is suggested by the work of Biot [ 1 ]. Instead of expressing 5(.*»t) as the linear combination (7.10), we write E pitted- */J(t)). 4Ct> is, *w> define o and use (7.1) to obtain 122 5 jf <Pi<ky^ Proceeding as before we obtain the initial value problem * n S*(t> CjiplO + ( Aa - iit>.sit> QA) Jit) * O -SU>= - $ -f £ ft.(t)v?(o) + hCt)7 C4pto)= J. .S<.o)-W where ^ 0 - 123 -It is interesting to note that by taking V»v\ (x^* - (,v%-^)TTi i uw-fC'n" j)TT xj 4 (7.15) becomes the system of equations obtained by V.G. Melamed (cf. Rubinstein £30] Chapter 8) Melamed reduces (0.3,a,b,c,d,e) to a denumerable system of differential equations by assuming that oo ^ (4) • UUA)* X A^lOeoifVA-xj-Tr . Moreover, he has shown that-if an approximate solution is obtained by considering the first - ,V A/ coefficients ] Ct) ( then as /V->co the approximation converges. In other words, he takes for an approximate basis the first A/ functions | cos j[(vv-l)'TrxJ | and obtains an approximation solution to (0.3,a,b ,c.,d ,e) by solving the appropriately truncated version of (7.15). Instead of proceeding as above, i.e. using, for an approximate basis functions which are global on £o,l], we propose a finite element formulation of (7.14) and (7.15). We will start with an approximate basis, consisting of finite elements. For convenience these basis functions will be labelled so that the systems we obtain will be appropriately truncated versions of (7.14) and (7.15). In particular we wish to express uCK,t) as a linear (4) T-This motivated us to write g(.x,"t) in the form given by (7.10). 124 -combination of piecewise cubic Hermite polynomials. To this end we define a partition 77" of the interval [o,l]: and let where 9^-tlx) is the cubic Hermite polynomial defined by with x % Jf, s VCx) = and we have introduced the notation - 125 * vs > (cf.schultz [ 32 ] Chapter 3). It. is. well known.that fp„rv fixed;. "fc a», fu.nct.ion• on £o,l], which has sufficiently well behaved derivatives, can be approximated arbitrarily well in mean^^ by a linear combination of functions in ^/^(Tf**), i.e. for a set of functions > G 'U,i>*r E 9^<*) ~> g<*,*^ in mean for fixed t as V\ —> o . Again, for ^t*>'^ sufficiently smooth it is known that, for fixed "t , the derivatives G^U,t} ; ( *, t)t G ? IK, * > tend to the derivatives ^t/x.-t'i , Q*« (*,i) , ^w (v,t) respectively in mean Hence taking jv.ort.VjU) | ^..Wj •• •. <^W*! the systems (7.14) and (7.15)7 become "* The sequence of functions V,*. (OG L*C°,»] -r\z\,3t, is said to converge in mean to a function •^"(,0 6 L* C«*, il provided S - V W}'<1* -> o as v\ -> «o • - 126 -and s'cof] Q(i) + (tf,-i<tmty3,) QCt) = o5 (7.16) (7.17) respectively. Here | > /^*"\ > ' " 1' * are the appropriately truncated versions of | C ^ , A« , | , i - 1> 2. respectively and - 127 -We see that and f"^ are nonsingular since they are Gram Matrices of linearly independent functions. Thus both (7.16) and (.7...17) can be solved as initial value problems. The question of convergence of the schemes outlined by (7.16) and (7.17) is unresolved. A Galerkin Scheme We note that, because the basis functions I ^Pw^**| each have support' C XK.,, xK + ,] (Xo-o, XW4 ?= i ) t H*\,..., A/+ / y j P,, a/,; ^ (| are block tridiagonal matrices, while ^ f\ } °it } f^*-1 ave matrices. However, the estimate for •uxci(t»,t) used in - 128 -(7.17) is better than the estimate used in (7.16), since in the former we use the approximation while in the latter we use tke estimate\ We would like to achieve a formulation which combines the better approximation of -ux(s(0,-r^ in (7.17) with the computational advantage of the sparse matrices of (7.16). How to proceed becomes apparent if we consider the first equations of the systems (7.14) and (7 .15) , ..r.espect ively , i.e.. and Mr1"1 Integrating by parts the middle term,then substituting the appropriate express ions for 5 ("i"^ > 5XK C><i't\, we obtain 129 and respectively. That is, the Lagrangian equations do no more than force the heat displacement field j>(x,"fc) to satisfy the heat equation in a weak sense. This suggests that we forego the introduction of S§C*i"t} and instead work directly with v<Cx,-fc) . That is we write where now is a basis for L*C«>(il with x>J(o) =• v. OV- o ia-i, a, • • . and the parameters, are to be determined by the Galerkin conditions 00 I 1 ^ ( S V< V-(X) dJt ) cl; <o) r ^ Vo(V>0 v In other words, we force vt(.x,-t) 'to satisfy the heat equation in a weak sense. The system (0.3,a,b,c,d,e) then becomes where - 130 -* iUt) C efct) 4 (A--i<k)i«t)Q)d(i)so o o M;> = W«>vj<*>^> (7 <2>: .«=.»,-a, Now to express this in terms of the finite element approximation we take Hence (7.18) becomes ict)5U)- «* (,&A/*»,»<t) - JU)K(il) ^ (7 TOCo):' Oc b where - 131 -0 ((SC, C^il, »"r>' - 132 -fckls* W*U) V°cU s ^ a,...,*/ i * o, D(i)~ (0,,u>, 0A,ct),..., D^c-t), D„4,,x(t»T = ( ., cl3A/(t))T and we have introduced the notation a - 133 -The approximation TJU,^). to V-(x,i) can be written as UC*,t> o„ct> Hence the variational principle (7.5) has led us to a semi-discrete or continuous Galerkin formulation for the Stefan Problem {0o3-,a-,b,c ,d,e). That is, the spatial variable of the system of equations•(0.3,a,b,c,d,e) has been discretized while the time variable remains continuous. We remark -that the system (7.19) combines t.h.e sparse matrices of (7.16) with the better approximation to M„ <t},-t) of (7.17). Moreover, unlike (7.16) and (7.17), the solution of (7.19) yields directly and hence most accurately approximations to the quantities of interest ^ u.Cx,-t^, U* U ,-k>^ . It should be noted that if "ue OO <j C'Oi^lj then an initial value of "OCtV can-he, obtained by interpolating the initial condition "U^Cx^ at the points \ X. { of the partition TT" , i.e., O^, (oV- "U/b*^ , Oj, loi= b -U0(bXj) ^=',* instead of projecting "U0(>o into j^^ilr") . The order of accuracy of the approximation U"lx,t) to \*. (x|t) is not affected. Since is a Gram Matrix of S. A/ linearly - 134 • independent functions over [u,l] it is invertible and hence the system (7.19) is solvable locally in time. In fact (7.19) has a unique solution for all i. such that •S(\t>c, To show this, we define and note that which upon integration by parts becomes ' ° o that is, Hence multiplying the first equation of (7.19) by DT(t) we obtain i[6T(t)TD(t>l 4<tW cU(t). fDT(t)^ 6(t)l - - 1 0(t) ' dt cli 1 J i(t) or cL ns<t)-"5T<t>^D(i>l - - z DTCt)6( 5<t). cli J JCi) Since ,S(M>© and «< is positive definite, we have - 135 -3<t) 0Tlt>Y*Dlt> $ biSJVo,, , We see that if (.*)-><*? as t increases then Dct) —> o . This leads to a contradiction since the second equation of (7.19) then implies that <3(t)$o' as t increases since V\ (A) >, o, Hence -S(i) remains bounded. If D It) -> °o then we must have that J ct) —> o . Hence we are able to conclude that (7.19) has a unique solution for all i. such that jtt)>o. A Collocation Scheme Because of the properties of \ > , (7.19) is a convenient system to analyze. However there is another (Collocation) formulation closely associated with the continuous Galerkin system (7.19) which is more convenient for numerical computation. To introduce the Collocation scheme we write the Galerkin Conditions (7.19) as and integrate the middle term by parts to obtain - 136 -- A/4/, rf= 2 , Using a two point Gaussian scheme to evaluate the integral over • C*v\, Xw+il we obtain A/ 2 O 'S, *t hHi) wt<'?«•t, - w>" f*»Si*>- ?«*«>*<t> w,(v£t>? ^W.jjv, where Now if - 137 -then the first equation of (7.19) -is satisfied to OC^M It has.been shown that a continuous Galerkin scheme, using cubic Hermite polynomials, produces at best O ( W¥) approximations to solutions of fixed boundary parabolic systems (cf. Douglas and Dupont [ 11 ] ) hence greater accuracy in evaluating the above- integral- is- in* a? s/en-s-e; "wasted". Thus we have the following Collocation scheme for (0.3,a,b,c,d,e) ittViCtV= ( wyOjti - J-lt)W itv) y (7. wi^,o) * uecb^) ^i,...^ *S<©)= b It is clear that (7.20) can be written as the system - 138 -•sVt)^J CK^ + + i<usct> B) Bet) 0 J 1 3(o^ - b where A-'f% e, \ ft, (6) For a collocation formulation based on Cubic Splines see Doedel [lO ]. 139 /V-A: a; 6" »K KM, - 140 -A K A/- ; (51 e: >?,K<,07l) \ To show that the system of equations (7.21) can be solved locally in time we must show that is invertible. The following argument, due to Douglas and Dupont £ 12 In [ 12 J it was demonstrated that collocation schemes such as (7.21) provide O t V\*) approximations to solutions of a large class of second order non-linear parabol differential equations. - 141 -demonstrates this. The proof proceeds by contradiction. Suppose there exists a nontrivial vector such that b - o . (7.22) Upon constructing the piecewise cubic polynomial " i .we-see ..that (7.22) .implies that 2<>?J>-0 , ^t^...^/. (7.23) Since 3J(l)*o we see that £ (x) vanishes at three points on L^A/>YA/-*»1 • Hence either I: ?(x)*0 on CX^I] , or II:' < o . Since the piecewise cubic 2(^)€C,C°/l] > it is clear that I together with (7.23) implies that j?(x) = 0 on [o,l]. Since la^o this is impossible and hence II holds. - 142 -From Z(X*) Z'(*~)<o and ? (J?,"*'^ £ C1?'') - O we have that the piecewise quadratic must vanish in l"') and ^f'j ' hence ifc follows that Continuing inductively we obtain that *!(O£Yo)<0 which is -a., con trad it Ion since . Hence the original assumption is false and is nonsingular. We remark that there is little difficulty in generalizing (7.21) to include more general one dimensional single phase Stefan problems. However even for the simplest ease of (7.-21) the convergence question is a diff icuIt -one. The difficulty arises from the extreme stiffness of the system of equations (7.21). Moreover, the non-1 inearity is such that the methods of Douglas and Dupont £ 12 J are not easily applicable. We conclude this section with the remark that the piecewise Cubic Hermite functions were chosen for convenience. No doubt a wealth of systems similar to (7.21) can be obtained by choosing bases constructed from other finite element funct ions. Numerical Results Here we give numerical results for the Collocation scheme (7.21) only. We do this since for parabolic equations on fixed spatial domains, Collocation and Galerkin schemes, based on Piecewise Cubic Hermite Polynomials, have the same order of convergence, provided we collocate at the Gaussian points ' Si.nce we have no reason to believe that the Galerkin scheme pro videos a., bj&tker. e;s4t,ima>t;e> fo.r "Uk^A-^'h , and hence itt) , than does the Collocation scheme, we adopt the latter on the basis of computational ease. Computationally the Collocation scheme (7.21) is much easier to implement than the Galerkin scheme (7.19), since the former requires function evaluations where the latter requ-ires quadratures. Moreover, the "bandwidth oJf the Collocation Matrices ^&><&J *s tour while that of the Galerkin Matrices ^V^^C^y^! is six. Hence the system (7.21) is computationally more efficient than the system (7.19). To solve the initial value problem (7.21), we must contend with its stiffness. That is, the condition number of the matrix "o^"'^ , arising in the Galerkin system (7.19), increases as , hence we see that the time constants present in the solution of (7.19) have radically different magnitudes for large N . The intimate relationship between the Galerkin system (7.19) and the Collocation system (7.21) leads us to - 144 -suspect that the system (7.21) is also stiff. This is substantiated by the following numerical experiment. We employ two numerical procedures to solve (7.21). The first is an Adams-Basford-Moulton Multistep Predictor-Corrector Method (cf. Isaacson and Keller [ 21 ] p. 388) while the second is a Multistep Predictor-Corrector Method due to Gear £ 18 J constructed specifically for stiff systems. In all cases tested, we have found that the time step required to maintain a given accuracy in the solution of (7.21) was much larger for Gear's Algorithm. Consequently Gear's Algorithm was •found to execute five to ten times faster than the "A'dams'-B'asTord Moulton Algorithm. We conclude from this that the system (7.21) is indeed stiff. In what follows we adopt the notation: where {M>0 is the solution of (0.3 ,a ,b ,c , d , e) , (U>?) is the - 145 -solution of (7.21) and it is understood that v. (.*,t) = MW)t (x,"D ' o for J< > *<.i,\ . To illustrate that the Collocation Scheme provides accurate approximations to the solution of (0.3,a,b,c,d , e) , we solve the inhomogeneous heat equation ••"M^ cx/O + *"lx»t) o < x < 4tt> y \ >o (7.24) together with the boundary conditions (0.3 ,a ,b ,c ,d , e) . We obtain our first example by setting b = ) j and picking "f"l*,i.)> W (i), and w0 (*) so that the solution becomes ( (7.25) We do this for A =10, 20 and 50. In each case we take uniform spacing and use the approximation to T = .4. Table 7.0 provides a summary of the results. To deal with the inhomogeneity, V***"!) > tne obvious changes are made to (7.21). - 146, -Table 7.0 Errors in lUx.i ) , "Uy (,x,-t), TAX« and JCt) (Exact Solution (7.25)) N e,(.4) e,(.4) e^4l A=10 3 •30(-l) .30 (-1) 5.2(0) •58(-3) 5 .20(-2) .95(-2) 3.5(0) .80(-5) 7 • 55(-3) .44 (-2.) 2.3(0) .60(-5) 9 .13(-3) .14(-2) 1.4(0) •12(-4) Observed Order of Convergence 4.9 2.7 1.2 A=20 3 .25(0) .30(0) 6.5(0) .40(-l) 5 •ll(-l) .22 (-1) 9.5(0) .19(-4) 7 .24(-2) .14(-1) 7.1(0) .20(-5) 9 .80(-3) .73(-2) 4.8(0) • 25(-4) Observed Order of Convergence 5.2 3.3 1.2 A=50 4 .90(0) 1.10(0) 18.3(0) .15(0) 5 .25(0) .38(0) 18.2(0) .41(-1) 7 .16C-1) •30(-l) 24.2(0) .62(-3) 8 .94(-2) •32(-l) 22.6(0) •36(-3) Observed Order of Convergence 6.9 6.5 --From Table 7.0 one can see that good accuracy is obtained inspite the large values of the spatial derivatives of the temperature distribution near x-o . - 147 -Next we consider (7.24) with b-l » a^*=l and •^tx,*.); W it) t u, lx> chosen so that the solution becomes VUx/k)* cosx ( x*-Approximations to the, above solutio.n- are obtained- fo,-r. S^. =100, 500 and 10,000. We are interested in the accuracy of the Collocation approximation when \j(t)| is largest, hence we employ the approximation until approximately 20% of the slab remains. In each case the 'partition Tf" is taken to *be •"uniform. The results are summarized by Table 7.1. Table 7.1 Errors in UCx.-tl^.tx.-t), u«,(x,t) and si-t)  (Exact Solution (7.26)) N eB(T) et(T) e?(T) e.,(T) B=100 3 .13(-3) .80(-3)~ .11(0) .96(-4) 5 .19(-4) .ll(-3) .43(-l) .16(-4) .14(-1) .60(-5) 1.9 (7.26) 9 .40(-5) .24(-4) Observed Order ^ 3 2 of Convergence - 148 -B=500 3 .14(- 3) • 14(- 2) .12(0) .97(-4) 4 .39(- 4) .48(- 3) •67(-l) •23(-4) 5 .13(- 4) .16(- 3) •43(-l) .80(-5) 6 .18(- 4) .90(- 4) •30(-l) .80(-5) Observed Order of Convergence 4.6 4.1 2.0 B=10,000 3 .17(- 3) .84(- 2) .17(0) • 80(-4) 4 ,46(- 4) .25(- 2) •67(-l) .22(-4) 5 .25(- 4) .11(- 2) .47(-l) .10(-4) 6 .16(- •4) .45 (- 3) .31(-1) ..60.(-.5) Observed Order of Convergence 3.4 4.2 2.4 Table 7.1 shows that good approximations for both J(i.) and \*<.*,*) are obtained although is relatively large during the periods of time under consideration. Tables 7.0 and 7.1 indicate that the Collocation Scheme is of high order, however, the results concerning the order of convergence are scattered and do not give a good estimate for the actual order of convergence. This is not surprising since the order of convergence is a characterization of the asymptotic behaviour of the error. Hence accurate estimates for the order of convergence can usually only be obtained when fine partitions 77^ of [o,i] are taken. - 149 -The above results indicate that the Collocation Scheme (7.21) can be used to obtain accurate approximations to the solution of (0.3,a,b,c,d ,e) . However, it should be emphasized that the stiffness of (7.21) makes it a numerically inefficient scheme. Hence until a method is devised to solve (7.21) efficiently, the utility of this scheme, is in dqu,bt. CHAPTER VIII CONCLUSIONS This- thes-is has; pr-esent.ed two- algorithms- for the numerical solution of the Stefan Problem (0.3 ,a ,b ,c ,d , e) . We have seen that the Similarity Algorithm provides us with a reasonably efficient method of obtaining "rough" approximations to the solution of (0.3, a ,b ,c ,d,e) . Moreover, for the practical situation of both a smooth initial telnp'ter'a'tu're distribution and a constant heat flux, the. Similarity Algorithm promises to be an efficient algorithm. Furthermore, the heat flux generated by the algorithm can be used both to improve the accuracy and to estimate the error of the approximation. We observe that if a very accurate approximation is required the Similarity Algorithm is not an efficient algorithm. The large number of terms and the small time increment required to achieve this accuracy results in long computation times. The Similarity Algorithm is the direct result of applying the Similarity Method to a system of differential -'150-- 151 -equations. The above procedure illustrates how the Similarity Method can be used effectively to obtain approximate numerical solutions of non-linear problems. While the Similarity Algorithm gives "rough" accuracy, the Collocation scheme is capable of achieving high accuracy. Although both, the. S.imilar.ity,. Al gpr;i£"hm> and .Lp.tltin-.'s . dif scheme execute faster than the Collocation scheme, we have found that for relatively coarse partitions, Tr" , the Collocation scheme achieves accuracies which the other schemes cannot "attain. We have seen that the apparent stiffness of the system of ordinary d iff erential equations (7.21") is 'the caus'e'of th'e inefficient performance of "the Collocation scheme. We conjecture, that the simple form of the non-linearity appearing in (7.21) will allow us to construct a scheme which deals with the stiff ness of the equations in an effective way. We also conjecture that this Collocation Method can be used to deal with Stefan Problems involving both more general boundary conditions and more general governing parabolic differential equations. BIBLIOGRAPHY M.A. Biot, Variational Principles In Heat Transfer, Oxford University Press, Ely House,London W.l, 1970. G.W. Bluman, Construction of Solutions to Partial Differential Equations by the Use of Transformation Groups, Ph.D. Thesis, California Institute of Technology, 1967. G.W. Bluman, Similarity Solutions of the one-dimensional Fokker-Planck equation, Int. J. Non-1 in. Mech., Vol. 6, pp. 143-153, 1971. G.W. Bluman, Applications of the General Similarity Solution of the Heat Equation to Boundary Value . Problems, Quart, of A. Ma., Vol. 31, No. 4, pp. 403-415, 1974. G.W. Bluman & J.D. Cole, The General Similarity Solution of the Heat Equation, J. of Math, and Mech., Vol. 18, No. 11, pp. 1025-1042, 1969. . G.W. Bluman & J.D. Cole, Similarity Methods for Ordinary and Partial Differential Equations, Springer-Verlag, New York, Inc., New York, 1974. J.R. Cannon & C. Denson Hill, Existence, Uniqueness, Stability, and Monotone Dependence in a Stefan Problem for the Heat Equation, J. Math. Mech., Vol. 17, No. 1, pp. 1-19, 1967. R. Courant & D. Hilbert, Methods of Mathematical Physics Vol. I, Interscience Publishers, Inc., New York, 1953. P.J. Davis & P. Rabinowitz, Numerical Integration, Blaisdell Publishing Co., Waltham,Mass., 1967. E.J. Doedel, A Collocation Method with Cubic Splines, Manuscript, University of British Columbia, Vancouver, B.C., 1973. - 152 -- 153 -Douglas, Jr., & T. Dupont, Galerkin Methods for Parabolic Equations, SIAM J, Numer. Anal., Vol. 7, No. 4, pp. 575-626, 1970. Douglas, Jr., & T. Dupont, A Finite Element Collocation Method for Non-linear Parabolic Equations, Manuscript, University of Chicago, Illinois, 1973. Douglas, Jr., & T.M. Gallie, Jr., On the Numerical Integration of a Parabolic Differential Equation SubjecJ: to aa Mo><g> B.ouj^da^y. Condition, Duke, Math. . J., Vol. 22,' pp. 557-571, 1955. Fasano & M. Primicerio, Convergence of Huber's Method for Heat Conduction with Change of Phase, ZAMM, Vol. 53, pp. 341-348, 1973. N.G. Filon, On a Quadrature Formula for Trigonometric Integrals, Proc. Roy. Soc. Edinburgh, Vol. 49, pp. 38-47, 1928. Friedman, Remarks on the Maximum Principle for ••Bar-abol ic •Equations and Its Applications, Pacific J. Math., Vol. 8, pp. 201-211, 1958. Friedman, Free Boundary Problems for Parabolic Equations I. Melting of Solids, J. Math. Mech., Vol. 8, No. 4, pp. 499-517, 1959. W. Gear, Numerical Initial Value Problems in Ordinary Differential Equations, Prentice-Hall, Inc., Englewood Cliffs, New Jersey, 1971. Hellwig, Partial Differential Equations An Introduction, Blaidell Publishing Co., New York, 1964. Huber, Hauptaufsatze liber das Fortschreiten der Schmelzgrenze in einem Linearen Leiter, ZAMM, Vol. 19, pp. 1-21, 1939. Isaacson & H.B. Keller, Analysis of Numerical Methods, John Wiley & Sons, Inc., New York, 1966. - 154 -H.G* Latidau, Heat Conduction in a Melting Solid, Quart. J. Appl. Math., Vol. 8, pp. 81-94, 1950. M. Lotkin, The Calculation of Heat Flow in Melting Solids, Quart. J. Appl. Math., Vol. 18, pp. 79-85, 1960. J.C. Mason & I. Farkas, Continuous Methods for Free Boundary Problems, Proceedings IFIP Congress 1971, Numerical Math. Section, pp. 61-65, Ljubljana, Yugoslavia, 1971. V.G. Melamed, Stefan's Problem Reduced to a System of Ordinary Differential Equations, Izv. Akad, Nauk SSSR Ser. Geofiz., pp. 848-869, 1958. E.A. Muller & K. Matschat, Uber das Auffinden von Ahnlichkeitslosungen partieller Differential-gleichungssysteme unter Benutzung von Transformationsgruppen, mit Anwendungen auf Probleme der Stromungsphysik, Miszellaneen der Angewandfcen Mechanik, Berlin, pp. 190-221, 1962. L. Nirenberg, A Strong Maximum Principle for Parabolic Equations, Comm. Pure and Applied Math., Vol. 6, pp. 167-177, 1953. L.V. Ovsjannikov, Gruppovye Svoystva Differentsialny Uravneni, Novosibirsk, 1962. (Group Properties of Differential Equations, translated by G.W. Bluman, 1967). L.I. Rubinstein, On the determination of the position of the boundary which separates two phases in the one-dimensional problem of Stefan, Doklady Akad. Nauk SSSR, Vol. 58, pp. 217-220, 1947. L.I. Rubinstein, The Stefan Problem, AMS Translations, Vol. 27, Providence, Rhode Island, 1971. R.W. Sanders, Transient Heat Conduction in a Melting Finite Slab: An Exact Solution, ARS J., Vol. 30, pp. 1030-1031, 1960. - 155 -M.H. Schultz, Spline Analysis, Prentice-Hall, Inc Englewood Cliffs, New Jersey, 1973. P.V. Soloviev, Fonctions de Green des Equations Paraboliques, Comptes Rendus (Doklady) de l'Academie des Science de l'URSS, Vol. 24, No pp. 107-109, 1939. ' APPENDIX A Lemma 1.1 (Friedman [l?])• Let £>Cfc) be continuous on the interval Co,1?! . In addition, let set) satisfy the following Lipschitz condition for all t,,t4€C°>°"3 and some constant M . Then defining \ /MT) CK(X,t: S(T),T^J d T r1 SK 3C-4Ct) Ct) KCx,ti 5<r),r) we have C>t -1 /5 (t) , X-S-SOW-O X Proof. Before proceeding we make the following definition: rt I(*U XVCT) t*'-*(Y» KC/tt^cy),r) dg - \ J\r) Lliill-iilll K (S(-l), t; 5Ct),t> dr J, 2Ct-t) - 156 -- 157 -where ^ is any continuous function and % e (o,t). Now if we can show that (A-l) r Ji—. | Kp^ + L + -L/>Ok)| <• as % ->o . where. i _ ^P(t) (x-j(y))K(x,t-.iWl^df UW.-3»-tt K (Jit), tj5C*VT)d2• J Hit-*) o then the statement of the lemma follows by allowing <B to tend to zero. Further note that since (A-2) (A-l) follows if we can show that To x->stt)-o 2 %->o (A-3) I IC«)| < OU) see this we form the expression p < *)i p Civ- (pti^-^o <r>) - 158 -and substitute it into (A-l), to obtain where To establish (A-3) we write 1(0 = I, + I. and consider c Since -SOf) is Lipschitz continuous and le , - .(g..." M .5_J for all real, we obtain P»l * ^ (A-4) To estimate I( define - 159 -and consider *-* (A-5) {lx -5tr»*-{x - 4 «»*j /V (t-If we take and "M" we see that the expression in the exponential can be bounded by one since <: K\ IX.-.SC-OI + M1 & $ ] 2 & Then using the inequalities I t - <c~*l < 3 ^ for \£) $ I and 160 Z ¥ (A-5) becomes-IT,-1,1* [S fA S*^. (A-6) 8 n''*. To evaluate ^T, as x-> ,stV>-o, we let Then where $> " /(x-SCtn • From (A-7) we conclude that cj T, = - -jj • (A_8) To complete the proof of (A-.3) we note that (A-6, 7, 8) imply - 161 -(A-9) Hence writing 1- 1 I, + ii I we obtain 11(01$ ft M 5V* ' +. JL . . (A-10) 2 Tr"* .2 Thus x->Jtt>-o X->JC«-e> Ti*a and the proof is complete. APPENDIX B Proposition^^ 1.1 The (weak) maximum principle. If vu*,-t) satisfies with .St/O a positive continuous function and /UCx.t) 6 C C-B) t -M^U.t), VtU,t) ^ C ( «£> U (g^) where is the closure of and. ST - I <*,T) : O < X < 6 C.T)| then wU,t) attains its maximum and minimum value on the data boundary - iQy . (See Fig. 1.0 Chapter I) Proof!" For any €>o define v<*,-t)= v. Ix,-U - •£ t • where vttXj-fc) satisfies the hypotheses of proposition 1.1. If vt*,t) assumes its maximum value at a point Cx.,i,> e °& U @T (1) See Hellwig [ 19 ] P- 47. <2> For more general results see L. Nirenberg [ 27 ]. - 162 -- 163 -then vx< (»,t) is defined and continuous on 5 C ^ 'Sj where for some positive h and K« Furthermore, C*, ,!,')»< o and hence (.x, .jt.Y 5,,- € . By the continuity of vt^l-*,-!') and hence V^Cx.-t) we establish the existence of a S€ CO,K) such that *.x,,\>$ - £/z for all Is <t,-fc,i,7 Hence and vtx.t) cannot attain its maximum value on ^ that is e do©~<g>r Now suppose v<C.xl't) attains its maximum value on odU*ST at t>*,t<,V such that i.e. The maximum value is not attained on the data boundary JJd'&j % Then - 164 -i Thus the maximum value of v(*,i) is attained at Cx,,!,), -t,<i.„ However, in particular *w.C*,,t,) >, -uCXo,!^ -«r(te-l,) and hence for £>o but otherwise arbitrary. Noting that t|<"t0 for f>© and 'allowing f to tend to zero We 'have which is a contradiction. Applying the foregoing argument to -u()r,-fc) we obtain the conclusion of the proposition. Proposition 1.2. If (v(.,,s) is a solution of the system of equations (0.3,a, b, c, d, e) then vtxcscu,t^ >o Proof. (3) For more general results see Friedman [l6J - 165 -If for some "t , -wx(stt>,ii < o then by (0.3b) there exists X0 6 Csuch that v\(xe)i)^ h > o for some & .» Now define i€D».T] t $ and note that if *c > o then the Maximum Principle (Proposition 1.1), together with conditions (0.3b, d) implies that there exists a "tv(S>) $ t! such that that is v^(o,t(6))- &• > © • To show that this is impossible, define for any £^£0,%] and note that -u. (o, "t (t)) - £ is the maximum value of -vi.(.*li.) on S£ - jCx.t): Of xj5CU, oj-t ii(tl Then since ( o + , •fc') - MX o we conclude that Ux^o+ , ttn) c M (ov.Ho) { O for all ^c-Co.J,! , Thus - 166 -that is "U. (o, t(.b)) .< O which is a contradiction. Hence the original assumption is false and thus v\x (sct»,ti z,o (Note that by extending the above arguments we can show u(o,i) *o for all L>A3 ). APPENDIX C To show that as «f-?o (1.1) becomes (1.2), we write (1.1) as where Vt =. ^ Vv(^( v-£) <5+(x(t;?) t-f) d£; o (c Then considering V,, va ,Vj separately we take the limit as £-> o . Since G+(x, t • has an integrable singularity at T=t , we see immediately that ofV3r - V V^SC^t) G^x.t; SC*),*) . (C o To find efCwv V, we write - 167 -- 168 -where o .sen V," - ^ j><^) -"U.tO] G*0,t:s-,c) d£ , V/M: - ^ G\x,t;^,o) . From which we have G^x.tjf.o) j J(f). Since for fixed t.>o, G+(x, ^, 7:) is continuous at r= O we see that 3l v.'-O . (C-3) £->0 It is easy to see that - 169 iv,"l 5 -*"P 1 "K(£,d> ~-U0(£)\ hence Caf V/' r O , (C-4) Finally for fixed i > o, V,'" can be written as where -S If) $ * £ b from which it can be easily seen that cj wv^ V,"' - o . (C-5) Combining (C-3, 4, 5) we have €->© o Now for V, we fix S f (o,x) n (o, i tt-n-and (C-6) wr ite > as £ tends to zero the first and last integral, as well as the - 170 -K(xj'tj-f -t-£) part of the second integral do not. contribute since -3T as £->© . Hence we are left with ^ -u.C;5;i-.£) l< C*.,. x-* After making the change of variable V- Villi) in (C-7), we consider fV"V c/vi—\\ \ e"Vu(.x-at'/'v>i-cUlvI - v<.W,i) . '«->© L TT '* ) -J Since A*Cx,"k) is continuous in and hence oQ£ , we can write (A-8) as for some x (.*-$>,x+ Letting S tend to zero, we concluc Combining (C-2, 6, 9) we see that as £->o (1.1) becomes 4 ^ vVl-S«t>.t) G^tXjt-.-K*),*) At APPENDIX D Lemma 1.4 (The equivalence of the differential and integral systems). •If *u(.t) is a solution of (1.4), where JC-tvis given by (1.5).., then C/M,S> ( ) defined by (1.2) and -sc-f) defined by (1.5)) forms a solution of (0.3,a,b,c,d ,e). Conversely if <x*,s) is a solution of (0. 3 , a,b , c , d , e) then v>(t^ = VK (•sci^t) is a solution of (1.4). Proof. By construction any solution ("u.s) of (0. 3 , a ,b , c , d, e) defines a solution V- u„i-sit\t) of (1.4). Conversely if viC-tV is a solution of (1.4), (1.5) then in the region oQ the integrals of (1.2) are regular. Hence they can be differentiated directly to show that va^t) (defined by (1.2)) satisfies (0.3) in <=& . Furthermore, since G+(.x,tj T} is an even function of x. we see that condition (0.3c) is satisfied and differentiating (1.5) we have that condition (0.3e) holds. ; Evaluating cb - 171 -- 172 -shows that condition (0.3d) is satisfied. To demonstrate that condition (0.3b) , holds , we note that we have shown that •u.c.x.'O (defined by (1.2)) satisfies (0.3,c,d,e) on JQ . Hence we can integrate Green's Identity with wU,t) (defined by (1.2)) over and use condit ions (0.3,c,d,e). After taking the limit as €->© and subtracting (1.2) from the resulting expression we are left with O 6 - \ G^- (x, tyJt*),T)d»-o using the relation Of* txY~ - Gx Cx;t; w.e have. O- 2- \ v<C5«),t) G"( x,t; 3tt»,T)cl2-dx ) o + ^ •uC.Strj.t) ic«r) G4Cx,t3 5ltj,r) d> Letting x->JCt>-o and using Lemma 1.1 (D-l) becomes O- vt-C-SCt^t) + X } u(str\,t) [G'CSCDjt^itt),^) o (D-l) (D-2) Since s$ Ct) is Lipschitz continuous, from (D-2) we deduce that Ivastt),-fc)| satisfies the inequality - 173 -VM\ $ • ^ l-ul * ' d* o where ^" has at most an integrable singularity. Using an inequality of the Gronwalltype we conclude that u. (ic-t^-fc) u.o and thus the lemma is proven. (1) If where ° then APPENDIX E THE FIXED BOUNDARY SOLUTION We wish to solve the following system of equations. x « 1 >o o To obtain a solution of (E-l) we let v.(<(t)= w < *,t W v (x, where wlx,Vi satisfies Wx0,"k^O- V >o - 174 -- 175 -and v(*,t) satisfies Vx (o,t^ o ^ > o \ . (E-3) The solution of (E-2) can be written as W tx where 4-<5 (E-4) To construct the solution of (E-3) we take the Laplace transform of the equations (E-3) and obtain -U„ (o,p) * e •^x 0,p) = HCp) (E-5) - 176 -where The solution of (Er5) can be written- as *U(*,p)- H (p) Co<V\Jp_x. 4 p 5w>.\\Tp and hence the solution of (E-3) is given by vff,tv- -i- \ M ( pVP Pt c P-S^^P K), d p (E-6) V-toa Finding the inverse Laplace Transform of \ cas^n) ) by expanding the same for large p we see that V(*,t) can be expressed as the convolution integral - £ ||e ( e"°+W/*(t-» + e-(,+aw+**cfc.»))| 6\ (E-7) It is easy to see that (E-7) satisfies the heat equation on the domain | UAV. x t (t>,A ,t >t)| as well as the conditions v^c^so and ••v„('b,t)-0 . To see that v„ U,t)= V\(fc) - 177 -we note that *:t"« For fixed f> C V - x > we can take \l-xl small enough so that (\-x) / c for any prescribed £>o. Hence x-Jt"« •n*k x->i 1 *f*i where € >'« , S> > o a"d MlWm-KU)| . If *t UG>,VI a continuity point of Y\(i) we see that - 178 -as Hence the solution u.(«,t) of (E-l) can be written o (the small time solution). APPENDIX F THE ASYMPTOTIC EVALUATION OF FOR SMALL t To obtain an asymptotic approximation of +1)^+ t ) we differentiate (3.14) with respect to x and evaluate v'„(x,t) = ^ 'W'CJt.At;) G„(«,tic;> el? where GO t) (F-l) - e asymptotically for small -t at x* -^-t. - 179 -- 180 -We have, up to exponential terms, Hence we must obtain an asymptotic expression for Upon making the substitution ,JC £ (•" (F~2) becomes Substituting the expression •f o into (F-3) and using Watson's Lemma we have / - 181 3¥ and hence + OCO. 


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



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"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items