LIBRARY OF THE UNIVERSITY OF ILLINOIS AT URBANA-CHAMPAIGN 510 .fc* TS
0. iii. The error p_(t, e) = u - Au in approximating the partial derivatives by finite differences is controlled at t. such J that AAt. I |a(t . , e) I I < e/2 where At . = t . n - t . . In J — J 2 j j+1 j practice, this is too much to expect. At the end of this section, we discuss the effect of relaxing this requirement to make it more realistic. Let A. denote A(t., e). In the event that A( t . , e) 4 A(t , e) define J J J J -*- A(t" e) = A(t e) = A Let v denote the exact solution to v. = XA(t, e)v (2) — t — We obtain approximations w. to v. = v(t., e) by numerical integration using tJ J J the trapezoidal rule ' *i+i = v -i + 1/2 Xit J !A(t j^ + "Vi'^n 1 = w + 1/2 XAt A(t )[w + v ] (3) J J J J J •*- Define t. to be the local truncation error of the trapezoidal method in (3). -J That is, if w, = v., then t. = v ! _ - w... Assume for the moment that we -J -J -J -J+l -J+l can control the stepsizes {At.} such that for all j £ , x. L. = e/2. Later, J ' ' - J ' ' 2 we will examine the effect of slightly weakening this assumption. Let z = v - ¥ be the difference between the exact and numerical solutions to (2). We will now derive a bound for z L. , then address ourselves to the question 1 '— n' ' 2 of bounding |u - vj | , and lastly combine the result of the two analyses to arrive at an estimate of the global error | |u - w| | for the program. In solving (2) by the trapezoidal rule, each time a step is taken, the current global error is multiplied by exp(At.XA.), and a new local J J truncation error is added. We thus have z.,_, = exp(At .XA. )z. + t. (h) -J+l . J J -J "J Assuming that w. - V-. then z = 0, and applying (k) , we have "0 — - "0 — n-1 n-1 z = i 1 - i • TLi exp(At .XA. )]t. Since the A. ' s are symmetric negative definite, I exp(At .XA. ) I I _ $ 1. Using l ' ' J l 2 this fact and the assumption that It.IL = e/2, we obtain N^ll 2 « °Jo e/2 " ne/2 (5) We now examine |u - vj | Q . Using the above definition of o_(t, e), we can rewrite (l) as u. = XA(t, e)u + Xa(t, e) (6) — t — — Subtracting (2) and integrating, we have t t |u(t ) - v(t )|L $ |/ n exp[X/ n A(s', e)ds']Xo(s, e)ds|L — nn^Q s — d We again use the fact that A(s', e) being symmetric negative definite implies t ||exp[X/ n A(s', e)ds']|| * x to S et s c - t t |u(t ) - v(t )| L $ / n I |exp[X/ n A(s», e)ds']||l |Xa(s, e) | Lds — n— n^Q s d — d t $ / n ||Xa(s, e)||_ds d We approximate / J | Xo_( s , e) | | ds by At . X | |o_(t . , e)|| . Assuming that the t ^ ^ J error in this approximation is small enough to "be neglected, we can use the assumption that XAt.||o_(t, e ) | | <: e/2 to get J <- n-1 t A +1 u(t ) - v(t I < e/2 .1. f J ds/At. = ne/2 (T) '- n — n '2 J=0 t J J Finally, we combine (5) and (7) to get a bound for |u - w| | . Using the triangle inequality, we get |u - w| | $ ne (8) Since we are using the trapezoidal rule to integrate (2) numerically, it is reasonable to expect that Ix-llp = 0(AtT > ). But, since | t_. | | = e/2, it 3 1/3 follows that At. = 0(e) and hence, that At. = 0(e ). If we now assume that J J -1/3 At. = 0(l/n) (as it ought to be), then we have n = 0(e ). Substituting J this relation in (8), we get |u - w|| 2 $ ke 2/3 (9) for some constant k. In practice, we cannot expect a program to actually control o_(t., e) J and x_. as strictly as was assumed in the above analysis. We will now look at the consequences of reducing these expectations to a realizable level. Since u is smooth and A(t, e) is piecewise smooth in t, we can express the truncation error of the trapezoidal method at t . , as a power series in At. J+l J 00 x. = e x. .At: -j 1=3 -tj Instead of requiring that ||t_JIo = £ /2 5 we require that I|t At.|| $ e/2. In practice, we do not know t_ exactly, but we will assume that we can estimate it accurately enough so that the following conclusions are not altered. Since | t_ A"t . | L $ e/2, we may infer that At. = 0(e 1/3), and -30 J u hence | 1 t_ | | = | |t At^| I + o(At^) * e/2 + o(e) (10) The (at least) piecewise smoothness of u and A(t, e) allows us to 00 "1 similarly write a(t, e) on [t., t., n ) as a(t, e) = .1 a.(t, e)Ax . We - j o+l - T=p. -l will no longer purport to control o_ such that XAt . | |£.(t., e)|| < e/2. J J *- Instead, we only require that P i, , XAt. o (t., e)Ax J L $ e/2 (ll) '-p. 2 J Once again, in practice, we only have an estimate of a (t., e), but we will assume that the estimate is sufficiently accurate for the discussion to 1/3 2/3p l remain unaltered. Since At. = 0(e ) , it follows from (ll) that Ax = 0(e J ) and hence, 1+2/3P, At. | |a(t, e)| | S e/2 + o(e J ) (12) J We can therefore, still expect the results of the previous analysis to hold, even when the requirements are less stringent. 3. Programming Application A version of the numerical procedure analyzed in the previous section was implemented in a FORTRAN IV program. This section "begins by giving some of the particulars of its operation. For a more detailed account of its functioning, see the block diagram which has been included in the appendix. The program was written specifically to solve the heat equation for various initial conditions. It approximates the second spatial derivatives using centered differences and uses the trapezoidal rule to integrate with respect to time. It was numerically verified that all of the centered difference formulas used give rise to symmetric matrices A( t , e) which are negative J definite. Prior to attempting an integration step, an estimate of the local spatial truncation error is made using centered differences on the current numerical solution to estimate the magnitude of the higher order derivative which appears in the leading term of the Taylor expansion for the error. If the estimated error exceeds e/2 , the current centered difference formula for the second derivatives is replaced by another centered difference formula which is two orders of accuracy in Ax higher. The local spatial truncation error for the new formula is estimated, and if its magnitude still exceeds e/2, the process repeats with the centered difference formula that is higher in order by yet two more powers of Ax. This continues until either a formula is found for which the estimated error is less than e/2 in magnitude or the order of the formula reaches a limit specified in the initialization section of the program. If, on the other hand, the estimated spatial truncation error is much less than e/2, an estimate is made of what the spatial truncation error would 8 be if the current centered difference formula were replaced by a centered difference formula whose accuracy is two orders lower. If the estimated error for this lower order formula still is less than e/2 in magnitude, the lower order formula replaces the current formula. An integration step in t is then taken, and the temporal truncation error is estimated using step-doubling. This procedure involves taking a step of size At and comparing the result with the solution obtained by taking two steps of size At/2. If the estimated error exceeds e/2, At is adjusted downward, and the integration step is reattempted with the adjusted At. If, instead the estimated error is less than e/2, the step is accepted, and At may possibly be increased for the next step. The program was used to solve the heat equation u, = Xu t xx with the diffusivity constant X = 1/10 subject to the initial condition u(x, 0) = f(x) We assumed f(x) was periodic with period 2 and antisymmetric with respect to x = 0, ±1, ±2, ... . Hence, we can restrict our attention to the interval [0, l]. This reduces the run time and also permits us to use centered differences at all mesh points in [0,l], Under these assumptions, with boundary conditions u(0, t) = u(l, t) = 0, the solution is oo 2 2 u(x, t) = Z- b exp(-n tt Xt) sin nnx n=l n where b = 2 / f(x) sin nirx dx n The following initial conditions were used: Case 1. f(x) = sin ttx, b = 1, n = 1, n 0, n > 1. Case 2. f(x) =1 < x < 1. b = U/(m0 , n odd, 0, n even. Here u(x, 0) is a square-wave. Case 3. f(x) - kx, $ x $ .25, 1, .25 S x * .75, h - kx, .75 $ x $ l. b = 8 (sin mr/U + sin 3nir A ) , n odd, n — nu 0, n even. Here u(x, 0) is an isoceles trapezoid on [0, l]. Case U. f(x) = x/a, < x £ a, = 1 - (x - a)/(l - a), a { x $' 1. 2 2 b = 2 sin (mrp)/[n tt p(l -a)]. . Case Ua. p = 1/2 Case Ub. p = 7/10 Case Uc. p = 9/10 Case Ud. p = 999/1000 Here u(x, 0) on [0, l] is a peak with its summit at x = a. As a approaches 1, the slope to the left of p slowly decreases, while 10 the slope to the right of y becomes steeper at an ever increasing rate. For y close to 1, u(x, 0) is very nearly a sawtooth wave. _2 The program was run for each of the above cases with e's of 10 , -3 -8 -9 10 , . . . , 10 except for CASE 1 where it was also run for e = 10 and e = 10 . For all cases except for CASE 1, the integration was begun at t = 1/10 rather than at t = in order to eliminate the difficulties arising from the discontinuities at zero. The initial number of points in the centered difference approximation to the second derivative was three. The initial settings for At and the number of internal mesh points were .01 and 19 respectively. Integration was halted as soon as t exceeded 2. 11 k. Results The results of the tests described above are summarized in the tables which appear in this section. An explanation of each column that appears in the table is given below: e_ - On each integration step, the magnitude of the estimated spatial and temporal truncation errors are both required to be less than or equal to e/2. #_ OATT.fl to RKSTEP - Each time that the subroutine RKSTEP is called, an integration is performed using the trapezoidal rule. Since the program uses step doubling, each time a step is attempted, the number of calls to RKSTEP is increased by three. #_ CALLS times #_ POINTS - Each time RKSTEP is called, the current number of internal mesh points is added to this number. It thus provides a rough measure of the work done by the program. LOG 10 DIFFERENCE in WORK - The difference between the log (base 10) of the entry from the current row of the tf CALLS TIMES # POINTS column and the log (base 10) of the entry from the preceding row in that column. This may give a better idea of the rate at which the work done by this program increases for each factor of ten decrease in e. INTERNAL MESH POINTS - The number of equispaced points {x. } (exclusive of x^ = and x ,_ = l) at which the solution is computed. If the n+1 program has difficulty in taking the first step successfully, it may attempt to raise the number of points . INITIAL COUPLING - The number of points used in the centered difference approximation to the second derivative when the first successful step was taken. 12 FINAL COUPLING - The number of points used in the centered difference approximation to the second derivative when the upper limit of integration was reached. 13 e H W « & O E ^ K M CC < O M K K ft; s « w a -p M < o M C3 O o En W Eh CO g O o s S H H O op* q E s t-3 l*H H H Q CO CO CO Eh H H S J S H W o s Ph i — 1 H o on CM ON t— CO H Pn 13 LTN H ir- CM en on o w o CM on -* -=t -=f 3 Pn H ■« CO CO EH t-3 CO a ON LTN on IT) LTN H _ht 3g H CO VO CO On on on t- O LTN VO o LTN ON On CM O H Ph H CM VO CO H B H LTN Sfc sfc a Ph LTN H -=f .si- On On t- 3g EH -3" LTN CO t- VO CO O CO in en t— C— o £ H * H CM on -5f LTN VO t- u H W H w H H W H H H H H H H CM Ps 15 x < o O H Fn -P < O H E p o o o H h3 o o En en EH H o PL, o ro ro H h-J CQ S sal g H O o H H Ph ^te 3fc CO Pn i-l M 3 o H CO on C\J on o CO o on H LT\ LTN VO t- o I 1 vo 1 VO vo 1 vo 1 on i -=r oo on on I J- J- CM I LTN I W o H CM LTN I w o on i w On o\ on on i w CVJ on CO I vo On J- I LTN H CVJ l/N I W o t- -3- o CM on on CM c— vo O a. on o o o\ H ON ON ON ON CM CM LTN I w LTN I w CM O D— H CO CM H LT\ H LTN LT\ On LTN CO VO rH c— H CM on H H H vo I VO H CM VO I VO H CM H vo o vo On H O o o t— vo on j- H vo on vo ltn VO On ON ON CM J" o h cm on j- j- LT\ H on on H CM H H rH on VO CM vo CM CM -3- ON on j- vo en CO On ON rH CM LT\ CM H ^t LTN VO C— 3 H H H on g CO Eh P CO on EH 16 s -p M <1 W H W O o K CO W g En S EH W O o as H pa o Ph ^ W s W M H Q a CO CO E-t i-5 co a d @ M < U o O M PL, en H OJ OJ no H CO j* LTN t~ LTN t- o on NO NO NO NO oo oo .j- I I i WWW un co j- H H oo I w On H ir\ no I I W W O CO J- CO 00 OO I I w w O OJ -=r cm -3- l w CO 1" w On H ltn I W o M3 I w CO CO LTN CM LTN CM On H O On O t— O H CM O OO ON ON ON -3- CM NO 00 CO O O H 00 LfN 00 t- H On H I CM W oo w -3- LiA w w NO w NO I w ON NO I w On H -3- oo CO o o On ON CM CO NO H H CM CM o-i O ON t— On H ltn 00 -J" _=f t— ON t— LTN t- NO O LTN LTN CO -H- CO H OO o CO H CM NO ON o OO H LTN o OO 00 NO CM NO LfN o H CM LTN CM H I w LTN • II •H IS Ph a3 W o Ph o w CO En 5 CO w Ph- w <; s k S M K fin W s « w 13 -P H <3 1*1 3 h H S O H O En 5 CO H O PL, O S ffi h g o o E s J fe H m a CO 3 s H O o H LL, ^fc 5»fc ro a< h-1 m !5| O EH H CO o IT OO t- H CO O -d- H oo vo l#\ t— t- o 1 CM 1 VO 1 VO 1 VO 1 VO 1 ro oo I I W H t- H CO CM ltn LTN I W ON J" OO VO I W oo I w o oo I w oo CM -3- i w oo CM I ON LTN I w o\ oo vo I w vo o LTN CM CM CM t- O CM O O ON o CM OO VO ON 00 t— ON H H On H ON ON On OO CM On CM VO o CM -=t H CM oo _=r H H H ^ VO I w oo LTN VO I w OO LTN H o c— VO oo C— VO LTN t— 0O H CM CM CM -3" _=T -3" H OO On -=t VO u^ LTN ON CM O ON On LTN OO t- LTN ON o H _=T H H CM 00 vo CM ON \s\ r-\ CO LTN ro -4" VO 00 ON H OO H CM VO OO H vo I~- I I H H II -P Pm O CO &H w Eh 18 H < P2 CM -=)• 00 CO _=r H H H 55 O J- -=r ^t ir\ CO CO K 55 H K H 00 \o vo VO vo hIhK P£ + I 1 1 1 1 hJ Ph ^ H *^! O ^ Ph o -=r 00 1 ^* -H- LfN VO VO >$^ H 1 Ph Ph w 1 pq 1 Ph 1 S 5? Ph O o CM t— H -=f o S M Ph VD oo O oo O CM oo ph W on H V£> rH OO vo H oo 00 J" J" LfN vo vo K Ph I I 1 1 1 I 1 38 £ H w W Ph Ph W H LTN OJ CM t— H -=t- o S K W C\J CM O OO O CM 00 W t— CM VO H 00 VX) H 3 VO ^H- C— vo O H CM O t- CO o 55 -P LTN -h" CM ON J" <-\ H H < CO CM CM o o o o P>H O S 33 on LTN LTN t- t— t- t- Hh£ h o o o %n H J LTN t- ON H 3 H H ^£ H H H s o M O 3 g w e Ph CO 55 ON On ON O LTN CM O w g H H H H CM CM 00 _=* Eh S 55 £ M 8 o s £ H H O VO CO O VO 00 OO Ph Es LTN 00 00 H _3" CM O Ph O CM oo -H- _=r J- J Ph R M fi CO CO EH >-3 CO 55 H c— 00 On J" LTN ON H LTN CM CO H J- LTN c— 5 < Q < a k S M K Pn W 3 a -P M < o a M § D o u 19 CO ro rH t- rH ON lo .* OJ On LTN t- rH oo C— ion NO NO i w OO NO OO CO w LTN OJ i w 00 NO US J- LA w W 00 CO o t— H C\J no p4 H OJ no OO I _H/ On NO ro i w H OJ OJ _H/ I 00 no LA I w o H LA I w OO t- OJ NO I w H OJ no o LA OJ OJ OO CO CO o _H/ o OJ o en NO w o CO NO I Pm o ro On o On o o t— o ^n M 1-3 LA t- H H H rH H EH Ph r-l H H rH H H P a O H O 3 CO a w EH K CO a On ON On OJ CO NO LTN W g M rH H H OJ OJ OO _=r s H g o a pq H w o OO NO OJ -=* LTN NO pq ^ LA On o OJ ro OJ CO W O P«H O H -d- -=f -=r _=r a j E H M Q J CO CO EH a ON LTN LT\ c— H OO ^t 53 H H CO NO ^* OO On NO OJ o LTN NO o NO On o CO O M ft H OJ NO On o EH H LTN =Jfc =»fc CO Ph t-q H US H H t- OJ OJ OJ 3£ EH J" LTN CO t- t— On Os CO H OO t- NO O K rH =tfc H OJ OO -=r us NO t— (J W w W H pq W w H H rH H H H H CO On On On -P •H o p^ o (in CO EH CO 3 20 FINAL At - The stepsize when the upper limit of integration was reached. MAX ERROR EVER - The magnitude of the worst error in any individual component of the numerical solution vector that was ever found. The numerical solution was compared to the exact solution on each step that t equalled or exceeded an integral multiple of 2/10. MAX FINAL ERROR - The magnitude of the worst error among the components of the final numerical solution vector. LOG 10 DIFFERENCE IN MAX FINAL ERROR - This column is to MAX FINAL ERROR as LOG 10 DIFFERENCE IN WORK is to § CALLS TIMES # POINTS. Results indicate that the idea of controlling the global error through separate control of the local spatial and temporal discretization errors is indeed quite feasible. Of particular interest is the LOG 10 DIFFERENCE IN MAX FINAL ERROR columns. These show that for each factor of ten decrease in e, 2/3 the global error decreased by a factor which appeared to approach (10) This is as predicted in the previous section. The LOG 10 DIFFERENCE IN WORK columns show that as e was successively decreased by factors of ten, the amount of work done increased by factors 1/3 1/2 which tended to stabilize in the range of (10) to (10) with the exact number depending on the initial condition. For CASE 1, where it was never found necessary to increase the number of internal mesh points, the number was 1/3 precisely (10) . This is what would be expected with the trapezoidal rule. Since the local truncation error is proportional to the cube of the stepsize, and the program attempts to keep the local errors on the order of e, if the -1/3 mesh is fixed, the number of integration steps performed should be 0(e ). Combining the relationships between e and the global error and between e and the amount of work done by the program, we can get a relationship 21 between the global error and the work. Denoting the global error by e and the work by q, we can write 2/3 limit e(e) = c e E -*■ limit q(e) = c e~ m 1/3 ^ m $ 1/2 e -*■ and hence, limit e(q) = c q P -2 £ p $ -k/3 e -> J where c. , c_ , and c are constants. While this may not seem like a particularly good return of accuracy for the computational effort expended, it compares quite favorably with most of the methods commonly used to solve partial differential equations. The well-known Crank-Nicolson method which for the heat equation gives rise to the system n+1 n . ,_. , / .2 n „2 n+l N u. = u. + 1/2 Ar( 6 u. + 6 u. ) i i x l x i where u. = u(x. , t), r = At/Ax , and 6 is the standard 3-point centered 11 x 3 2 difference, has a local truncation error of 0(At + At Ax ) . The Douglas formula is the most accurate formula involving the same six points to approx- imate u (x., t ._). It generates the system xx i n+1 n+1 _ n , , x 2 n , x 2 n+lx 1 ,,.2 n x 2 n+l x u. = u. + 1/2 Xr[& u. + 6 u. ) + r-r X( 6 u. - 6 u. ) i i xi x l 12 xi xi for the heat equation which results in a local error- of 0(At + At Ax ) . If, as At approaches zero, r is kept fixed, the local error is 0(At ) . The global 2 error should therefore be 0(At ). This means that limit e(At) = 0(At 2 ) At -> just as for the trapezoidal rule. In advancing the solution one time step, Douglas' method must solve a 22 tridiagonal system. The amount of work necessary to solve such a system is proportional to the number of equations in the system and hence to the number of internal mesh points. But since the number of internal mesh points is always within one of Ax , we may take the work proportional to Ax . The number of steps that must be taken and hence, the number of times the system must be solved is proportional to At . Hence, for Douglas' method, Ax At may be taken as a measure of work analogous to the measure # CALLS TIMES # POINTS described above for the program. If, as before, r is assumed to remain fixed as At varies, for the Douglas method we have limit q(At) = 0(Ax _1 At" 1 ) = 0(At~ 3 ' 2 ) At ■*■ and hence, limit e(q) = 0(q /3 ) At + which is no better than the worst case observed for the program. 23 These include X, e, At . , At mm max H me:.h points max and coupling c Sturt 3 Read in various parameters ffor programs Coupling * Coupling + S T.etting starting values, etc. Various .initializations) Coupling is tne number of points in the difference approximation to the second spatial derivative Coupling * '', Uart « n Xo is the magnitude of the largest component in the estimated spatial truncation error Calculate Start Set up A matrix for current coupling 'Calculate LU ' decompositions for I T x ' and T + fit , £ C Take step of size At and store result in varray YWH0LE y Stop 3 V 2U n mesh points * Biech points Calculate HSAT, the smallest num- ber of internal mesh points r. ,t. Xito < e/2 mesh points NSAT Subscripts max or "min" indicate maximum or minimum per- missible values <3 Take 2 steps of size At/2 and store final | result in array YNEW Calculate t(At) = 1/ 3*Max | YNEW( i ) -r.mOLEC i ) i T * T + it YES ,/Print message^ \ and STOP J SAT Set it SAT Max(4t SAT' A W At/5) 25 Print computed solution in YNEW, the exact solution, t and the errors printout time * printout time + .2 YOLD ■*■ YNEW Calculate FACTO R#l= factor by which At could be in- creased and rse/2 •still hold I FACTO PJl «- .8 * FACTO R#l FACTO R#l FACTO R#l FACT0PJ1 Min(FACTOR#l, 5) Max(FACTOH#l, l) Min( FACTORS!, At max /At) Calculate At \o \ F FACTOR #2 «- e/(2 At Xt) 26 Calculate At Xa for coupling - 2 Coupling ■*- coupling - 2 FACTOR «- Min( FACTO R#l, FACT0R#2 * .8) ir At + At y FACTOR \ 1 c FormAEC-427 U.S. ATOMIC ENERGY COMMISSION < 6/68 > UNIVERSITY-TYPE CONTRACTOR'S RECOMMENDATION FOR AECM3201 DISPOSITION OF SCIENTIF:C AND TECHNICAL DOCUMENT ( See Instructions on Reverm Side ) 1. AEC REPORT NO. COO-2383-0025 2. TITLE AUTOMATIC INTEGRATION OF THE HEAT EQUATION 3. TYPE OF DOCUMENT (Check one): tf l a. Scientific and technical report ] b. Conference paper not to be published in a journal: Title of conference Date of conference _^_ Exact location of conference _ Sponsoring organization □ c. Other (Specify) 4. RECOMMENDED ANNOUNCEMENT AND DISTRIBUTION (Check one): ET) a. AEC's normal announcement and distribution procedures may be followed. ~2 b. Make available only within AEC and to AEC contractors and other U.S. Government agencies and their contractors. 3 c. Make no announcement or distrubution. 5. REASON FOR RECOMMENDED RESTRICTIONS: 6. SUBMITTED BY: NAME AND POSITION (Please print or type) C. W. Gear Professor and Principal Investigator Organization University of Illinois Department of Computer Science Urbana, IL 6l801 Signature ^i\j2l t0 Cft-r Date October 3, 1975 FOR AEC USE ONLY 7. AEC CONTRACT ADMINISTRATOR'S COMMENTS, IF ANY, ON ABOVE ANNOUNCEMENT AND DISTRIBUTION RECOMMENDATION: 8. PATENT CLEARANCE: I I a. AEC patent clearance has been granted by responsible AEC patent group. I I b. Report has been sent to responsible AEC patent group for clearance. I I c. Patent clearance not required. BIBLIOGRAPHIC DATA SHEET 1. Report No. UIUCDCS-R-75-75 1 * 3. Ret ipienl '■• Ai i i isii in No. 4. Tn It- and Subt it 1 1- AUTOMATIC INTEGRATION OF THE HEAT EQUATION 5. Report i). iti- October 1975 7. Autrtor(s) Michael H. Ostrar 8. Performing Organization Rept. No " UIUCDCS-R-75-75 1 * 9. Performing Organization Name and Address Department of Computer Science University of Illinois at Urbana-Champaign Urbana, IL 6l801 10. Project/Task/Work Unit Nc 11. Contract /Grant No. US AEC AT (11-1) 2 383 12. Sponsoring Organization Name and Address US AEC Chicago Operations Office 9800 South Cass Avenue Argonne, IL 60^39 13. Type of Keport & Period Covered 14. 15. Supplementary Notes 16. Abstracts A scheme for automatically integrating the heat equation is discussed. A program based on this scheme is explained, and the results of integrating the heat equation for a number of different initial conditions are presented and interpreted. 17. Key Words and Document Analysis. 17a. Descriptors 17b. Identifiers /Open-Ended Terms 17c. COSATI Field/Group 18. Availability Statement Unlimited 19. Security Class (This Report) UNCLASSIFIED 20. Security Class (Thi Page UNCLASSIFIED 21. No. of Pages 33 22. Price FORM NTIS-35 ( 10-701 USCOMM-DC 40329-P7 1 *«« % 1* ^ i& wMm^MMmuuDtiuiuumi