LIBRARY OF THE UNIVERSITY OF ILLINOIS AT URBANA-CHAMPAIGN 510.84 IJifcr T\p.54l-M6 Digitized by the Internet Archive in 2013 http://archive.org/details/automaticsolutio546lars 'JlLflS UIUCDCS-R-72-5^6 ~p/UZtf COO-1U69-0212 1 October 1972 AUTOMATIC SOLUTIONS OF PARTIAL DIFFERENTIAL EQUATIONS Leonard Andrew Lars en THE LIBRARY op XHe DEC I- 1972 UIUCDCS-R-72-5U6 AUTOMATIC SOLUTIONS OF PARTIAL DIFFERENTIAL EQUATIONS* by Leonard Andrew Lars en October 1972 DEPARTMENT OF COMPUTER SCIENCE UNIVERSITY OF ILLINOIS AT URBANA- CHAMPAIGN URBANA, ILLINOIS 6l801 Supported in part by the Department of Computer Science and in part by the Atomic Energy Commission under contract US AEC AT(ll-l)lU69 and submitted in partial fulfillment of the requirements of the Graduate College for the degree of Doctor of Philosophy in Computer Science, October 1972. Ill ACKNOWLEDGMENT I am indebted to Professor C, W. Gear for advising this thesis; I appreciate his suggestions and guidance. I am also grateful to the Atomic Energy Commission for supplying the computer time needed to complete this research. TABLE OF CONTENTS Page 1. INTRODUCTION , 1 2. PROCEDURES 1+ 2.1 Spatial Discretization 1+ 2.2 Error Estimates 18 2.3 Efficiency 30 2.k Initial Partition 35 2.5 Changes in Integration 38 2.6 The User's Equations ko 3. PROGRAMS 1*3 3.1 PARDIF k3 3.2 XSTP k6 3.3 DERIV U8 3.U DIFSUB 52 3.5 DEC0MP and S0LVE 5U 3.6 DIFFUN 5U 3.7 DIVERR 55 k. RESULTS 56 LIST OF REFERENCES fk APPENDICES A. THE SYSTEM OF PROGRAMS 75 B. EXAMPLE 1 119 C. EXAMPLE 2 122 D. EXAMPLE 3 125 VITA 129 1 . INTRODUCTION Define an automatic program as. one that can 1. take as input members of some broad class of problems, 2. try to control errors within specified limits, 3. try to minimize the work required for the solution. The problem of the present paper is to look for an automatic method that can be applied to some subset of partial differential equations. The primary value of a program involving an automatic method would be to provide the person who needs only a few runs with a particular type of equation the opportunity to find a solution, within acceptable tolerances, without having to write and debug a specialized program. In this environ- ment, the program will have no information about the solution, or how errors might accumulate. While maintaining sufficiently small errors, an attempt needs to be made to keep the cost of such a solution as low as possible. To accomplish these goals, decisions need to be made during the integration for balancing the two conditions above based on some sort of error estimate and the user's error tolerance. In the present study, the subclass of partial differential equations was chosen to be the set of second order mixed initial boundary value problems in two variables (x and t) of the form F(u, u , u , x, t) = u . The general X XX T> method used involves integrating along lines which are parallel to the time axis and which pass through a finite set of points on the spatial axis. To obtain an idea of this method, consider the equation 9u _ afu 9t ~ 2 dX with initial condition u(x,0) = sin(Trx) x in [0,l] and boundary conditions u(0,t) =u(l,t) = To begin, select a partition P = {0 = x Q , x^ x 2 ,.. ., x n _ 1 , x n = 1} and then define a set of functions u. (t) such that u.(0) = u(x. ,0) i = 1, 2,..., n. By the use of some interpolating polynomial, a discrete analog is obtained for the spatial derivative, which will be represented by in 8x 2 X. 1 2 W Using this analog, we obtain the discrete system du. -~ = A (u. ) i = 1,. .. , n dt n i This system differs from the system 9u(x,t) 3t 3 u(x,t) 3x x. which is obtained from the true solutions by an amount 6(i,t). The analog A (u. ) and the partition P are chosen in an attempt n i to keep u. (t) - u(x. , t) within the user's desired tolerance. Once the partition and analog are chosen, one then integrates the system in the time direction along the lines implied by P from t = t to t = t If at any time during this integration it appears that the error is too large, the step is repeated with the possibility of changing either A Q (u i ) or P. A change in the partition would mean continuing the integra- tion along a new set of lines. There is also the possibility that the error is much smaller than is required. Again, modifications in A^iu) or P might be attempted so that the work of obtaining the desired error is minimized. These decisions will need to be made on the basis of some estimate of the error. The basis for the integration of the ordinary differential equation system above is a program (DIFSUB) developed by Gear [l] for automatically solving systems of stiff ordinary differential equations. As an automatic program, it attempts to keep the errors introduced in integrating in the t direction within a given tolerance . In integrating along the possibly changing set of lines, there are thus three ways of trying to keep the errors in the solution at a desired value. First, the partition can be changed. Second, the nature of the analog used to discretize the system can be altered. And third, the amount of the error parameter given to the integration routine can be modified. In order to achieve an automatic program for the entire system, several items will need to be considered: 1. the forms of the analogs to be used, 2. the methods of obtaining error estimates, 3. the methods of obtaining and modifying partitions, h. the implementation of the user's equations in an easily usable form. 2 . PROCEDURES 2.1 Spatial Discretization There are several problems to be considered in looking for a spatial analog. These include the number of equations required, the ease of changing the partition, and the situation at the boundaries. Three methods were initially considered: 1. Lagrangian methods, 2. Hermite type methods, 3. Linear transformations of the Lagrangian methods. These will then be used on equations of the form — = f(u, u , u , x, t) 3t x xx The equations are then put into a discrete form by replacing u by a set of functions u. (t) and replacing the right side by an analog A (u. ) . This l n i then gives the system of equations du. ar=VV w In the case of the Lagrangian polynomials, there are several possibilities. First, there is the class of polynomials using equidistant points. If the number of points (n) in the interpolating polynomial is odd, then it is possible to take advantage of symmetry properties to eliminate one of the terms in the error. Given the approximation n y(x) = £ I (x) f(a ) (2) i=l to a function f(x) for interpolating points a.,..., a , the error term is given by -Sr- *m where p^x) = (x-^ ) . . . ( x-a Q ) and n is some point in the interval containing these points. Taking the derivative of the interpolating polynomial, an approximation is obtained for f'(.x). n y'Cx) = z i ;Cx) f( a .) to) 1=1 x x To obtain ideas about the error, take the derivative of the error term of (2) (see 2.2 for discussion and development). For equation (3), the error term would be given by n! f (ri l )+ InTnT f V (4) Since p n (x) = at tabular points, where all such evaluations will be made, (k) reduces to _n An), v ^T^ f (ri i } (5) If this process is repeated, the error for the second derivative is ^ f(n » K)+ ^ f(n+1)(n2) + ^ f(n+2)(n3) (g) If the interpolating polynomial uses an odd number of points, and if the derivative is estimated for the center point, then the first term is zero. Consider n n Pn (x > P"(x) = I I , "L i=l j=l (x-a^tx-a ) Since evaluation is at the center point, x - -^ = , and all terms containing this factor will vanish. As a result, P ;'(x) can be rewritten as (where the subscript nl is 2^) 2 ' n 2-p U) iftexl Further, since the points are symmetrically located around a CT) nl (a -a +i)=(a n -a,-i). Noting this symmetry, (T) can nl nl nl nl he written hi 2-p (x) 2 'pJ x ) P n U) = !=1 '^V^V^TT + U-a nl )U-a nl = D = Since the last term of (6) is also zero, the expression reduces to 2p'(x) , ._ v ^TT^^In) (8) Since p (x) is a function of the spacing h, the order could he defined as the lowest power of h contained in the error. In (5) and (8) these orders are the same, so that the truncation errors should react similarly to changes of h. In the event the symmetry is lost, the order associated with the second derivative is reduced hy one, and P^(x) is no longer zero. The error term for the second derivative would prohably then he the determining factor in selection of an h. One of the main difficulties with this class of methods results from the equidistant spacing. If the solution is "non -smooth" , the spacing for the partition will have to he set so as to reduce the error at the worst point. This could lead to a much larger number of equations than needed, for in the smoother regions the error would he expected to be much less than is required. Ideally, one is aiming to have partition points 7 such that the total error introduced at each is at the level desired by the user. In an attempt to keep the error closer to the user's request, one could divide the interval into sub intervals with spacings that are locally equidistant. If the analog uses five or more points, such subdivisions have problems near the end of each subinterval. The derivative values near these subinterval end points would either have to use a nonsymmetric form or a point from an adjacent interval with possibly a different step size. When using a three-point analog, these problems are reduced. Each analog could use spacings that are multiples of some h. In order to circumvent the problem of the end point, one could consider it as the center point of another subinterval overlapping those on either side. For example, given the partition p 6 P3PUP5 and an h ' 15 3 71 ITIS8I6T rr, the analogs could be set up in the following way derivative at points used v v y 2 , p 6 P 2> p 3' V h P 3 » P^, P 5 v k , p 5 , p 6 P r P 6 , P ? 8 In practice, enough points would be added so that several of the sub- intervals would not be contained within another interval. Another way to try to keep all errors somewhat close to the user's request would be to use unequal spacings throughout the interval. There are several disadvantages to this approach. First, the order of the error term would be lower than those associated with symmetric- equispaced ones. Second, there exists the problem of calculating the derivatives. With equidistant points, a set dependence exists between values of the function and the derivative analog at any partition point . In the case of unequal spacings, either the analog would have to be calculated each time it was. needed or the dependence of A (u. ) on each n 1 point would have to be calculated and then stored for a given partition. Third, the initial partition would be somewhat more difficult to develop because of the freedom of changing the spacings between the points. As a result, it could happen that the spacings between any two pairs of adjacent points could be greatly different. In this case, it would be desirable to go back and change other points so as to avoid widely differing spacings for neighboring points. A second class of analogs involves those methods that use more than just the functional value at each point. The primary subclass would be the Hermite interpolating polynomials. With polynomials of this form, it is necessary to save the first derivative at each point as well as the functional value. In order to develop this additional value at each time step, it would be necessary to solve an additional equation at each point. Such an equation can be found by taking the partial derivative with respect to the space variable of both sides of the given equation U t = f(u ' V U xx' x ' t) (9) This gives 3l U t = 3x" f(u ' V u xx' X ' t} ("J Since the solution u is assumed to be continuous and to have as many- continuous partials as may be required, the order of differentiation is changed to give 3t U x = ^ f(u ' V U xx> x > t} (id At this point, u is replaced by a new variable v. and the x l right side of (11) is replaced by a discrete analog B (u. , v.), giving n' i' i 3 3? V i = VV v i } (12) The solution v. of this addition equation is not expected to be ill-behaved provided the solution is continuous. Consider the Taylor series expansion of u u(x,t) = u(x ,t Q ) + u x (v ,t )(x-x ) + u t (x ,t )(t-t Q ) + .,.. The presence of u x in this expansion indicates that any unusual behavior that might be reflected in u x would be present in u. As a result, if the original system of equations (9) can be solved, then it is expected that the combined system of (9) and (ll) can also be solved. Just as in the Lagrangian case, there is an advantage in using the symmetry properties to eliminate terms in the expression for the error. Further, all of the problems associated with the Lagrangian polynomials — when more than three points are involved— are compounded by the fact that the first derivative must be calculated and saved for each point. A new problem exists in that additional values are needed at the boundary. In 10 the case of the original equations, the boundary conditions provide additional equations that can he solved with the differential equations . With Vn and v at the boundaries, it will be necessary to develop other equations, or at least provide values, for these variables. Another possible member of this second class of interpolating polynomials would be those that match first and second derivatives as well as the functional values at the interpolating points . The use of such polynomials (trimite) and the problems involved will be very similar to those encountered in the Hermite case. Here, two additional equations will need to be supplied for each, point. The first of these would be the same as the additional equation in the Hermite case. The second would be found by differentiating both sides of (ll) with respect to x giving 2 3 9 3 „/ . \ T" TT u — T f(u, u , x , x t) 3x 3t x _ 2 x xx , 3x ' dv. Interchanging the order again, introducing an additional variable w. =— — l dx and an analog C (u. , v. , w. ) , the discrete version becomes n l l l dw. jr =c n ( v v w i> (13 > Another problem, that was also present in the Hermite but becomes more evident in the trimite polynomials, is the presence of higher ordered derivatives. The error term for the trimite will contain the ninth and tenth derivatives. Although accuracy will not be required, a reasonable approximation is desirable. The one additional possibility that was considered for obtaining an analog was to use a linear transformation of a Lagrangian method which would save different quantities, but which would possibly save work. Suppose 11 eight interior points and two boundary points (0 = x ,..., x = 1} are used with the equation u t = u^ x e[0,l]. The discretized equations du would be represented by _£ = Au where A is given by Q Q 1 -2 1 Q 1 -2 1 1 -2 1 1 -2 1 1 -2 1 Q 1 -2 1 1 -2 1 1 -2 1 n, it would be easier t< 3 save l u_ . l] . l i i l" » 1' l' "3 s "3' "5' u 5 ' u T , u^ instead of the functional values at the eight interior points. The transformation T would then be given by 1 0. 1 1 -2 1 1 -2 1 1 -2 1 1 1 -2 1 1 1 -2 1 1 12 -1 The transformed matrix B = TAT " then become 1 2 -3 l a l a 1 -2 Q 1 1 1 -2 1 a 1 1 -2 1 0. J 2 For 3u = l_u the trans formation gives results that would perhaps he 3t 3x worth using. However, if the original equation contained the first derivative as well as the second, a second matrix D = TCT -1 needs to he calculated where 2'C is 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 L' 13 and 2*D is i 2 -2 -1 -k 5 2 -1 -2 4 2 -2 -1 k -7 -1+ 1+ 2 -1 2 -It _P k 2 -2 -1 -1+ 8 It -1 -It It 2 -1 -2 U 2 -k -2 1+ 2 -2 -1 1+ -8 J* 8 1+ -7 -1+ It 2 -1 Since the user's equation is unknown, the need for both matrices is extremely likely. Therefore, the idea of saving the second derivative loses its appeal, If the function and the first derivative are saved, then the results are similar. For the first derivative, the resulting matrix is (B = TAT" 1 ) . 1 1 "It 1 It 1 B = 1 IT 1 2 1 IT 1 1 IT 1 ~2 1 1 1 it 1 2 1 IT lU ,-!• and that for the second derivative is (D = TCT ) D = G a 2 -2 -2 Q Q 1 1 2 2 1 2 2 -k -2 -2 1 2 1 "2 2 -k 1 2 -u -2 -2 -2 1 2 2 -i+ a -k 1 2 -u -2 -2 -2 1 2 Since the problems the procedures are meant to solve could involve both derivatives, the -whole idea of using linear transformations does not appear to provide anything of interest. As a result of the disadvantages of going to more than three points it was decided that three methods would be used. These were the Lagrangian, Hermite, and trimite polynomials using three points. This would provide methods of three different orders. Although the Hermite and trimite forms have additional problems, the flexibility of having the different ordered methods and the possibility of a fewer number of equations makes consideratic of these methods worth undertaking. If the polynomials for Lagrangian, Hermite, and trimite are called Y , Y , Y , respectively, the resulting formulas can be written L H T 15 n Y = I Z.(x) f(x.) L i=l X X n Y H = = L {h (x) f(x. ) + h,(x) f'( x .) i-1 -L 1 1 1 n T - 2 (t (x) f(x ) + t (x) f»(x.) + t.(x) f"(x.)) i=l j. x ii x where t.(x), £(x), and t ± (x) are polynomials associated with the trimite method. The development of t, t, t parallels the development of the Hermite polynomials. Let \ = (a ik ^ + a 2k X + V *k (x) \ = (b lk ^ + a 2k x + a 3 k } *k (x) \ = (c ik ^ + C 2k X + C 3k } £ k (x) There are three conditions for each of the equations ? k ( Xj ) = o t.( Xj ) = o !»(«, > . « k . Solving each row gives 16 b lk ■ -*£ A n.aCu t ), x. , t) ( 15) The analogs A^tu,) and A^^fo.) produce a truncation error d x ( Xl ,t). An estimate of d^x^tl can be obtained by perturbing the right hand side of equation (.15) by the error estimates of A (u ) and n,x i A n,xx (u i ) (e i and e 2' res -P e c-tively). This would give the expression f( v A n,x ( V + e r V (u i } + v x i> *> Subtracting the expression (lU) from this, the result d x (x is t) = f(u., A njx (u.) s A njxx (u.) 5 x., t) - f(u., A n)X (u.) + e r A^Cu.) + e 2 > x i' ^ (16) is obtained. In order to be able to obtain this estimate, the values of e and e 2 will have to be obtained. If the analogs involved were just those of the Lagrangian polynomials, then several approaches are available for obtaining the error terms e ± and e g . Since the present approach is to use Hermite type polynomials as well, the approach of Wendroff [2] will be followed. If y(x) is any polynomial approximation to f (x), then f(x) = y(x) + p( x ) G(x) (17) where P(x). G(x) is the error term, G(x) is an unknown function, m. m m P(x) = (x-a^ X (x-a 2 ) 2 ...(x-a k )^ : and m is the number of values obtained at a . x i In order to find the error estimate of f^'(x) (the trimite method would require that j=U), take the j-th derivative of both sides of (17) 20 For the first four derivatives (assuming that G is continuous and has j continuous derivatives) the results obtained would he f'(x) = y'(x) + P'(x) G(x) + P(x) G'(x) f"(x) = y"(x) + P"(x) G(x) + 2P'(x) G'(x) + P(x) G"(x) (iQ) f*(x) = y"'(x) + P m (x) G(x) + 3P"(x) G'(x) + 3P'(x) G"(x) + P(x) G"'(x) fH(x) = ymt( x ) + j*( x ) G(x) + UF»(x) G«(x) + 6P"(x) G»(x) + UP'(x) G"'(x) + P(x) G""(x) The first problem is to determine the validity of the assumption that G G' G", G ,M , G m ' are actually continuous functions. To do this, for G, solve (IT) for G giving G(x) = f(x) - yU) (19) Glxi P(x) Since f(x) is assumed continuous, G(x) is continuous whenever P(x) t 0. Consider now a point a. , p (<3) (a.) = for <_ J < m. but P x (a.) t 0. It is also noted that [f(x) - y(x)] = for <_ J < m.. Therefore, G^) is defined by m. m. f X (a,) - y 1 (a.) G( a. ) = m P X (a.) which by L'Hospital's rule shows the continuity of G. To show that G' (x) is continuous, Wendroff [2] takes the derivative of both sides of (19) giving ,, v P(xUf'(x) -t'(x)1 -P'( xUf(x) -y(x)] (20) [pun 2 He then defines m, h(x)- U^flxf^ (x-a ± ) E^ 21 where k is the number of interpolating points. Substituting this into (20) and using L' Hospital's rule again, he arrives at the continuity of G'(x) (m i +l) provided f (x) is continuous for the maximum m . i The same approach can be used to show that each of G" , G m , G"" is also continuous provided f has a sufficient number of continuous derivatives. Since the only major difference is that in each succeeding case the computations become much more involved, only the proof of the continuity of G" will be shown. To this end, take the derivative of both sides of (20) giving _ 2 (gf) 2 (f(x)-y(x))-|^f (f (x) - y(x) ) u W -^ 2( PTxT' ) < f, ( x )-y'(*)) + (F"(x)-y" (x)) P(x) (21) Again, the problem of continuity exists only at points where P(x) = 0, since otherwise all functions are continuous. So consider a point a i where P(x) has m i zeros. Multiply the numerators and denominators of (21) o by (x-a^ giving (let f e f( x ), y E y( x ), P E p( x )) _ ^*~\) 2 (f-) 2 (f-y)- (x-a,) 2 p(f- y )-( x -a.) 2 . 2 .£l (f'-y>) (x-a.) 2 P (x-a ) 2 (f"-y") + =- (x-a.) 2 P (22) Since the denominator has m. +2 zeros at a. and the (m. +2)-nd derivative is nonzero, the limit of G'^a..) exists— provided the numerator 22 vanishes for the first m.+l derivatives. The continuity of G'(a i ) and definition of the derivative then would allow the conclusion that G"(a i ) is continuous. It is assumed that f(x) has a sufficient number of continuous derivatives . Define pi k m h(x) = (x-a.)-= (x-a.) Z ^j and d(x) = (x-a ) 2 f-= (x-a ) 2 [ Z ( ? ) + Z J \ 2 1 (23) 1 P X j,£=l U a j n V J*l (x-a.T It is noted that h'(x) and d'(x) are then given by k m. (x-a. )m. j=l (x - a o } (x-a/ and 2 k 2 (x-a.) m.m, 2 (x-a. ) m.m d'(x) = Z [-, 1, i \ - i$ ^-H j^=l ^V^V (x-a,) 2 (x-a) k m.m. (x-a. ) m.m. j=l ^^V (x-a.) 2 o k 2(x-a. ) m. (m.-l) (x-a.) m. (m -1) + z [ 1 ,1 ,1 _ 1 J— J ] (2U) J-l (X " a J } (x-a.) 2 If (23) and (2k) are evaluated at a. , they become 23 h(a.) = nu, d(a i ) = m. (m. -l), km km h'(a ) = Z t — J-^r, d'(aj = 2m. Z 7 — J-_ Substituting these into (22), the numerator G"(x) becomes 2h (x) (f-y) - d(x) (f-y) - 2{x-e, ± ) h(x) (f'-y') + (x-a. } 2 (f"-y") (25) The derivatives of each of the four terms (call them A, B, C, and D) in (25) are now evaluated using Leibnitz' rule. This gives 1. A = 2^-(h 2 (x)(f-y)) = Z 0(h 2 (x)) (r) (f-y) (a - r) dx r=0 This term is now evaluated at a . Since (f-y) has m. zeros, A is identically zero for < a < m — 1 (m. ) 9 (m. ) for a = m . , A = 2«ti(x)(f-y) X = 2m? (f-y) x for a = m i+l» 2 (m.+l) (m. ) A = 2-h (x)(f-y) + Uh(x)-h'(x)(m.+l )(f-y) X 2 (m.+l) fo) k m. = 2m? (f-y) x + Urn. (m.+l )-(f-y) X Z — J- 1 11 . , x-a . j=l J 2. B 5 -^_ (d(x)(f-y)) = I ( a )(d(x) (r) (f-y) (a - r) dx r=0 r B is now evaluated at a . For 0_1, and if f(x) has a zero of at least m at a , where a < a, : a <...< a, < b, and if x 1 — l— 2— — k— ' k 2 m. ^n+1, then there exists an n e[a,b] such that r n ^(n) = 0. 26 In the case at hand, H(z) has at least n+5 zeros including five at z = x, and n associated with P(z). This would mean that there is an n such that H^ (n) = 0. But on taking the derivatives, after multiplying out the right side of (26) and noting that P(z) G(z) = f(z)-y(z), the function "becomes H (n + M (n) =f (*+D (n) + in±iil G m< x ) = Solving this expression for G""(x) gives The general form would then "be o (J) (x)= I ^ TT f (n+J) (n) In order to know the value of the error estimates given in (l8), the values P, P' , P", P"', P""need also to be obtained. Since the polynomials in the present case are the Lagrangian, Hermite and trimite polynomials, P(x) takes on the form p(x), p 2 (x), p (x) where p(x) = (x-a^ (x-a 2 ) . . .(x-a n ) . The derivatives of these polynomials are given in Table 1. 27 Lagrangian Hermite P 'W P*(x) 2p(x)p»(x) P "( x ) P n (x) 2[(p'(x)) 2 +p(x)p"(x)] F"(x) not needed 2[3p' (x)p"(x)+p(x)p"'(x) ] P"'(x) not needed not needed Trimite p '(x) 3p 2 (x)p'(x) P"U) 6p(x)( P '(x)) 2 +3p 2 (x)p"(x) F"(x) 6(p'(x) 3 4-i8 P ( x )p'(x)p"(x)+3p 2 (x)p"(x)) P"^) 36(p'(x) 2 p"(x)+l8p(x)(p» , (x)) 2 +2^p(x). P'(x)p , "(x)+3(p(x) 2 p"'(x) Table 1. Expressions for Error Terms Considering that p(x) = and p"(x) = 0, since the calculations are done at the center point of an odd number of interpolating points (in this case n=3), the substitution of the values of Table 1 into (l8) yields 1. for Lagrangian polynomials first derivative error estimate is P 3 U) f <3), , , , "51 — f (n) (27a) second derivative error estimate is 2P 3 U) Jk). , , , —IT! — f (n) < 2 Tb) 28 2. for Hermite polynomials second derivative error estimate is 2(.pi(x)) 2 ,„ — ^ f (6, (n) (28a) third derivative error estimate is 6CpUx)) 2 m ~ f (n) (28b) 3. for trimite polynomials third derivative error estimate is «'^» 3 „<9> f^'(n) (29a) fourth derivative error estimate is f (l0) (n) (29b) 2U ^3 (x))3 .(10) 10! The values of these estimates are then put into equations like (l6) in place of e and e to determine d (x. ,t). In the case of the Hermite or trimite, d (x ,t) would be found by using equations (12) and (13) for XI the development of d (x. ,t). In the time direction, d. (x. ,t) can be found from the automatic program used to solve the system of equations generated by the discretization, Once these error estimates have been determined, there are two ways that the overall error could be approached. First, the sum of the two could be used by asking that during any time step this sum would be less than the user's requested error. This would be saying to the user that the equation being solved is close to his. However, the actual error in the solution will depend on the nature of the equation. If, therefore, the 29 user desired a particular error in the solution, he would need to have some idea about the nature of the solution so that a value of the error request can be obtained. This procedure would take less computer time but would involve more thought for the user. The second approach would be to attempt to obtain an estimate of the error in the solution. This could be done by numerically solving a system of asymptotic error equations along with the discrete system of differential equations. The system of equations being solved is of the form du. The system differs from the desired system in that the errors d and d x have been introduced. This would lead to a system 3u(x. ,t) ~ =f(u(x.,t), A nfX (u(x.,t)),A njXx (u(x.,t)), x., t) + d t + d x Subtracting and letting e.(t) = u.(t) - u(x. ,t), the equations obtained are de. dT^V A n,x (u i } ' A n ,xx (u i } ' V V e i " d t - d (30) Since the solution to this equation need not be accurate, it may be possible to use a lower order method to solve this last system provided f (u(t)) is li- the same for both systems. The advantages of the latter method led to solving the system (30) along with the original system. 30 2.3 Efficiency Once methods are available for solving the problem and for estimating the errors, attention needs to be given to producing an efficient system. That is, the solution needs to be found in such a way that it has certain error properties and such that the work involved in obtaining that solution satisfies some sort of minimization of computer time. If greater accuracy is desired, it is expected that more computer time will be used. If there is a decrease in computer time, then provided the overhead is always kept at a minimum, a loss of accuracy is expected. The goal would be to have the error in the result less than, but as close as possible to, the amount requested by the user in such a way that no more computer time is used than is absolutely necessary. The primary emphasis needs to be that of keeping any estimated error less than the user's request. If for some reason the estimated error becomes too large , over- head will have to be expended to rectify the problem. If the error becomes much smaller than the user's request, it might be desirable to modify the method — provided the overhead is less than the expected savings in time . There are three ways to satisfy these goals: 1. modify d (x. ,t) by changing the error request given to the subprogram solving the system of equations (DIFSUB), 2. modify d (x. ,t) by changing the spacing in the partition, 3. modify d (x. ,t) by changing the order of the spatial analog. .A. J- 31 When "beginning the integration, it would be desirable to determine first (if possible) those items involving the greatest amount of overhead. After the integration has proceeded a number of steps, the item involving the least amount of overhead or the item least apt to adversely affect the solution should be changed first. Care must also be taken that changes do not occur too frequently since modifications will, in general, have a tendency to introduce additional errors. Of the three items above, the first has the most desirable properties. The change of this parameter could cause the vectors in DIFSUB to be multiplied by some factor depending on the change required in step size. On the other hand, the last two involve a greater amount of overhead and also involve interpolation that could introduce larger errors than those involved with changing the time parameter. Consider the problem of obtaining an initial partition or of changing to a new partition. The problem is to determine the order of the method and the points to include in the partition so that the integration may be accomplished most efficiently and yet have the desired accuracy. If it is assumed that given a partition, the majority of work is involved with the function evaluation or calculation of derivatives , and in the solution of a system of equations required by DIFSUB, it would be possible to develop a sort of work estimate based on the number of multiplications and divisions present. In the case of the derivative calculation, the two things to consider are the calculation of the derivative and the calculation of its error estimate. In the first case, we are using the Lagrangian, Hermite, and trimite three-point formulas, so that the coefficients are known as socn as 32 the step size is given. In the case of the error estimate, the coefficients for its interpolating polynomial are calculated and then stored for a given partition. When this estimate is. calculated, there -will he four inter- polating points for each derivative (-with the exception of the second derivative using the Lagrangian method which requires five points). The following table lists the number of operations used in calculating the derivatives for each method. The derivatives requiring calculation are the first and second for the Lagrangian, second and third for Hermite, and third and fourth for the trimite. Multiplications/Divisions In Derivative Multiplications/Divisions In Error Estimate Lagrangian Hermite trimite lower higher h 5 8 8 12 12 lower higher 1 2 k 5 7 8 Table 2 . Number of Operations in Equation Evaluation To aid in determining which method would be most efficient, the user is to indicate the derivatives required for each of the methods used. In solving the system of equations, the procedure is divided into two steps : 1. decomposing the matrix into upper and lower triangular matrices, 2. using the decomposed matrix to solve the system of equations. In each of these steps , the number of nonzero multiplicative operations required for each equation is proportional to one-half the number of of f -diagonals . In the case of decomposing the matrix, this proportionality 33 constant is two. In the case of the solving step, the constant is three. The number of off -diagonals is two, five, and eight for the Lagrangian, Hermite, and trimite, respectively. The average number of operations per equation for each method with each step is given in Table 3. Decomposition Solving Lagrangian 2 o Hermite 5 trimite 8 12 Table 3. Number of Operations in Decomposing and Solving The problem in trying to bring the three sets of numbers together is that the relative frequency of the number of function ( Nf ) , decomposition (N d ), and solving (H q ) calls depends on the nature of the problem. The number of function evaluations can be given directly as a function of N d and N bv s J N = n *N + i. N f m d 2 s where r^ is the number of equations in the system. The reason for the factor of \ is that for each call to obtain the function evaluation, the solving routine is called for use by the system of differential equations and for use by the system of asymptotic error equations. The first term arises from the fact that DIFSUB calculates a Jacobian by a numerical differencing method that calls for the function evaluation for each equation in the system. The ratio (R) of N d to N g is variable and depends on the solution to the system of equations. In systems where the Jacobian varies slowly, the ratio ranged from about .033 to .0^5. Since those equations were very well behaved and since an increased ratio indicates an increased 3^ amount of work, a value of .05 will he used initially. The ratio will then he modified as the integration progresses . At the start of the integration, the value N, could he written as N, = .05 N d s An estimate of the work required for each method can now he obtained as a function of the number of calls to the solving routine. Such an estimate is given by W = n {^7-W.. + R(W •-— + W.) + W }«N m m 2M fl f2 M d s s where n is the number of equations, W, and W are the number of operations m as. in the decomposition and solving subroutines. M is the number of values stoi at each point, and W , W are the operations required in a function call fc a step in the solution, and a function call in determining the Jacohian. The latter distinction is brought about because the process of calculating the Jacobian does not require any error estimate calculations. In order to decide which method should he used, the values of W /N for each of the methods are compared. Although one is interested in m s using the method with the smallest value of W /N , there is an additional m s complication in trying to decide which method to use. This involves the fact that if an attempt were made to change the partition during the integration to a higher ordered method, the interpolation polynomials would be such that error estimates would be too large to work with. As a result, it might be appropriate to bias the selection slightly in favor of the highei ordered methods during the initial partitions selection. After integration has begun, if a change is to be made in the partition, the simpler lower ordered methods would be preferred. 35 Once the partition has been selected, it would then be possible to adjust d t to some extent so as to keep the error close to the user's request. If DIFSUB is asked to produce a large d t then it will increase the step size or lower the order of integration. If asked to produce a smaller error, it will decrease the step size. The modification of d t and d x should not occur too frequently since the introduced errors could become large enough that a satisfactory solution is impossible. 2.k Initial Partiti on In order to start the integration, it is necessary to have an initial partition. Since any modification in the partition will introduce additional errors, it is important to find a partition that will be the basis of an efficiently obtained solution that will satisfy the user's request. One way to develop such a partition would be to control the amount of error present in the derivative. This would lead to trying to obtain a spatial error parameter that would be used to calculate this initial partition. There are three basic approaches that could be followed. One, the user could be asked to provide this parameter. Two, a set function of the user's overall error request could be used to develop such a parameter. Three, an arbitrary set of partitions using equal spacings could be set up and then tested to see what error results in the asymptotic error formula The first method of the three has the least desirable properties from the standpoint of the user. It requires that he make some sort of 36 judgment about the equation. If his judgment gives poor results, then the parameter would have to be modified. Since the program would eventually have to modify this parameter, it would be appropriate to generate the parameter initially. The second method has. the advantage that the function could be an easily calculated value. For example, one could have this parameter be one-tenth the value of the user's request. If this generates a partition that creates, valid errors, then no unnecessary work has been done . If the errors generated are too large , then the partition would have to be recalculated based on some formula of parameter modification. The third method would start with a partition using a step size a— b of , where a and b are the end points. Using this partition, the 2 m integration would be performed over one step. The error in the solution would then be compared with the errors in the derivatives. The process would be continued until the errors in the one step of integration correspond with the user's request. The maximum error estimate in the derivatives would then be used to generate the working partition. The disadvantage with this approach is that overhead will always exist. If the first partition yields a value for the spatial error parameter, the working partition would still need to be generated. The exception to this would be the case where the working partition contains equal spacings. Since the second method is simpler to work with and since it can, if required, perform somewhat like the third method, it was decided that the second approach would be used. 37 Once an error parameter has been determined, the working partition needs to be generated. Two approaches were considered. One, it would be possible to start at one of the end points and then choose each successive point so that the error condition is met, Two, it would be possible to start with three intervals and then subdivide each interval in which the error condition is not met. In the first method, one would start by adding two points to the partition, a center point and a second end point of the first subinterval. Thes.e two points, would be chosen so that the error condition would be met. The next step would be to use the second point above as the center of a new subinterval, and if necessary, add points on either side of this point so that the symmetric three-point formula would again satisfy the error condition. The process is continued until no new subintervals need to be added to completely cover the entire interval [a,b]. Some of the disadvantages of this procedure are 1. each time a new center point is added, it may be necessary to backtrack in the interval to handle newly added points, 2. even if all step sizes used in the subinterval are powers of two, there may develop problems terminating at the second endpoint , 3. the possible backtracking could create programming complications In the second method, start with the three points P = a + ^^ P = n + ^ZiL 1 2 ' 2 a k ' b-a p 3 = a + 3' -jp. (see Figure l) a P I I 2 P 5 ) 1 F b 3 Figure 1 38 The subinterval [a,P ] is next checked to see if the error is reasonable. b-a If not, points are added on either side of P 2 with step size -g— . If the error was acceptable, then this fact is recorded so that on future passes through the interval, this subinterval will be bypassed. The next step would be to check the subinterval [P ,P ] to see if the analog for P 1 satisfies the error condition. At the same time, the newly added points P. and P q are checked to make sure that they also satisfy the error condition. This process is repeated until all subintervals satisfy the error condition . In order to find out if a given derivative satisfies the error requirements, it is necessary to generate some sort of error estimate. An attempt is made to duplicate the type of error estimate that would be made in the actual integration. In an attempt to guard against unusual initial conditions, the b— a ., . . , testing at each point is done twice — once with an h = — — and once with 2 h = a— —-where .9 < a < 1. 2 n 2.5. Changes in Integration Once the integration parameters have been set, nothing needs to be changed until one of two conditions occur. If the estimated error ever becom larger than the desired value, then some sort of corrective action will be required immediately. If the estimated error becomes much smaller than the requested value, it may indicate that a change is desirable. Consider the case where the estimated error is larger than desired. There are two possible actions — either the error parameter to DIFSUB or partition could be changed. 39 The first action to be taken would "be to attempt to change d . This would be done provided the difference between the estimated error and the error request is smaller than the error reduction possible from DIFSUB by changing d^. If this condition cannot be met, or if several changes in this parameter have already occurred, then it will be necessary to change the partition. Although it would be possible to add selected points, it was felt that since the changing of the partition is being done as a last resort, the desirable action would be to change the entire partition to one with smaller step sizes. There are two other advantages in this course of action. First, the program used to change the partition could be basically the same as that needed to modify it. Second, since the entire partition now is satisfying a smaller value for d , there is the possibility that the partition would not need to be changed as frequently. If the conditions above occur during the first several steps, then it is assumed that a finer partition should have been used initially. Rather than create a new partition, it was felt that a better solution would be obtained if the process were restarted with a smaller value in the error parameter. If this process of restarting must be repeated a few times, then it is assumed that the problem cannot be solved with this method. If the estimated error is much smaller than that required, then a similar sequence of steps would occur as those just discussed. The first step would be to see if a change in the time error parameter is appropriate. The maximum size allowed for the time error parameter would be approximately nine-tenths of the user's request. If the error estimate continues to remain small, then consideration will be given to changing the partition. ko Since going to new partitions can introduce errors, there is an attempt to delay any changes to coarser partitions unless there is ample indication that such a partition will not have to be changed back. Just as in going to finer partitions, it would be possible to set up the program to delete selected points where the error estimates indicate it could be done. However, since precautions are taken in not changing too soon, it is felt that when a change is made to a coarser partition, the entire partition should be changed. Also, it would be appropriate to see if a lower ordered method could be used. The ideal situation in trying to go to a lower ordered method would be to be able to keep the same partition but just reduce the number of items saved at each point. Normally, when changing to a lower ordered method, one would expect to have to interpolate some in order to make up for the loss of values at the partition points. 2.6 The User's Equations The communication between the user and this automatic system of subprograms used to solve the user's equation (called the system) should be kept as simple as possible. Besides a program to drive the system, the user would be expected to provide subprograms that give the differential equation, initial values, and boundary values. In the case of the differential equation, the user would be expected to provide an equation of the form 2 3u „/ 3u 8 u . v W = f(u > 3? 72' X ' t} ' 9x The system would provide the values for the arguments of the function f , and his program would return the time derivative. Ul If the higher ordered methods are also used, then two other equations would also be included. These equations would have the forms 9v _ , 9v 9 2 v ^ 9x and 3w _ , , 9w 9 2 w , 9t~~ n(u » V > W > ii"» — * x > t ) 9x where and 9x - iJi 9 2 u V " 9x' V '" 2 9x Again, the system would provide the values for the variables of each function, and the user would return the time derivatives. A similar approach is used for the initial values. Considering that all ordered methods might be used, initial values would be required for u, v, and w. If u(x,0) ■ r(x), then v(x,0) = s(x) = dr ^ x) and dx 2 w(x,0) = t(x) = — ^ — = —^ — . If the derivatives of one of these do not exist, then either an approximation with continuous derivatives will need to be developed or the method will have to be restricted to those lower ordered methods for which the initial values exist. The problem of the boundary conditions is more difficult to handle The general form of the boundary condition is u(0,t) + u (0,t) - C(t) =0 at x = and u(l,t) + u x (l,t) - d(t) =0 at x = 1 k2 The user supplies these in terms of algebraic equations of the form r = u. + A(u. , x. , t) - d(t) 1 1 i where A(u. x. , t) is some representation of u x and r will he the error in satisfying the boundary condition. If the Lagrangian method is used, the function A will have to be determined by some form of interpolation. In the case of the two higher ordered methods, u = v. so that A(u ± , x ± , t) = v ± could be used. Since the higher ordered methods require additional values at the boundaries, it will be necessary for the user to provide additional equations involving these variables. It would be ideal if these boundary value equations could be obtained from the boundary condition of the original equation. The simple example of u = u , u(x.O) = sin(nx), u(0,t) = u(l,t) = . t xx shows that this will not work. The solution of this equation is given by 2 u(x,t) = e _7r sin( x) . From the definition of v(x,t), the solution is found to be v(x,t) = TTe" cos(ttx) Thus, -ir 2 t v(0,t) = -v(l,t) = Tre which cannot be obtained from u(0,t) = u(l,t) = One possibility for handling the problem would be to use an equation of the form r = v. - B(u. , v. , x. , t) l i' l l where B is an approximation of v. which is developed from interior points. 1*3 3. PROGRAMS The purpose of this section is to provide a description of the programs used in developing the automatic solution. There are eleven such programs. Eight of these programs (PARDIF, XSTP, DIFFUN, DERIV, DIVERR, DEC0MP, S0LVE, DIFSUB) are set programs, that make up the method. The remaining three (BDRY, INIT, UFUN) are user-supplied and provide information about the user's equation, initial values, and boundary conditions . The user also needs a driving program which would handle the printing. In order to use the system, the driving program calls PARDIF for each step in the time direction. As a result, the user can determine when and what information he desires to print out, and when to stop the integration. These programs are in Appendix A. 3.1 PARDIF PARDIF is the main program of the method. It has the basic purpose of integrating one step for each call made to it. In performing this step, it must cause an initial partition to be generated, create functional values for this partition, integrate the step, decide whether the step was acceptable, and decide whether modifications should be made. The creation of a partition is done in two steps. First, it calls XSTP which provides as output a pointer to indicate the proper order of points, and a vector which contains partition values as well as an indication of the other points which are used in analog. The second step consists of placing values in the proper order in the array x and placing two values in a second array (INT). The one value indicates the relative kk position of one of the other points in the analog. This will provide sufficient information since the partition is generated in such a way that two of the three points of the analog are adjacent. The pointer, therefore, is used to indicate the third point. The other value in INT is the starting position for the interpolating polynomial of the error estimate. Once the partition is well-defined, the value of the function at each of these points needs to be found. If this is the initial partition, then the values are obtained from an initial value subprogram (INIT). If the partition is a modification of an existing one, the program tries to find the values in the old partition. If this attempt is unsuccessful, then the values are found by interpolation by calling a subprogram (DERIV) . . Since there are several different vector components that would use the same interpolation, DERIV is asked to return the coefficients of this polynomial. The interpolation of all components is then accomplished in PARDIF. Every time a new partition is generated, the program calculates — by calling DERIV — the coefficients used to calculate the error estimates. These coefficients are then stored in an array (DIVER). If all the values associated with the partition have been defined, or if this step does not modify any of the parameters, the program is ready to integrate and check the results. The integration is done by calling DIFSUB. On returning from DIFSUB, the program first checks the return code to find out if the step was completed satisfactorily. If not, the program first tries to decrease d by decreasing the time based error parameter (TPS). If this procedure is unsuccessful, an attempt is made to change the partition. ^5 Whenever the return code from DIFSUB indicates that the integration was successful, PARDIF determines whether any additional action is to he taken. To accomplish this, the program determines the maximum error (RMAX) introduced in the step by looking at the results from the asymptotic error equation. First, consider the case where RMAX is larger than the user's request (EPS). First of all, the number of steps is checked to see if it is within ten steps of t = 0. If so, the time based error parameter (TPS), and the space based error parameter (SPS) are reduced by a factor of ten, and the integration is restarted from t = 0. This process is repeated three times before taking a terminal exit. When RMAX is larger than EPS, an attempt is made to decrease TPS until it has a value less than 10~ 5 . Once TPS has reached this value, the space parameter (SPS) is decreased by a factor of . 5 in an attempt to change the partition. The desire is to add enough points to cause the error estimate to be less than EPS. In the event that RMAX is less than EPS, consideration can be given to modifying parameters so as to decrease the work. The program DIFSUB will automatically change the order of the method it uses and the step size associated with that method to keep the error close to parameter TPS. No changes will be attempted in TPS or SPS unless two conditions are satisfied. One, the last modification has not occurred in the previous ten steps. Two, the error introduced in the step must be less than some ot times EPS. If the two conditions have been met, TPS is increased by the minimum of U6 logTPS + lofiEPS gmd Tps + i. (EPS.RMAX) . 2 2 Although many other methods of changing are possible, this choice appears to change TPS rapidly without causing the errors to become too large. After three such changes are made in TPS, consideration is given to changing the partition. If SPS is changed (and, therefore, the partition), an additional restriction is imposed. This restriction is that RMAX be less than a 2 -EPS. If this condition is met, then SPS is given the value of the minimum of 2-SPS, SPS + ,5(EPS-RMAX) and .U-EPS. If a change does take place in SPS, then TPS will be given a new value by TPS = . 5 • ( EPS-SPS ) . Again, the changes represent one choice of many possibilities. This time it is desirable to increase the parameter rather slowly so that the partition does not change too much. Consideration will also be given to changing the partition without changing SPS provided the following set of conditions is satisfied, 1. at least 30 steps have occurred since the previous modification, 2. RMAX has remained less than EPS, 3. each component of the spatial truncation error is less than a 2 -EPS. 3.2 XSTP XSTP has the function of generating the positional values of the partition under two different situations . The first occurs when the initial partition is generated and the second whenever an existing partition is modif. hi In both cases, the program begins with partition points at .5, .25, .75, and a step size h = £. On each succeeding pass through the interval, h is divided by 2. Until h <_ HMAX (the maximum step size the user desires), the program places points at (2i-l)«h, 1 < i < — . After — — 2h" h > HMAX, the program places points only when the error estimate in the analog does not satisfy the error condition. These additional points are placed on either side of the test point. Because of the way the partition is generated, each point that is added vd.ll be halfway between two existing points . If the analog at a particular point satisfies the error request, this fact is recorded in the vector NCK. Before trying to add points, the program first checks NCK to see if the testing can be ignored. When all values of NCK are set to +1, then the procedure is terminated. The differences between the two uses of XSTP take place in obtaining the functional values and the error estimates. In obtaining the functional values for an initial partition, the program calls the user's initial value program (INIT). The case of modifying an existing partition is accomplished by the interpolation of values of the previous partition. For estimating the error of an initial partition, the program adds two points whose positions relative to the points of the analog, are fixed. Tliis allow-, the coefficients to be calculated for all cases. When the partition is changed, the functional values are not as easily determined, so the error estimate is obtained by using values determined for the partition and calling the subprogram DERIV to interpolate. Once the partition for a particular order has been developed, the work estimate is obtained using the ideas of section 2.3. The entire process is then repeated for the next ordered method. 1+8 When the partitions for all possible .orders have been developed, the work estimates are compared to determine which method should be used. This is done by asking that between two adjacent ordered methods, the smaller of aw w (where m is the lower order) indicates the order m ]. m 2 1 preferred. The a is used to set the biases desired in the selection of a method. Once the method has been chosen, the program then checks each analog to make sure that the smallest available step size is being used. It also prepares output vectors indicating the natural order of points and the points that will be used in each analog. 3.3 DERIV DERIV is a double purpose subprogram involved with interpolation. The one task is to calculate the coefficients for the error estimation polynomial of a given partition. The other task is to provide interpolation polynomials for calculating new values of a partition. The error estimates needed are given by (27), (28), and (29). Each of these involves the factor P'(x), which can be determined as a function of the spacing h. Since n = 3, P (x) = (x-a^ (x-a 2 ) (x-a 3 ) , and x = & 2 , the value of P' (x) is given by P'(x) = (x-a )(x-a g ) + (x-a 1 )(x-a 3 ) + (x-a 2 )(x-a 3 ) = (x-a )( x-a _) = -h 2 U9 Thus, the error terms become -T- f (n) and -rr f (n.) for the Lagrangian |gg- f^(n) and 1^- r T) (n) for the Hermite and ZiMo f(9)(ri) md 151200 f(l ° } ^) for the trimite The problem remaining in determining the estimate is the evaluation of f (n). In "the case of the Lagrangian polynomials f (j) (x) = E *f j) (x) f(x.) (31) i=l X X In order to obtain an estimate of f (x) , and f (x), the &.(x)'s were chosen so that the proper derivatives existed. In the case of the third and fourth derivatives, k and 5 points were required, respectively. Using these polynomials gives f (3) (x) ~ E -t-^ f(a.) i=1 T (x - a j } J=l and f U) (x) " I — 2i f(a.) 1-1 — (x- &J ) The program then generates the coefficients of f(a.) for use by other programs In calls from XSTP, the actual value is desired, so the coefficients are multiplied by f (a. ) and the value returned. 50 In the case of the Hermite polynomials, the problem is to obtain estimates of f (x) and f (x) when given that n f(x) = 2 h„(x) f(a, ) + h.(x) f (a.) (32) i=1 i l i i For this the derivatives of h. (x) and h.(x) are desired. The results are given by (assuming four points) h[ 6) (x) = 20(l-a!U.)(x-a i ))(il" , (x)) 2 + 120(-2£.!(a i ))U"(x) V"U)) hf T) (x) = lU0(-2£!(a.))(£"'(x)) 2 i 'ii and hj 6) (x) = 20(x-a.)(.r(x)) 2 + 120^(x) ^'.'(x) hj T) (x) = 1U0(J?'.'(x)) 2 In the case of the trimite polynomials, the desired estimates are f (x) and f (x). The polynomial f(x) is given by n f(x) = l (t.(x) f(a.) + t.(x) f'(a.) + t.(x) f"(a.) (33) i=l so that the ninth and tenth derivatives are desired for t. , t. , and t. . Each ii l of these has the form V = A(£(x)) as developed in section 2.2 where V is one of t. , t. , or t.. V^ 9j (x) and V (10 '(x) are then given by V (9) (x) = (l68o(£*"(x) 3 )-A + 15120•Ji"(x)(r"(x)) 2 •A , + 36(630U"(x)) 2 V"(x) + U20£'(x)(£"»(x)) 2 )-A n V (l0) (x) =16800(2, (x)) 3 -A» + T560QH"(x)U (x)) 2 -A" It will be noted that in all of the cases above the estimates turned out to be functions of the Lagrangian polynomials. Therefore, the program 51 generates ^(x). £ «(x), £'(a.), £»(x), £»(a.), *"( x ) as needed ^ then combines them to form the desired estimate. The second task of this subprogram is to interpolate to find values at new partition points. Since these interpolating polynomials are also functions of the Lagrangian polynomials, much of the code used above can also be used here . Since new partitions will only involve the same or lower ordered methods., the Lagrangian method only requires (31) with j = 0. For the Hermite polynomials, the first derivative and the functional values are needed. The functional values are obtained from (32). Since h. and h can be i i written in the form V = A£ 2 (x) where A is [1-2SL ! (a. ) (x-a. ) ] and (x-a. ), respectively, the first derivatives can be written using V'(x) = A'£ 2 ( x ) + 2M(x) £»(x) Finally, the trimite method will require the first and second derivatives as well as the functional values. The functional values are obtained using (33). Using the form V = A 3 (x) for the fs (see page 15 for values), the first and second derivatives can be obtained by V'(x) = A'£ 3 (x) + 3M 2 (x) I »( x ) and V"(x) - A"£ 3 (x) + 6A'£ 2 (x) IHx) + 3A[2£(x)U'(x)) 2 + * 2 ( x ) *"( X )] Since calls to DERIV for purposes of interpolating new partition values may require the coefficients in some cases and values in others, both of these options are made available. 52 3 . h DIFSUB This program is a modification of a program by Gear [3]. The modifications basically deal with four areas: 1. the presence of the asymptotic error equation, 2. the use of DEC0MP and S0LVE, 3. the use of algebraic equations for the boundary conditions, h . the presence of externally changed parameters . Probably the largest number of modifications deals with the area of the asymptotic error equations. The main item here is that every time something is done with the solution vector (Y), the same thing is also done with the asymptotic error system. Since the solution to the asymptotic error system (HMY) does not need to be as accurate as Y, coefficients were developed by Gear [k] so that the first two coefficients would be the same in all methods, and thus, allow the asymptotic error system to use a lower ordered method. This can be accomplished since the solution involves the Newton iteration z 1 j.->\ = z / \ + & [-TT— £] -1 F(z , J (3*0 n,(m+l) n,(m) — 3a. — n,(m)' VJ ' where F( z / \) = hf(u) - hu' in this formula, n,(.mj — W L 3a - J lZ l M 03u J where I and SL are the first two coefficients of the method. Since W is the same for all ordered methods , the two systems can be solved using different orders. 53 Another slight modification required hy the asymptotic error system is that in order to start integrating, initial time derivatives are required. This means that the evaluation of the Jacobian must be present earlier in the program. The need for the second modification is discussed in Gear [3]. He noted that the break point for using DEC0MP and S0LVE occurred when about five or more equations are in the system. Since the number of equations solved here is in general eight or more, the need is apparent. As an added reason, the presence of the asymptotic error equations doubles the number of calls to S0LVE whereas the number of calls to DEC0MP will remain about the same. Two changes are required, by the presence of algebraic equations for the boundary conditions. First, in equation (31) the function F is replaced by G( Uj .) which is the residual of the algebraic equation. Second, the calculation of those rows of W involved with algebraic equations will need to be modified. In this case, the component is given by da — 3a where g is the algebraic equation. The absence of I is due to the fact that derivatives of the variables are not present in G (see Gear [5]). In the fourth area, there are two conditions that take place: 1. change in the error parameter 2. change in the partition. In the first of these, the only requirement is that various error condition values be recalculated. As a result, a jump is made to that portion of the program. 5h The second is more involved. If the partition has been changed, then it will "be necessary to change those variables dealing with the number of equations in the system, and to re-evaluate the Jacobian. 3.5 DEC0MP and S0LVE The purpose of these programs is to solve Y = Ax for x where A is a banded matrix. The procedure is carried out in two steps. (The programs are modification of ones given by Forsythe and Moler [6].) The first step is to decompose the matrix in such a way that the lower triangular portion will represent the operations done to the vector Y in the first half of a Gaussian elimination procedure, and the upper triangular portion will represent the coefficients of the variables used in the back substitution process. This step is done by DEC0MP and is called each time the matrix A is recalculated. The second step consists of two parts. The first takes the lower triangular portion and carries out the operations of the first half of the Gaussian elimination on the vector Y. The second part performs the back substitution. In these subprograms, the matrix A is stored in a double precision array (UL) diagonal by diagonal, with the main diagonal stored in UL( K ^ +1 , I) where KT is the number of diagonals and is assumed to be odd. 3.6 DIFFUN DIFFUN calculates the time derivatives of the system based on the differential equation. It does this by calculating the various space 55 derivatives required and then supplies these values to the differential equation subprogram (UFUN). At the same time, it also calculates the error estimate of the space derivatives by using the stored coefficients of the interpolating polynomials for such estimates. These results are then added to the derivatives and given again to UFUN. The difference between the results of the two calls to UFUN then provides an estimate of d x (x jL ,t). A second use of DIFFUN is to help determine the Jacobian used in DIFSUB in the predictor-corrector method and in the differential equation for the asymptotic error. In this environment, the error estimates are not used so a variable is set which causes this calculation to be bypassed, 3.7 DIVERR The asymptotic error equation was given by (section 2.2) e»(t) = f (y(t))-e(t) -d-d y t x This program calculates the time derivative e'(t) for each variable using f y (y(t)) which is stored in the array PEQ, d which is given by ERR0R(IA ) PERTST ' and d x which is st °red in ERR(IA). 56 k . RESULTS The system was used successfully to examine three separate examples. In each, case, the problem was solved using the system and an alternative numerical technique. The three problems were: 1. u = u , u(0,t) = u(l,t) = 0, t XX u(x,0) = sin(ut) (35) 2. u x = u + u + u, u(0,t) = u(l,t) = 0, t XX X u(x,0) = e" X sin(.TTt) (36) 3. u = u , u(0,t) - u (0,t) = t XX X u(l,t) + u (l,t) = 0, u(x,0) = sin (ttx) (37) Example 1 In this example, the solution is quite simple to calculate and is given by -TT 2 t u(x,t) = e~ sin(iTx) (38) The independent solution uses a method developed by Crank and Nicholson [Tl The programming of the method is based on information in Richtmyer and Morton [8]. The program is listed in Appendix B. In solving the problem with the latter method, an answer was obtained using a step size of Ax = g and At = Ax. The process was then repeated with step sizes of -7-, r^", ZJ7»-" until the error was less than .001. Results are in Table h. 57 00 VO CM -=*• CO VO 00 VO CJ LTN H CM OJ hIoj ii X o •H -P H O en oo l -3- VO On OJ 00 I w oo CO O -3" I w on i OJ vo I w o ON I o H I ON -3- I t— LT\ CO I w OJ o OJ oo CO 0- o 1 H Ol w 1 1XN H W CM t— on O -p C o co H O o ooloo 5h o u on i On ltn oo i on r -3- i CO LTN. CO o H I w CO CO H I w oo On 3 03 O CO o •H -P H O CO CO OJ CO o OJ I H vo VO £ CO I w OJ t— On H I w H H O O O h|j- 00 I w VO OJ I I w f VO I w 00 vo I o H I W VO CO H I pa O OJ I OJ hIco II § •H -P H O to G O ■H -P H O to LPv O OJ 00 I w ITS OJ H H H O VO O 00 I W j- H OJ 00 o OJ I w o LTN I w -3- OJ OJ OJ I w VO 00 VO I w -4- oo ON H CO w OJ I w 0O ON I H I w CO H I w oo t— I OJ o OJ I 0J o OJ I w OJ o OJ w OJ H •H CO Ov H H w I 1 1 -p ON O -3- _=J- OJ w w w VI n CVJ H vo ir\ CM vo H LTV H H o 1 1 1 1 CO II CVI CM CM CO CO J- CO U O I 1 1 I 1 1 1 ^ !h ON vo CO OJ o ir\ H w Sh CVJ CO CM LTN CM LTN H -d- n CD 1 1 1 1 1 1 1 H CU H a CM CO CO -3" P< o 1 1 1 1 ■J* vo •H vo H H w w H I 1 CO P> c— O ON o O LTN 3 w X 3 OJ H l/A LTN CM LTV -3- W H o 1 1 1 1 ^ CO O

w W H LTN 3 CVJ o -4- CO H -3- CO CM H CD o l 1 1 1 H CO •a CVI CM CO CO -rt EH fi 1 | 1 1 1 J- WN °? O w w W w W 1 1 u CVJ ir\ vo H _=r w H w u H H ON CM oO CM _=r -3- CD 1 1 1 1 1 i G CM CO -3" O 1 1 1 -4- L/A CO ■H -=r CM m w H 1 | 1 P H -^ -x H _3- w w w o H O CM CM CO CM _=f -3- | 1 1 CO H CVI OJ CO ON VO LTN ir\ ITN H OJ J- O Lf> -3" ^t -=f CO vo H 6o -p

.9 +° CO CM CQ CM " J .TV ON M CM oo H 53 «H O • -O CO cu O* P cu cd O ro CM H CM CM O ro H VO o -=r VO OO LTN CO UN CO LTN vo un H vo CO OO o 6 cu > •p =tt cu en CM CM l oo -=t- VO t- o H Ov CD cu -P •H hIcm o u u 4) H pq H H r oo vo J- LTV * * 1 1 OO LTN f CM r pq ON • pq ii X a CM .-=*• | | VO 1 t— o H ON CQ o •H -P 3 O oo ON o 0O H \o oo oo LTN pq CM 1 pq ON pq H H O H • • * | 1 • O CO II 00 CM oo -=f 1 1 vo I t— Ov o H CQ o U 1 pq CO ON 1 o H pq W O CM _=f LTN m ON pq CM • 1 9 pq co ft oolco CU | r r i" i r l H H CU X s CM ^ VO tr ^ o H H & o ■ri P 2 H CO CM o ON Pel W CO ON LT\ CM pq ON 1 pq CM i 3 pq CO pq H • | 1 f O U CQ o CO oo OO ^ vo I t— o o H CH 8 •H P U O U pq VO t— H CO oo -=r pq CO 00 1 3 pq 1 pq LTN H|-=f CU i I i 1 1 + 2 H || o X a CM ■=* VO 1 t— o H O H CO o •H P 2 t— H CM ON vo o H pq LTN CM J- CM CO oo 1 1 pq J- pq LTN VO .H • • * * 1 i CU o H en ■s oo 00 1 OO -3" VO 1 t- t— t- o u u 1 pq CM pq -=t- pq H vo H H CM pq H pq CM pq oo pq oo hIco cu • • 1 1 1 + II X a CM -4 I 1 VO t- Cr ty o •H vo H H 00 o J- CM CM H t— H pq CM • I pq oo pq 00 H • • • • 1 O en O CM H LTN 00 CM LTN O H O oo o 00 CM t— CM CO UN H 6l CM en C\J OJ CM CM CO U ft CD CD g CO a <* o CM NO CO CO o H VO CM H ON CM CO CD o< -p =»fc CD hIcm o •H p H o CO LTN -J" CM CM CO CM On On ltn £ CO I ON CO o o CO on on CM i W CO CM a CO CO I on CO I CO CM -3- on 0\ CO I NO CO I w -3- o On H O CM H H I CM O NO o H CM H W CO i CO CO CM ON o CM H I CM H on a cu CO CD p bO C •H co p o O onlco I CO CM I LTN i w t— NO I* LTN I w NO CO I w H I w CM I CM NO CO Oh m o •H -p H O CO H NO CM CM CO O CM I w on t— w NO CM CO I on l CM I H H CM H I W NO o u Ld- CD )|cO (3 O •H P H O CO NO CM O 3 -J- I w CM CM ON LTN on NO o -3- I w LTN m I -3- I w CM LTN CM I w NO LTN CO CM LTN I CO I o CM LTN I w on co I 9 On W ON On on NO CM H w NO CM NO H I w NO CM H I w NO CM H I W NO H H LTN CM H I W CO CM CO H H W -=J- -=t LTN t— o 1 I w 1 H S W 1 o u OO o 00 VO ON LPv VO •H u CD LTN ON LTN oo I CM 1 H CM co 1 H o o •H CM I i" -3" 1 c— CO -p H H W w m W m w • ^ ON LTN CM 00 ON 00 LTN VO ON H O w H O -=t- -3- CM H CM 1 CO cu H ■s _=j- _=}- l^ -3- -=f ITS t- EH o 1 -=f W W 1 a W 1 *H o W VO CM LTN H -=f M LTN t- _=r 00 CM H CM t— CD 1 1 1 1 i O CM -4" _3" J- LTN i •H ■p LTN _=f- 1 1 w W W w 2 vo J- t— t— LTN Pi H -=f H H O oo 00 CM CM t— O en CM LTN CO 00 1 1 CM LTN o VO VO VO VO vo H CM L/N ON ON ON On On 00 t— H o u > 00 cu ,£5 P CU 1 1 1 (30 § •H -P 2 o c— OJ 1 -3- 1 i w LA LA 1 vo 1 •H CO o OJ OJ OJ LA VO O H A H O CO H o OJ OJ H 1 VO 1 OO H H O h O 1 LA 1 J- 1 LA 1 LA 1 VO II CO ?H LA vo H -4- t— OJ oo ao H CU LA ao LA 00 1 OJ H 1 vo + 1 ft OJ PS O •H OJ -d- "T LA LA | VO 1 CO H & -P 3 C\J OO On 3 3 OJ oo ao H O CO H O J- _=r OJ 1 H vo + 1 U M -f "f -T LA i -* j- LA o O m OO o £ VO ON OO OO o LA ON LA oo 1 OJ 1 H H + H CI O •H OJ | -J- -3- | J- -3" LA o CO « H ON H LA OJ $ I w ON 1 00 1 -3- CO H O co H O -3- -3- OJ 1 H H + H O H cu H -3- LA _=!- -3- -=!• LA o 1 1 w VO 1 w OJ LA 1 3 1 -3- 1 W LA h cu LA b- -4- 00 1 OJ H 1 OJ + OJ a o ■H OJ -3- -^- J- J* LA -P 3 LA VO -3- c— {a ■g i LA H o H o OO oo OJ H OJ OJ CO CVJ LA OO oo 1 1 + CVJ LA o vo VO VO vo vo H OJ LA ON ON ON ON On H LA H •p •H O CO a 1 CD oo on oo en ro CO po oo 66 to ooU li o =j O CO -=f c— CM -=J- VO J" VO ON VO _=T LfN vo o oo o H i w o o -=t- o -3- I CM o -=r vo 00 us I On t— ON ON H CO VO OO o H VO CO CM Po CM -J" I LfN CM ON VO CO o oo LfN o H VO I vo VD I VD H O O CO O H CO I H VD CO VD s cu -P CO !>> CO CD ,S -P bO fl •H CO CVJ g •H H O CO LA 00 CO VD co VD O H in LfN I w LfN 00 LTN I oo VD t- LfN H CM I LfN LfN I LfN 00 LfN I LfN H w w o H 00 H VD t— w H o H 00 H oo|co II X hUj- II C5 o •H H O CO o •H I o co LfN ON H VD VD H CM LfN O CO t— LfN o i 00 LfN oo -1 VD CM I W CO OO LfN t— OO if CM 00 -3" I w H -3- LfN I VD OO VD 00 LfN I 00 CM oo LfN I VD H LTN I VD H 00 H LfN I m 00 H VD OO VD I OO VD OJ "9 CM I CM H I CM I o I o H H CU H •3 EH H CM LfN CM VD O LfN LfN LfN ON VD H VD oo t— VD VD LfN H 67 o u CD CO -p C •H O CO CD 00 00 00 00 oo no oo oo to U ft CU CD >9 ^ oo ir\ VO oo H LfN 00 On S CO OJ CM CM OJ 00 00 J- LA C «h o s CD CO >> • T3 CQ CO CD o< -P oo H u\ oo LfN H oo f— CD +3 CD CCJ OJ -3- ON H ON OO vo OJ 2 t— VO VO vo 00 ■^f CO [•— O 03 !> =tfc CD vo t— c— oo o H H OJ w H H H H d •H CO ^> H|CM II oo|co hI-=i- h|co oo 1 VO II LfN OJ vo OO O 1 1 C~- CQ OO t— vo CO LTN ON oo W P-i -=|- CM LfN. o H CM ^t oo w -3" 00 OJ H O _=t- 00 oo 1 VO oo 1 it 00 CD H -* CO ft On CO 1 W t— W t~- CM VO -=f -=c oo H m 00 CM LfN O r-i CM LfN OJ o -3- OO CM H O -=1- H CM 1 cm CJ O •H O CO 00 1 VO -3- ON O VO CM 1 1 °? ON ON 00 OJ -4- ON LTN H , H O -=}■ O H O On CM 3 -3- OO CM H O -=t- OJ OJ CD H ,Q ct) EH 00 vo OJ CO -3" LfN 00 I W [*- O CO VO LfN oo oo ON W ON CO CM ON H CO O CO 00 CM CM O O oo t- OJ 00 H OO H 00 LfN LTN -3" 00 LfN ON CO VO O CM ON o o CO t— H vo H 68 CO p o -H o CD CO oo oo m oo CO CO u ft CD CD £> P -3- CO =*s » CO CD Xi P CO t— t— VO OO CO CO vo -d" _=r bD LfN ON o t- o fl H ?l -3- -3" LTN •H H H H H CO ►—i CM O o VO J- CM CO LA LA CM CM CO ON CM W J- t— o ON H J- O O VO LTN I w LT\ OO H colco II H CO LA co VO vo CO LA CM CM CM CO ON O H O CO W CM CO vo H LA H I I _=f 00 hI-3- H CO 00 VO LA CO VO CO H CM ON On CO O CO VO CO t— o VO LA I H I H CO H CD H ■3 EH l|0O H O VO CM CO o On On H CO CO OO o CM LA H O CO LA VO LA CO H O H I 3 LA CO H H O CM H CO LA CO o VO 00 o ON LA On t— H LA t— H LA H CM CO t— LA H 69 The independent solution was again obtained by using the Crank-Nicolson method. The program used to solve this system of equations was based on information from Ozisik [9]. The program listing is given in Appendix D. The results are given in Table lk . To obtain a solution using the CrankJJicholson method requires approximately six seconds. A similar solution requires about 3U seconds using the system. In examples 2 and 3 above, it was found that results were easier to obtain if the auxiliary algebraic equations at the boundaries were bypassed. This is done by setting up equations of the form \ =B(u., v., x., t.) r= v. -B(u., v., x., t.) where B represents., an interpolation polynomial. ' In order that the Jacobian will remain invertible , a 1 is inserted in the main diagonal of the equations bypassed. The dependence that interior points would have had on these boundary values, as present in the Jacobian, is now transferred to the interior points used to obtain that value. Although the error estimate used for the spatial truncation error appears to be accurate enough for the previous examples, the presence of a more accurate one may be necessary in more complicated problems. The system was used to try to integrate a problem where the initial value involved a peak in the center of the interval . Although the system was able to produce values of adequate precision, which would appear to indicate an analog of sufficient accuracy, the error estimates were large enough that the integration would have been rejected. (See Table 15) 70 oo VO CM hIoj I 1 H J- o CM _* _=f W H 1 ii 00 LT\ oo H CO LTN LTN W Ti Lf> VD oo o H O ir\ t- O X _=* OO CM H O 00 CM -P o ooloo 1 1 H W ON H o O CM w w 1 H ii VD O LTN O CO t- vo H O J- VD OO O H On -=r vo & X -=f CO CM H O ltn VD t— o 00 VD CM H |-=f 1 | H >H 00 t— VD H -=f W W 1 o l| t— LT\ ir\ V£> t— CO o H CM J- CM On H t— CM oo «■» Xi _=r oo CM O O ltn VO t- 00 R< rHlCO II ON OO ITS CM O CM CM OO LTN CM CM O CM U"\ CM ON VO CO H o o o CM 00 OO LTN O M3 I CO t— o CO CM H I t- vo O vo £ -3- H cu H ■§ EH 71 Hcvj ii X I |0O u o LA | LA | LA 1 LA LA LA w w 1 H 1 1 u -=1- OJ 00 OJ H ON H 1 CM 1 OJ 1 OJ 1 OJ 1 H 1 o CO 00 VO VO VO LA 3 •H vo LA CO oo 00 OJ rH -P LA VO t— -=f OJ On 0} 2 VO LA VO On oo t— > H LA -3" oo OJ OJ H O CO H H H H H H H a} •H -P ■H f-l CO t— 1 VO 1 VO 1 VO t— o W W w 1 1 H £ u 00 VO oo LA 00 LA •H u o\ CO H H H ON a> -a s H 00 -=»■ -3- J- -3- -3- O H H H H H H CD to H ! )L| t— 1 t— LA LA 1 -3- -J- w o H w w 1 1 W «H *-l On OJ t— VO t- t— o N -3- OJ H t— H OJ 0) • I 1 1 • 1 C O •H P 2 a o J- J- On t— H O o OJ t— 3 00 VO o CO •H on VO VO LA H P _=t- vo t— VO LA _d- 3 vo VO VO VO VO vo . H t- t— t— t— t— C— LA O o o o o O o H CU H ,a cd EH m i 00 1 00 1 oo OJ i OJ w w w w 1 On LA LA CO H oo rH -d- VO CO H H O o o o O o CM -=f VO CO O CM 72 A similar situation occurred with the following example: U t = 87 (l+u) fp u ^°» t) = u ^ 1 » t ) = ° u(x,0) = sin 3 (irx) x e[0,l] For use with the system, the equation was rewritten U t " ( 3^ + (1+U) 72 dX An independent solution was obtained by using a predictor- corrector modification of the Crank-Nicholson method that has been developed by Douglas and Jones [10], The solution was also obtained for a few values of t using the system with the error control portions bypassed. The results agreed with those obtained with the alternative method. However, the error estimates for the space analogs would have caused the integration to cease. (See Table l6) The results obtained from the first three examples appear to indicate the feasibility of using the method of lineswith the subprogram DIFSUB to produce an automatic program when the asymptotic error equations are integrated along with the original system. The two additional examples seem to provide the possibility of integrating other types of problems including those of a nonlinear nature. 73 S cd -p CO !>> CO o CD t— CO CO L/N 00 ,£] .1 — . o CO o ir\ J- H p HlOJ IT\ CD -* -=r t— J- t— -3- o t— -=r o (90 H o o H o o fi o a) p CO •H CO a CD -3" no t— !>> H a oolco OJ OJ oo co OJ ON oo CO 00 H O o no -3- o CD NO -3- o •H -P H o o P H o o H O CO S u t— a a •H -P O t— co d) OJ CD « Pi hUt -3- CO CO OJ o 3 O t— CO 00 On OJ i C\J 00 o Tj OJ 00 o CJ H o o O H o o O • A d • cd c o 0 C SAVE STORES Y PDFU60I C MSTRT, PDFUTOi C MAXAR LIMITS ON ANALOG TYPE PDF1180S C 1 - LAGRANGIAN PDF119Q) C 2 - HERMITE PDF1200? C 3 - TRIMITE PDF1210( C IDER DERIVATIVES USED PDF1220( C Y ANSWER VECTOR PDF123QI C X PARTITION PDF1240I C HMIN, PDF125Q* C HMAX LIMITS ON TIME STEP SIZE PDFl260t C YMAX USED FOR RELATIVE ERROR WHEN Y IS GREATER THAN 1 PDF1270C C ERROR TIME ERROR PDF1280( C KFLAG RETURN CODE PDF1290( C HMY ASYMPTOTIC ERROR VALUE PDF1300C C EROR SPACE ERROR PDFi310< C HSAV SAVES HMY PDF1320( C PEQ JACOBIAN PDF1330C C ERORA ESTIMATED ERROR INTRODUCED IN THIS STEP PDF1340( C PDF1350C IMPLICIT REAL * 8(A-H,P-Z) PDF1360C DIMENSION ERORAdl PDF137QC DIMENSION KL8(3I,RINT(2),NCK(1I ,D I VER (9, 1 ) , NPSAVE (1 I, PSAVE (1) ,SX< 3PDF1330r 1),NEQ(3»,IDER<4,3), INT( 2, 1 ) , X{ i I , Y( 5 , 1) ,YT(3I ,YTA ( 3 , 5 ) , XT ( 5 I , PDF139K 2ERROR(l>,EROR(l ) , YMAX { 1 ) , HMY ( 2, 1 > ,GAM( 3 ) , SAVE (7 , 1 I ,HSAV( 3 ,1 ) , PDF1400C 3PEQ(1 »,MAXM(3} PDF1410C COMMON /DVR/ DVY( 41 ,COEF( 12 , 3 » , ALPHAi 2, 4,3 » , DI V( 16 , 5 I PDF1420C COMMON /CNT/ ICSTPS ,NUMDE ,MUMSO PDF1430C DATA PVAL/.5/ PDF14'?0C DATA RVAL,SRVAL/. ID0,.01D0/ PDFl4iJi DATA KL8 /9,19,29/ PDF1460C PDF147GC AND CHECKS FOR CHANGE CONDI TIONSPDF148D0 POF149O0 PDF15000 PDF15100 PDF15200 PDF15300 c c c 2 THIS SECTION INITIALIZES COUNTERS IFUSTART .NE. 1 GO TO 10 MAXSV = MAXAR ISPSCT = TPS = .4 * EPS 77 SPS » .1 * EPS DO 15 IA*1,3 INP * 1 SEPS = EPS JSTART » KTPSCT = ICPART = KSPSCT = ITUPS = 1 10 IF(SEPS .EQ. EPS) GO TO 23 TPS = .4 * EPS INP=2 J=NEQ(MO) PDF15400 VEPS * DLOGCEPSI POF15500 Nl«7*MAXN*l POFISftOO PDF15800 15 MAXMCIA) * MAXN/IA -2 pnF „;; NUMA«MAXN PDF15900 P0F16000 PDF16100 PDF16200 POF1630O PDF16400 P0F16500 P0F16600 PDF16700 J^V X P0F16800 iC U PDF16900 NCNGE ' ° pnPWAon ITPSC - 0F17000 ICSTPS * NUMDE « NUMSO « T » 0.000 MAXER = MAXSV GO TO 200 PDF17100 PDF17200 PDF17300 P0F17400 PDF17500 PDF17600 20 IFCNCNGE .EQ. 0» GO TO 255 PDF17800 NCNGE " ° P0F17900 GO TO 12 PDF18000 PDF18100 ^' S = ? L ° G 1^ S> PDF18200 PDF18300 SPS=.1*EPS pnF1 " 12 IFilCSTPS.LT. 101 GO TO 5 PDF18500 200 CALL XSTPCRINT.NCK, DIVER, NPSAVE I Nl I , NPS AVE ( N2 ) , SPS, SMAX,SMIN,SX PDF18600 1»KFLAG,NEQ,MAXER,MAXM,IDER,MSTRT,PSAVE,INP,NUMA,M0LD,X,Y,M0) PDF18 700 IFCKFLAG .LT. C .AND. MAXSV .EQ. 1> RETJRN PDF18flOO IFCKFLAG .LT. 0) GO TO 299 PDF18900 IFCICK .EQ. 1 .AND, MO*NEQ(MOI . GT. NUMI GO TO 208 PDF19000 MAXER = M0 PDF19100 THIS SECTION PUTS THE NEW PARTITION INTO USABLE FORM AND PDF19300 OBTAINS THE FUNCTIONAL VALUES FOR THIS NEW PARTITION. PDF19V30 PDF19500 PDF19620 ,t k ,T.,» PDF19700 52-Soi? li7?I? PDF19800 JA-MO*( J+l»*l PDF19900 SAVE(1,JA I-RINTC2I JA=M0-1 IFCJA.LE .01 JA=3 PDF20000 PDF2C100 IFCJSTART.EQ.OI INT(2,l)«l PDF23?0ft DO 222 IA=1,J PDF20500 N3= M0*3*(IA-1)4N2-1 PDFPD600 IB=NPSAVECN3I PDF20600 id i»ri«vi:i(N3i PDF2070n INT(1,IA*U = DIVERCM0,IB) PDF2C800 78 222 C C c 250 208 300 IAA»l*MO*(IAI SAVEl 1, IAA) -DIVER I MOO, IB I N3»N3-MOOA IF( IA.EQ.l»NPSAVE(N3»»l N3»N30 IFUA.EQ.l)NPSAVElN3)»l IFl IA .EQ. 2) NPSAVEIN3I » 2 IFIIA .GT. 2 .AND. IA .LE. Jl NPSAVECN3) IF( IA.EQ. J. AND.MO.EQ.il NPS AVEl N3 )»J-2 IFUSTART.EQ.OJ INT (2, I A+l ) »NPSAVE(N3> CONTINUE IFUSTART.NE.OI NPS AVE! N30 l»NPSAVE(N3) NUMA=NEQ(MOJ IFl JSTART.GE.l) GO TO 300 NEW VALUES FOR INITIAL PARTITION. IQQQOO J2»JO DO 250 IA*1,J2 IAA1 « MO*( IA-1 I IAA = IAA1 +1 XI IA»=SAVEI1,IAA» CALL INITIYT,M0,XUAI> 00 250 IB=l,MO IAA » IAA1 ♦ IB HMY(1,IAA) » 0.0 Yll,IAA)=YT(IBI GO TO 258 MO = MOLD J * NUM NEOIMOI = J GO TO 258 310 10101 350 C C IA -I /MOLD THIS SECTION FIRST CHECKS TO SEE IF THE POINT PARTITION. IF IT WAS, THE VALUES ARE REUSEO. VALUES ARE OBTAINED BY INTERPOLATION. IW=JSTART+1 NUM=NUM/MOLO J2=J +2 DO 3 90 IB=1,J2 IA = 1 IAA=l+M0*(IB-ll XI=SAVEll,IAA» IFIXI.LE.X1 IA I I IA=IA+1 IFl IA.LE.NUM+2) KFLAG=-30 RETURN KFLAG=-41 RETURN IFl XI.EQ.XC IAII WAS IN THE EXISTING IF NOT, THE NEW GO TO 350 GO TO 310 GO TO 370 INTERPOLATION SECTION 79 KG»0 PDF26400 351 IF< IA.LT.NUM*2HC = INT(2,IA) PDF26500 IF( IA. EQ.NUM+2) IC=NUM-1 P0F2660O IW1 « I PDF26700 352 DO 355 IE-1,4 PDF26800 ID«IC-1*IE PDF2690O 355 XT(IEI«X(IOI PDF27000 DO AL 35fir= X l T :MO A ' IDER,l,MO,,MO,1 ' l ' KG ' MC,LD » XI ' GAM ' XHI pSfItISS IAB=IE4M0*( IB-ll PDF2730O IF(IW1 .EQ. II SAVE(IW2,IAB> » 0.0 oncflc?? IFdWl .EO. 1) HSAV(l.IAB) » 0.0 n ° ? HSAV<2,IAB) = 0.0 POF27600 DO 358 ID = lUlftM PDF27700 SAVE! ID*ltIABI»0. P0F2r800 DO 358 IF=l,MOLD PDF27900 DO 358 IG=l,4 PDF28000 IH-IC-1+IG PDF28100 IAA=IG+4*(IF-1» PDF28200 IAC=IF+M0LD*(IH-1» PDF28300 358 > D :^ IFIIC .ST. I .OR. IC .LT. NUM-1, GO TO 390 £™ ° IF(IC .EQ. 1» IC = 2 PDF29000 IFIIC .EQ. NUM-ll IC = IC -1 PDF29100 GO TO 352 PDF29200 : PDF29300 USING THE OLD PARTITION VALUES an™™ 370 DO 371 IJ = l t MO PDF29600 IAA = I J + M0*( IB-1 I PDF29700 IA3=IJ+M0LD*(IA-1» PDF29800 SAVE(IW*2,IAAJ=YMAX(IAB) PDF29900 HSAV<1,IAA| = HMYU.IABI PDF30000 HSAV(2,IAAJ = HMY(2,IA8) PDF30100 DO 371 10=1, IW PDF30200 SAVE< ID-H,IAA» = Y< ID.IAB) PDF30300 CONTINUE PDF30400 PDF30500 UPDATES THE VECTORS FOR THE NEW PARTITION PDFSO^OO J2=J*2 PDF3CB0O DO 410 IA=1,J2 PDF30900 N3=N2-l*JA+3*( IA-ll PDF31000 INT(2,IAI=NPSAVE(N3) P0F31100 IAB=1*M0MIA-1> PDF31200 X(IA»=SAVE{1,IAB) PDF31300 DO 410 IC=1,M0 PDF31400 IAA = IOM0*( IA-1 I PDF31500 YMAX(IAA»=SAVE( IW2.IAAJ oS«fJ2! DO 41C IBM.IW PDF31700 PDF31800 8o HMY(1,IAAI » HSAV(1,IAAI PDF3190( HMY(2,IAAI » HSAV(2,IAA) P0F3200C 410 Y(IB,IAAI=SAVE( IB+1 ,IAA) PDF3210( DO 430 IA * 1,4 PDF3220C 430 XT(IA) » XllA + ll POF3230C XI « X(l> PDF3240C CALL DERIV(XT,YTA,IDER(l,MOI,l,l,l,0,l,XI,GAM,XHI POF3250C DO 435 IA = ltMO PDF3260C DO 435 IB * 2,IW PDF3270C YUBfIA) =0.0 PDF3283C IF( IB .EQ. 21 HMY(2,IA) » 0.0 PDF3290C DO 435 IF * l f 4 PDF3300C IE = IA ♦ MO *( IF-1 I PDF3310C Y(IB,IAI = Y(IB,IA» + COEF< IF,1)*Y(IB,IEI PDF3320C IF( IB .EQ. 21 HMY(2,IA» ~ HMY(2,IA) ♦ COEF ( I F, 1 l*HMY ( 2, IE I PDF3330C 435 CONTINUE PDF33400 DO 440 IA * 1,4 PDF3350G 440 XT(IA) = X(J2 - IA> PDF3360C XI = XU2I PDF3370C CALL DERIV » 0.0 PDF3430Q DO 445 IF = 1,4 PDF34400 IE ■ IA ♦ MO * (J2 - IF -II PDF34500 YU8,IG> = Y(IB,IG) ♦ C0EF(IF,1» * Y(IB,IEI PDF34600 IF(IB .EQ. 21 HMY(2,IG) = HMY<2,IGI* COEF( IF, 1 )*HMY(2 , I El PDF34700 445 CONTINUE PDF34800 IDOUB=IH PDF34900 IQQQC=1 PDF35000 C PDF35100 C CALCULATES THE COEFFICIENTS FOR THE ERROR ESTIMATING INTERPOLATINGPDF35200 C POLYNOMIAL PDF35300 C PDF35400 258 DO 262 IA=1,J PDF35500 ID=MO*IA«-l PDF35600, XI=X(IA*1I PDF35700 DO 260 IB=1,5 PDF35800 IC=INT(2,IA + 1I-1+IB PDF35900 XT(IB>=X(ICI PDF36000 260 CONTINUE PDF36100 IAA=INT(1,IA*1 I* I A PDF36200 XH=XI~X(IAA<-1 I PDF36300 CALL DERIV(XT,YTA,IDER(1,M0|, MO, 1,0,1, MOLD, XI, GAM, XHI PDF364Q0 DO 262 IB=1,M0 PDF36500 IAA1 = 9*( IB-ll PDF36600 DO 262 IC=1,9 PDF36700 IAA = IAA1 ♦ IC PDF36800 262 DIVER( IAA,IDI=COEF( IC,IBI PDF36900 K7=KL8(M0I PDF37000 N=(NEQ(M0»*2I*M0 PDF37100 C PDF37200 C INTEGRATES AND DETERMINES THE MAGNITUDE OF THE ERROR POF37300 81 255 270 275 288 291 CALLDIFSUB(N,T,Y, SAVE, H,HMIN,HMAX, TPS, YMAX, ERROR, KFLAG.JSTART. 2hSq!S ICSTPS * ICSTPS ♦ 1 ICK » IQQQC=0 NUM=NEQ(M0)*M0 MOLD* MO IWJSTART + 1 IF(KFLAG .LE. -3) GO TO 290 IF(KFLAG.LE.O) RETURN NUM1 = (MUM*M0 NUM2=M0*1 RMAX = 0.0 DMAX = 0.0 00 27C IA=NUM2,NUM1 DBStA = DABS(E RORA(tAi) IF(DBSEA.GT.RMAX) RMAX = DBSEA IF{ ICPART .LT. 31 GO TO 270 DMAXT = DABSIEROR(IA) I IF(DMAXT .GT. DMAX) DMAX = DMAXT CONTINUE CHECKS ON ACTIONS BASED ON ERRORS AND COUNTERS IF(RMAX.GE.EPS) GO TO 290 ERROR IS SMALL ENOUGH SEPS = EPS IF(RMAX .GE. RVAL*EPS . OR. ( ICSTPS-I RECSP > .LT. 10) RETURN ICPART = ICPART ♦ I IRECSP = ICSTPS IF( ICPART .LE. 31 GO TO 275 ICPART = IFIDMAX .LT. SRVAL*EPS) GO TO 292 IF (ITUPS .GT. 3) GO TO 288 CHANGES TIME PARAMETER VTPS = DLOG ( TPS ) VTPS = ( VTPS ♦ VEPS I / 2. TTPS = DEXP (VTPS) TPS = DMIN1 (TTPS, TPS*. 5 DD* ( EPS-RMAX > ) ITUPS = ITUPS ♦ 1 RETURN IF(RMAX .GT. SRVAL*EPS) RETURN CHANGES SPACE PARAMETER IF(SX(MO)*2.0DC .GE. SMAX .AND. MO .EQ. MSTRT) RETURN ICK = 1 TSPS = 2. * SPS SPS = DMINKTSPS, .5 * EPS) PDF37400 PDF37500 PDF37600 PDF37700 PDF37800 PDF37900 PDF38000 PDF38100 PDF38200 PDF38300 PDF38400 PDF38500 PDF38600 P0F38700 PDF38800 PDF38900 PDF39000 PDF39100 P0F39200 PDF39300 PDF39400 PDF39500 PDF39600 PDF39700 PDF39800 PDF39900 PDF40000 PDF40100 PDF40200 PDF40300 PDFA0400 PDF40500 PDF40600 PDF40730 PDF40800 P0F40900 PDF41000 P0F41100 P0F41200 PDF41300 PDF41400 PDF41500 PDF41600 PDF4170Q PDF41800 PDF41900 PDF42000 PDF42100 PDF42200 PDF42300 PDF42400 PDF42500 PDF42600 PDF42 700 PDF428Q0 82 292 290 296 298 299 29 3 295 294 TPS » .6 * EPS -SPS IF( SX(M0I*2.000 .GE. SMAX .AND, ITUPS » I ICPART * NCNGE « I RETURN ERROR TOO LARGE ITUPS = I ICPART « ICSTPS = ICSTPS - 1 IRECSP = ICSTPS IF( ICSTPS .LT. 101 GO TO 298 IF(KFLAG .LE. -3» GO TO 293 IF(TPS .LT. .00001 .OR. KTPSCT < GO TO 294 KSPSCT = KSPSCT ♦ I KTPSCT »1 IF(KSPSCT .LE. 101 GO TO 293 KFLAG = -42 RETURN IF THE NUMBER OF STEPS IS SMALL ENOUGH PARAMETERS SPS = .1 * SPS TPS = .1 * TPS NCNGE = 1 ISPSCT = ISPSCT ♦ 1 IF( ISPSCT .LT.4 I GO TO 10 KFLAG = -35 IF(MAXSV .EQ. 1 J RETURN MAXSV = 1 IF( ICSTPS .GT. 10 I RETURN GO TO 2 DECREASES SPACE PARAMETER MO .EQ. MSTRTI RETURN GE. 5 I GO TO 296 - RESTARTS WITH SMALLER GO TO 295 IF(TPS .LT. .00001) TPS = .000001 J START = -I GO TO 255 SPS = PVAL * SPS GO TO 400 DECREASES TIME PARAMETER VTPS = DLOG(TPS) VTPS = (VTPS ♦ DL0G(i.0-6M/2.00 KTPSCT = KTPSCT ♦ 1 TPS =DEXP { VTPS) ITPSC = i JSTART = -1 GO TO 255 PDF42900 PDF43000 PDF43100 PDF43200 PDF43300 PDFA3400 PDF43500 PDF43600 PDF43700 PDF43800 PDF43900 PDF44000 PDF44100 PDF44200 PDF44300 PDF44400 PDF44500 PDF44600 P0F447C0 PDF44800 PDF44900 PDF45000 PDF45100 PDF45200 PDF45300 PDF45400 PDF45500 PDF45600 PDF45700 PDF45800 PDF45900 PDF46000 PDF46100 PDF46200 PDF46300 PDF46400 PDF46500 PDF46600 PDF4670Q PDF46300 PDF46900 PDF47000 PDF47100 PDF47200 PDF473G0 POF47400 PDF47500 PDF47600 PDF47700 PDF47800 PDF47900 PDF48000 PDF48100 PDF48200 PDF483D0 83 c OBTAINS THE OLD PARTITION c c 400 CALCULATED CONTINUE DO 405 KF1»1,N HMY(l,KFl» * HSAV ,NPTER (3 ,1 ) , A( 3, 2 , 5) 1 AFC 3) ,AD(3) ♦ GAM 1 3 ) , AVA ( 5 » , A A ( 4) , AE( 4) ,HQ( 21 »HSQ< 2 ) , IDER <4 , 3 ) , 2YT(3,5),XT(5),YTEMP(3,1),L<5) OIMENSION SX<3),NEQ(3) DIMENSION ALFA<2) , X ( 1 ) , Y< 5, MOLD, I » ,MAXM( 31 DIMENSION I WD 1(3,4) ,R WD( 3 ,2 ) , WVAR ( 3 ) COMMON /CNT/ ISTEP, NUMDE , NUMSO COMMON /DVR/ DVY( 4 > ,COEF ( 4, 3 , 31 , ALPHA ( 2 ,4 ,3 ) , DI V( 16, 5 J DATA IWD1 ,R WD/ 8, 9, 18, 18, 2 6, 2 8, 2, 2, 5, 5, 7, 8, 2., 6. ,10. ,3., 9. ,15./ DATA ALFA/1. DO, .9500/ DATA AVA/-1. ,-.50000000,0.0,. 250000000,1. / AT = ALFAUNPI RATIO = .10000000 IF(INP.NE. 2 .OR. NUMSO .LE.40) GO TO 10 RATIO = FLOAT(NUMDE)/FLOAT{ NUMSO) CONTINUE DHMIN=HMIN/2. K*l CHECK INPUT PARAMETERS IFIMSTRT .GT. MAXO .OR. MSTRT .LE. 0) GO TO 2000 MO = MSTRT RNT=RINT(2)-RINT(1) IF(RNT .LE. O.DO) GO TO 2010 IF(MAXO .GT. 3 .OR. MAXO .LT. 1) GO TO 2300 INITILIZE PARTITION VALUES FOR NEW ORDER MAXN=MAXM( DO 30 1=1, NCK( I )=-l H=RNT/2. DP0S(M0,2, H=H/2. N=3 SX(MO)=H NEQ(M0)=3 NPT(l,l)=2 NPT(1,2)=1 NPT(1,3I=3 DP0S(M0,2, DP0S(M0,2, DP0S(M0,1, DP0S(M0,2, DP0S(M0,1, DP0S(M0,3, DP0S(M0,1, DP0S(M0,2, DP0S(M0,3, DP0S(M0,3, IF( INP .NE MO) MAXN 1)=RINT(1)+H MAXN + MAXN* 2) = RI 2)=RI 1 ) = DP 2)=DP 3)=DP 3) = RI 3) = RI 1) = DP • 1) 1)=RINT(1 I 2)=RINT(2) NT(1 ) NTI1 )+H OS(MO,2,2) 0S(M0,2,1) 0S(M0,2,1) NT(2)-H NT (2 ) OS(MO,2,3) GO TO 1800 STP139( STP140( STP14K ,STP142( STP143C STP1440I STP1450( STP1460( STP1470(j STP14800 STP1490C STP1500C STP1510C STP1520G STP15300 STP1540C STP1550C STP1560C STP1570G STP1580C STP15900 STP16000 STP1610Q STP1620C STP1630C STP16400 STP16500 STP16600 STP16700 STP1680C STP16900 STP17000 STP17100 STP17200 STP17300 STP17400 STP17500 STP1760Q STP17700 STP178D0 STP17900 STP18000 STP13100 STP18200 STP18300 STP18400 STP18500 STP18600 STP18700 STP13800 STP18900 STP 19000 STP19100 STP19200 STP19300 85 INITILIZES ITEMS FOR THIS PASS OF INTERVAL 100 110 C c c 120 130 210 H=H/2. ILD»1 INW»2 JLD*1 JNW=*1 NMCK=0 HH=H HHP=HH*. 941 59265 HQ(1)*HH HQ(2I=HHP HSQ(1)=HH*HH HSQ(2»=HHP*HHP SETS UP THE VALUES FOR ADDITIONAL POINTS I=NPT( ILD,JLD) IFINCKII I.GT.Ol GO TO 600 FWD=DP0S(M0,2f I )+H BWD=DP0S(M0,2,I »-H CHECKS TO SEE IF THE POINT IS ALREADY IN THE PARTITION KA*-1 DO 130 J=1,N IF = NPT( INW,JNW-3I STP2500d 1(51 = NPT(ILD,JLD+1 ) STP25100 IF( JLD.GE.NEQ(MO) I L(5l=MAXN+2 STP25200 IF1JLD.LE.2I L(1I*MAXN+1 STP25300 IF(JLD .GT. 1) GO TO 400 STP2540Q L(2) = MAXN + 1 STP2550(J L(1) = NPT< ILD,JLDU) STP25600J l(5l=NPT( ILD,JLD*2I STP25700 GO TO 400 STP25800 C STP25900 C CHECKS THE PRESENT PARTITION POINT AND PUTS POINT TO RIGHT STP26000I C INTO THE PARTITION STP26100: C STP26200 220 NPT(INW,JNWI=I STP26300 JNW = JNW-H STP26400 REC=DP0S(M0,3,I» STP26500 DP0S(M0,lfI )=BWD STP26600 0P0S{M0,3tI )=FWO STP26700 IF(N+1.GT.MAXNI GO TO 2030 STP26800 N=N*1 STP26900 DP0S(M0,2 f N»=FW0 STP27000 NPTUNW,JNWl=N STP27100 JNW=JNW+1 STP27200 DP0S(M0,l,NI=DP0S(M0,2tII STP27 300 DP0S(M0,3tN)=REC STP27400 JC=N STP27500 JA = 3 STP27600 ERAR = EROR STP27700 IFUNP .EQ. II GO TO 300 STP27800 KLA=5 STP27900 XI=DP0S(M0,2,NI STP28000 GO TO 1850 STP28100 230 L(1)=NPT( INW,JNW-3I STP28200 L(2)=NPT( INW.JNW-21 STP28300 L(4I=NPT< ILD,JLD+1> STP284CG 1 L(5) = NPTULD t JLD+2) STP28500 IF( JLD.GE.NEQ(M0)-1 I l(5l=MAXN*2 STP28600' IF( JLD.LT.NEQ = NPT( ILD,JLD-HI STP29900 L (2>=NPT( INW.JNW-31 STP3D000 L(4»=NPT( INW f JNW-ll STP30100 L(l » = NPT( INW,JNW-4) STP30200 IFULD.EQ.ll L(1I = MAXN+1 STP30300 87 IF(JLD.GE.NEQ*HH STP31000 IF(DABS(V1-RINT(1H.LT. 1.1*HH .OR. DABS ( Vl-RINT (2)1 . LT. STP3110Q 1 1.1*HHJ ERAR = .5 * ER3R zll.Z.iz.Z. STP31200 1.1*HH) ERAR ■ .5 * EROR CALL INIT(A(1,1,JB),M0,V1J STP^lino 310 CALL INIT(A(1,2,JB>,M0,DP0S(M0,2,JC>*AVA(JB»*HHP» STP31400 GO TO (320, 350, 380), MO STP31500 C EVALUATES USING LAGRANGIAN - INITIAL STP31700 C320 JB '° STP31900 1 330 JB . JB ♦ I iTP320oS IF(IDER(1,1I.EQ.0) GO TO 340 STP3710O AA(ll = (A(l,JB,5l-A(l,JB,in/HQ(JB» STP 32200 ,J C r!f , /!inl , ,o? , ^^^? < ? J0Q0 STP^nn 340 IF(IDER(2,1>.EQ.0> GO TO 510 STP32600 AA(2I = (A<1,JB,5|-2.*AU,JB,3)*A(1,JB,1U/HSQ( JB » STP 32 700 AC=( 10. 666666666666667' All, J 8,5 I -136. 53333333 3333* A (1,JB, 4 » +192. 00 STP32800 1000*A(1,JB,3J-85.3333333333333*A<1,JB,2)*19.200000*A{1,JB,1I)/*28STP34400 280.0JCC*A(2,JB,3l*2560.00O0*A(2,JB,2N /HQ(JB» STP34500 AE<2>=AA(2)*AC/32 0. 00000 stpiZa™ 37C IF(IDER(3,2I.E0.CI GO TO 510 STP34700 /A<3j=<(A(l t JB.5J-A(l,JB,in*7.5000000/HQ "242 6. 66666666667* A(1,JB,5> I /HQ< JB » * ( 5040. O0OC*A (2 , JB , 1 1 +3584STP35100 10.000*A(2,JB,2I*20160.000*A(2,JB,3I*560.00000*A(2,JB,5M>/HSO(JB» STP35200 rnVn ^a** 3 ' * AC/ 84O.OOO000 STP35300 c G ° T ° 51 ° STP35400 C EVALUATES USING TRIMITE - INITIAL STP35600 380 JB - O STP35700 380 JB " ° STP35800 C 88 390 JB » JB ♦ 1 STP3590 AA(1)=A(2,JB,3I STP3600 AE(1I=AA(1I STP3610 AA(2I=A(3,JB,3I STP3620 AE(2I*AA(2I STP3630 IF( IDER(3,3).EQ.O) GO TO 395 STP3640 AA(3»»ll3.1250L0/HSQ(JB)*CA(l,JB,5»-A(l.JB,l) 1-1. /HQ( JB I * (4. 125000STP3650 10*A(2,JB,5»+18.00 0000*A(2,JB,3H-4.1250000*A(2,JB,1) K. 375 00000* ( A ( STP3660 23t JBt 5I-A(3 f JBtim/HQ( JB) STP3670 AC=( (6441 12 00. *A( 1, JB,H- 263782400. *A( I, JB, 2) +199584000. *A(1,JB, 3) STP3680 l-212 8uC.DO*A(l,JB,5 I » /HSQ( J B » ♦( 1065960C. * A( 2, JB, 1 1-344C6400. *A( 2, JSTP3690 2B, 2 ) -43545600. *A( 2, JB, 3 J + 25200. 000* A( 2 , JB ,5 I » /HQ( JB )♦ ( 498960. 00* A ( STP3700 33, J B,l) -9461 760, 0*A ( 3, JB, 2) +3265920. 0*A( 3 , JB, 3H-1680. 0000*A (3 , JB, 5STP3710 4)>>/HQ(JB) STP3720 AE( 3»=AA(3»+AC/60480.000 STP3730 395 IF( IDER(4,3».EQ.0» GO TO 510 STP3740 AA(4l=(72.000000/HSQ(JB)*(A(l,JB,5>-2.*AU,JB,3>+A(l,JB,l»)-19.500STP3750 10000/HQ(JBI STP3760 1*(A(2,JB,5)-A(2,JB,1)»+1. 5000000* (A(3,JB, 51-24. OOOOOQ*A( 3 ,JB, 3 J+A( STP3770 33,JB,1 m/HSQ(JB> STP3780 AC=(( 1146880000. *A( 1, JB , 2 »-l 5240960C. *A( 1 , JB, 1 » -1001548800, *A ( 1 , J BSTP3790 1,3»+7078400.0*A(1,J8,5I J/HSQl JB >+( 206438400.* A( 2, J B, 2 1-22680000. *ASTP380Q 3 ( 2, JB, 11+2 32243200. *A(2 , J B, 3l-204963D.0*A< 2, JB,5 I » /HQ( JB )+ ( 344C64CSTP3810 30.*A(3,JB,2»-907200.00*A(3,JB,1>-21772800.*A(3,JB,3»+168000.00*A(3STP3820 4,JB,5)n/HSU(JB) STP3830 AE(M=AA(4)+AC/151200.00 STP3840 GO TO 510 STP3850 C STP3860 C PREPARE FOR EVALUATION OF THE DERIVATIVE STP3870 C STP3880I 400 L(3»=JC STP3890 DO 410 IB=1,5 STP3900i LA=L(IB) STP3910 XT( IB)=DP0S(M0,2,LAI STP3920 DO 41C IA=1,M0 STP3930 410 YT( IA,IB»=YTEMP(IA,LA» STP3940 CALL DERIV(XT,YT, IDERd ,M0) , MO, 0, ,D , MOLD,X I , GAM, HQ( 1 >> STP3950 GO TO (420,440,4601 , MO STP3960 C STP3970' C EVALUATES USING LAGRANGIAN - MODIFYING STP3980 C STP3990 420 IF( IDERd, D.EQ.OI GO TO 430 STP4000 A A( 1 I = (YT(1,4)-YT( 1,2)1/(2. 0*HQ(1)> STP4010' AE( 1)=AA(1)+DVY<1 I STP402CI 430 IF( IDER(2,l).EQ.O) GO TO 500 STP4030I AA(2»=(YT(1 ,4)-2.*YT(l, 3)+YT(l ,2))/HSQ(l) STP4040( AE(2J=AA(2)+DVY(2 » STP4050I GO TO 500 STP4060 C STP4070I C EVALUATES USING HERMITE - MODIFYING STP4080I C STP4D90( 440 AA(1)=YT(2,3» STP4100I AE(1)=AA(1> STP4U0 IF( IDER(2,2I.EQ.0I GO TO 450 STP4120* AA(2) = 2.*(YT(1 ,4I-2.*YT(1 ,3)+YT(l ,2 )) /HSQ ( 1 »-( YT( 2 ,4I-YT(2 , STP4130I 89 450 460 470 500 510 520 530 540 1211/(2. 0*HQ(1N AE(2»=AA(2»+DVY(2) IF( IDER(3,2).EQ.O) GO TO 500 l"i?( , 2 , :3IUT;2 , ;II!}/Hsi:n• 50oooooo/Ha(l, - l ■ 5o3oool> • |YT,2 AE(31=AA(3»*DVY(3I GO TO 500 EVALUATES USING TRIMITE - MODIFYING AA(1)=YT(2,3I AE(1I*AA(1) AA(2)=YT(3,3I AE(2)=AA(2> IF< IDER(3,3).EQ.OJ GO TO 470 AA(3» =( 13.1250uC/HS0(l»«-{ YT(1 2YT(2,2n+18.CC0000*YT(2,3M/ 1HQ(1»+. 37500000*1 YT(3 ,4>-YT(3 AE(3>=AA(3I*DVY(3» IF(IDER(4,3».EQ.0I GO TO 500 AA(4>=(72.0000GO/HSQ(l»*(YT(l ,4)-2.*YT(l 1HQ(1>»(YT(2 ,4)-YT(2 ,2 ) I *1. 5000000* ( YT( 3 2TC3 ,2MI/HSQ(ll AE(4)=AA(4»+DVY(4I MAKES USE OF THE USER'S EQUATION CALL UFUN(M0,DP0S(M0,2,JC),YT(1,3»,T,AA,AD,1I CALL UFUN(M0,0P0SCM0,2,JCI,YTC1,3J,T,AE,AF,1I GO TO 520 CALL UFUN(M0,DP0S(M0,2,JC»,A(1, JB,3>,0. ,AA,AD,1I CALL UFUN(MO,DPOS(MO,2,JC),A(l,JB,3»,0.,AE,AF,l| CHECKS THE ERROR DO 540 IL =* 1,M0 AF( IL) = AF{ IL I - ADCILI IF(0ABS(AD( ID » .GT. 1.D0I AD( IL) = AD(ID/AD( ID IF(DABS(AF(IL») .GT. ERARI CONTINUE IF{ INP .EQ. 2 .OR. JB .GE. 2) GO TO (330,360,390), MO ,4>-YT(l ,2))-(4.1250000*(YT(2 ,2>))/HQ(l> ,3»+YT(l ,211-19.5 ,4)-24. 000000*YT(3 GO TO 533 GO TO 560 GO TO 550 !??ini TES ™ AT ™ E ERR0R IS ALL RIGHT, SO THAT THE POINT IS IGNORED ON FUTURE PASSES 550 560 NCK(JC) » 1 NMCK = NMCK ♦ 1 GO TO (220,610, 2401, JA CHECKS FOR COMPLETION OF PASSES. FOR NEW PASS 600 NPT(INW,JNWI»I IF NOT DONE, PREPARES STP41400 STP41500 STP41600 ,41+8 STP41700 STP41800 STP41900 STP42000 STP42100 STP42200 STP42300 STP42400 STP42500 STP42600 STP42700 STP42800 »4) + STP42900 STP43000 STP43100 STP43200 STP43300 000GD/STP4340O ,3I*YSTP43500 STP43600 STP43700 STP43800 STP43900 STP44000 STP44100 STP44200 STP44300 STP44400 STP44500 STP44600 STP44700 STP44800 STP44900 STP45000 STP45100 STP45200 STP45300 STP45400 STP45500 STP456T)0 STP45700 STP45800 STP45900 STP46000 STP46100 STP46200 STP46 300 STP46400 STP46500 STP46600 STP46700 STP46800 90 JNW*JNW*1 NMCK=NMCK*1 610 IF( JLD.GE.NEQ(MO) I GOTO 620 JLD=JLD+1 GO TO 123 620 IF(NMCK .GE. Nl GO TO 790 ILD=3-IL0 INW=3-INW NEQ(MO>»N H=H/2. IF(H .LT. HMIN) GO TO 2090 SXl MO)»H GO TO 110 : C PREPARE THE WORK ESTIMATE : 700 IF(NEQ(MOI .GT. 1 GO TO 710 WVARJMOI = 1.0050 GO TO 720 710 WVTP = IDERIM0,M0)*IWD1(M0,U ♦ IDER( M0+1,M0I *IWD1 (M0,2) TMO = 2*M0 WVTP = WVTP/TMO WVTR = IDERJM0,M0I*IWD1(M0,3> ♦ I OER ( MO + 1 ,M0 ) *I WOK MO f 4 1 RN = N WVAR{MO> = RN*(WVTP ♦ RAT 10* i WVTR*RN+RWD( M0,1 » ) ♦ RWDI MO = 3 74-0 MONW = MO IF(WVAR(MO» .GE. 1.0D45) GO TO 2040 RETURN CHECKS THE PARTITION TO INSURE SMALLEST AVAILABLE STEP SIZE IS USED 790 JNW = 1 NEQ(MO)=N 800 JA=NPT( INW,JNW» IF< JNW.EQ.l) GO TO 805 IF( JNW.NE.NI GO TO 810 H=RINT(2»-DPOS(MO,2,JAI JLD=-1 GO TO 830 805 H=DP0S(M0,2,JAI-RINT(1I JLD=1 STP4690I STP4700< STP4710I STP4720I STP4730< STP4740I STP4750< STP4760( STP4770C STP4780< STP4790( STP4800C STP4810C STP4820C STP4830C STP4840C STP4850C STP4860C STP4870C STP4880C STP4890C STP4900C STP4910C STP4920C STP4930C STP4940Q STP49500 STP4960C STP49700 STP4980Q STP49900 STP50000 STP50100 STP5D200 STP5O3O0 STP50400 STP50500 STP5C60C' STP5070C STP538D0 STP50900 STP51000 STP51100 STP51200 STP51300 STP51400 STP51500 STP51600 STP51700 STP51800 STP51900 STP52000 STP52100 STP52200 STP52300 91 810 820 830 840 850 860 GO TO 830 JB»NPT(INW,JNW-1» JC=NPT< INW,JNW*1) JLD = 1 VA=DPOS(MO,2,JA)-DPOS XT( I11)=X(I12I CALL DERIV(XT,YT,IDER(1,M0),M0,0,1, 0. MOLD, XI, GAM, HQ(1 II DO 1920 113 = 1,M0 YTEMP(I13,JC)»GAM( 1 13 I GO TO 1950 DO 1940 113 * l f MO YTEMP(U3,JC»*Y(1,I13,N3) GO TO (1820, 1830,100, 210,230, 1840, 18101 *KLA K » -201 RETURN K » -202 RETURN K = -205 RETURN RETURN COOES ARE GENERATED K = -206 NEQ(MOI»-3 GO TO 700 K * -209 RETURN K = -208 RETURN K = -207 RETURN K = -203 RETURN K = -204 RETURN NEQ(MO) = -213 GO TO 7C0 END STP57900 STP58000 STP58100 STP58200 STP58300 STP58400 STP58500 STP58600 STP58700 STP58800 STP58900 STP59000 STP59100 STP59200 STP59300 STP59400 STP59500 STP59600 STP59700 STP59800 STP59900 STP60000 STP60100 STP60200 STP6D300 STP60400 STP60 500 STP60600 STP60700 STP60800 STP60900 STP61000 STP61100 STP61200 STP61300 STP61400 STP61500 STP61600 STP61700 STP61800 STP61900 STP62000 STP62100 STP62200 STP62300 STP62400 STP62 500 STP62600 93 c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c SUBROUTINE DERIVU, FN, IDER,MO,KB, KF,KG, MOLD, XI. GAM, HI THIS PROGRAM IS USED TO FORM INTERPOLATION POLYNnMTAK cno LAGRANGIAN, HERMITE, AND TRIMITE INTERPOUtJS POLWOHuS FOR FUNCTIONAL VALUES AS WILL AS SEVERAL DERIVATIVES li IS DIVIDED INTO TWO SECTIONS. ONE SECTION CALCULATES ill INTERPOLATING POLYNOMIALS NEEDED IN EST MATING ERROR I THf OTHER SECTION IS USED TO CALCULATE NEW VALUES FOR fSnCTiJnI DESlRll 0ClAJED DERIVATIVES WHEN NEW PARTITION POINTS ARE THE VARIABLES X Y IDER MO MOLD ALPHA DIV COEF KB IN THIS SUBPROGRAM ARE CONTAINS INTERPOLATING POINTS CONTAINS FUNCTION AND DERIVATIVE VALUES TELLS WHICH DERIVATIVES ARE REQUIRED ORDER OF METHOD (POSSIBLY NEW) OLD MO CONTAINS FOR EACH CONTAINS (1,1) (2,11 (3,1) C4,IJ (5,1) (6, I) (7,1) (8,11 (9, II (10,1 I (11,1 I (12,1) (13,1) (14,1 ) (15,1) (16,1 ) CONTAINS THE FACTORS X-A(I) AND A(J)-A(II POINT AND COEFFICIENT VARIOUS FUNCTIONAL VALUES FOR EACH POINT L"MXI LMX) LMIMX(II) L'MXI L'MIIIXIIII DENOMINATOR OF L(X) A* (A BAR)* (A DOUBLE BAR)* A* • (A BAR)* • DOUBLE BAR)** (A A A A BAR DOUBLE BAR L(X) THE COEFFICIENTS THOSE FROM 7 TO 15 ARE ALSO USED TO STORE TEMPORARY RESULTS IN INTERPOLATION OF THE RESULTING OTHERWISE KG KF XI GAM DVY H 1 IF COEFFICIENTS DESIRED, INTERPOLATING POLYNOMIALS 1-INTERPOLATE 0-ERROR ESTIMATE 1-AT PARTITION POINT 0-NOT AT EVALUATION POINT RESULTS OF INTERPOLATION RESULTS OF ERROR ESTIMATES DISTANCE USED IN ANALOF. HSQ,HSQQ ,&HSQC rur occr nc -,IDER(4|,GAM(3I HSQ = - H ♦ H ARE DRV10000 DRV10100 DRV10200 DRV10300 DRV104OO DRV10500 DRV10600 DRV10700 DRV10800 DRV1C930 DRV11000 DRV11100 DRV11200 DRV11300 DRV11400 DRV11500 DRV11600 DRV11700 DRV1180O DRV11900 DRV12000 DRV12100 DRV12200 DRV12300 DRV12400 DRV12500 DRV12600 DRV12700 DRV12800 DRV12900 DRV13000 DRV13100 DRV1320O DRV1330O DRV13400 DRV13500 DRV13600 DRV13700 DRV13800 DRV13900 DRV14000 DRV14100 DRV14200 DRV14300 DRV14400 DRV14500 DRV14600 DRV14700 DRV14800 DRVK900 DRV15000 DRV15100 DRV15200 9k 10 c c c c c c 20 C c c 30 40 C C C 45 50 C C C IF(M0 .LT. 2 » GO TO 5 HSOQ*HSQ*HSQ IFCMO .LT, 3 I GO TO 5 HSQC»HSQQ*HSQ TEST IF DERIVATIVES ARE USED KA * DO 10 IA»1,4 DVY( IA»=0. KA=KA*IDER{ IA I IF(KA. EQ.O ) RETURN SETS FOR HIGHEST LAGRANGIAN DERIVATIVE IF(M0.NE.1.0R.IDER<2).NE.l) GO TO 20 Nl»5 K=2 VA*12. V=24. GO TO 30 SETS FOR THE LOWER LAGRANGIAN DERIVATIVE IF(M0.EQ.1.AND.IDERU).EQ. 0) RETURN K»l Nl=4 VA=6. V*6. OBTAINS DENOMINATOR OF LAGRANGIAN POLYNOMIALS DO 40 IA=1,N1 DIV(6,IA»-1. DO 40 IB=1,N1 IF{ IA.EQ. 181 GO TO 40 DIV(6,IA»=DIV(6,IAI*(X(IAI-X(IBI» CONTINUE IF(KF.NE.l) GO TO 45 KSAVE=KA GO TO 70 THE FOLLOWING FORMS L^'MXI DO 50 IA=1,N1 DIV(1,IAJ=V/DIV(6,IA) IF(M0,NE.1I GO TO 70 CALCULATES ERROR ESTIMATES FOR LAGRANGIAN DVY(K>*0. KA = 4*(K-1) DO 60 IA*1,N1 IF(KB .EO. I GO TO 58 IB = KA ♦ IA DRV1530I DRV1540I DRV1550I DRV1560! DRV1570* DRV1580< DRV1590<; DRV1600I DRV1610< DRV16204 DRV1630I DRV1640 GO TO 240 KA=1 GO TO 90 CALCULATES A FOR TRIMITE DO 230 IA=1,N1 DIV(15,IA)=.5*(XI-X(IAH*{XI-X(IAII DIV(13,IA)=1.+DIV(15,IA)*DIV(10,IAI-3.0000000*(XI-X(IA|I It IAI DRV2630I DRV2 640 DRV2650 DRV2660I DRV2670 DRV2680I DRV2690I DRV2700I DRV2710J DRV2720I ORV27301 DRV27401 DRV2750t DRV2760( DRV2770( 0RV2780( DRV2790( DRV2800( DRV2810I DRV2820( DRV2830( DRV2840( DRV2850( DRV2860( DRV2870C DRV2880( DRV2890( DRV2900C DRV2910C 0RV2920( DRV2930C DRV2940C DRV2950C DRV2960( DRV2970C DRV2980C DRV2990(; DRV3000C DRV3010C DRV3020C ORV3030C DRV3040C DRV3050C 0RV3060C DRV3070C DRV30800 DRV3090C DRV31000 DRV3U0G DRV31200 ORV31300 DRV31400 DRV31500 ♦DIVC30RV31600 DRV31700 97 DIV(14,IA|»(XI -XUAII ♦DIV(15,IA)*DIV(11.IAI IF{KF.EQ.1> GO TO 230 ' J OIVlll.IAl DRV31800 DIV(6,1)«DIV(1,IA|*DIV(1,IA| DRV31900 C DRV32000 C CALCULATES THIRD DERIVATIVE FOR TRIMITE 0RV32100 V*l 680.0000*0 1 V( 1,1 A)*DIV( 6, 11 DRV32300 VA=15120.000*DIV(4,IAI*DIV(6,1| DRV32400 JB.226e0.0OO.DIVU f U,,D,vu fIM , OIV(4fI „^^ DO 225 KFL =1,3 DRV33100 2 " ^'^•< F '-»^ = VA*DIV(9,IA)*VB*DIV I A S!X«?S2 IF(KB.EO.OJ GO TO 250 ' ' DRV35100 DO 245 KFL=1,3 DRV35200 245 j° E ;s i ;j;^ L «w«".2tMLi/isi»».oo *hsqc ski*™ 251 CONTINUE DRV35700 IF ( KB .EQ.1» RETURN DRV35800 DVYU)=DVY(4,/151200.00 *HSQC DRV35900 RETURN DRV36000 260 IF(MOLD .GE. MO .AND. KG .EQ.l. GO TO 520 DRV362^ Nl«* DRV36300 GO TO 30 DRV36400 270 KA=1 DRV36500 GO TO (290, 90, 1201, MO DRV36600 C DRV36700 C CALCULATES LUI DRV36800 C DRV36900 280 IF(KA.E0.2| GO TO 270 DRV37000 290 DO 300 IA=1,N1 DRV37100 DRV37200 98 300 310 320 330 3*0 350 360 370 380 390 400 C c c DIV(16,IA)»ALPHA(1, I A , 1) * Al PHA{ 1 , I A , 2 )* ALPHA( 1, I A , 3 I /DI V( 6, I A ) GO TO (310, 380, 200), MOLD GO TO (360,340,320) ,M0 INTERPOLATES THE SECONO DERIVATIVE USING LAGRANGIAN GAM(3)«0. DO 330 IA*1,N1 COEF( IA,1,3)=DIV(4,IA) IF(K8 .EQ. 1) GO TO 330 GAM(3) = GAM(3H-DIV(4,IA)*FN(l,U) CONTINUE INTERPOLATES THE FIRST DERIVATIVE USING LAGRANGIAN GAM (2 ) = 0. DO 350 IA=1,N1 COEF(IA,l,2)=DIV(2,IA I IF(KB .EQ. 1) GO TO 350 GAM(2)» GAM(2) +D I V(2 , I A)*FN( 1, IA ) CONTINUE INTERPOLATES THE FUNCTIONAL VALUE USING LAGRANGIAN GAM(l) * 0. DO 370 IA=1,N1 COEF( IA,1,1)=DIV( 16, IA) IF( KB .EQ. 1 ) GO TO 370 GAM(l) = GAM(l) ♦ DIV(16,IA) * FN(1,IA) CONTINUE RETURN DEVELOPS COEFFICIENTS *OIV( 7,IA) DO 390 IA=1,N1 DIV(7 ,IA)=-2.»DIV(3,IA) DIV(13, IA) = 1.+(XI-X(IA) ) DIV(8,IA)=1. DIV(14,IA)=XI-X(IAI GO TO (43C, 420,400, MO INTERPOLATES SECOND DERIVATIVE USING HERMITE 410 GAM DO V = 2 VA = COE COE IF GAM 1 IA) CON (3) = 0. 410 I •*(DIV 4.*DIV F( IA,1 F( IA,2 (KG .E (3)=GA TINUE A=1,N1 (2,IA>*DIV(2,IA)+DIV(16,IA)*DIV(4,IA>) (16,IA)*DIV(2, IA) ,3)=V*DIV(13,IA)+VA*DIV(7,IA> ,3) = V*DIV( 14,1 AKVA Q. 1) GO TO 410 M(3) +C0EF(IA,i,3)*FN(l,IA)+C0EF(IA,2,3)*FN( INTERPOLATES FIRST DERIVATIVE USING HERMITE DRV3730( DRV3740( DRV3750( 0RV3760(j DRV3770C DRV3780d DRV379O0 DRV38000 DRV3810(J DRV38200 ORV38300 DRV38400 DRV38500 ORV3860C DRV 38700 DRV38800 DRV38900 DRV39000 DRV39100 DRV39200 DRV39300 DRV39400 DRV39500 DRV39600 DRV39700 DRV39830 DRV39900 DRV40000 DRV40100 DRV4020O DRV40300 DRV40400 DRV40500 DRV40600| DRV40700 DRV40800;! DRV4090Q DRV41000, DRV41100 DRV41200 DRV41300 DRV41400 DRV41500 DRV416Q0 DRV41700 DRV41800 DRV41900 DRV42000 DRV42100 DRV42200 2,DRV42300 DRV42400 DRV42500 DRV42600 DRV42700 99 420 GAM(2I«0. DO 425 IA»l,Nl V*2.*DIV(16,IAI*DIV<2,IAI VA«DIV(16,IAI»DIV(16,IA» C0EF(IA,1,2»*V*DIV(13,IAH-VA*DIV<7.IAI C0EF(IA,2,2»xV*DIV *FN(l,IA) 425 ♦ C0EF(IA,2,2»*FN(2,IA| 43 440 450 C C C INTERPOLATES THE FUNCTIONAL VALUE USING HERMITE GAM(1J=0. DO 440 IA"1,N1 VA=DIV(16,IA)*DIV(16,IA» COEFt IA,i,l|=VA*DIV(13, IA I COEF( IA,2,1 J=VA *DIV(14,IAI IF(K8 . EQ. 1 I GO TO 44C CONTINUE GAM*FN(2,IA) 460 470 INTERPOLATES SECOND DERIVATIVE USING TRIMITE GAM(3)=0. DO 470 IA*1,N1 VB=DIV(16,IAJ*DIV(16,IAI VA=6. 0000000* VB*D IV (2, 1 A) ve=^im:TA? ,V,4 ' U,t6 - Coooaoc * 0,v(l6 ' IA, ' OIV '^'*»*°'v, 2 ..A, cnlt "•!!•?!"'•<>"" ".'A. ♦vA.ome A MHo V 11 A S«8"M:n , "S To l "i u,w " wy " ,,u, * vwoM, » : »» CONTINUE ♦ C0EF(IA,2,3)*FN(2, IAI 480 490 INTERPOLATES FIRST DERIVATIVE USING TRIMITE GAM(2>=0. DO 490 IA=1,N1 VB=DIV(16,IAI*DIV(16,IA» V=3. 0000000* VB*DI V(2,IAI VA=0IV(16,IA)*VB COEF( IA,1,2I = V»DIV(13,IAH-VA*DIV(7,IAI C0EF( IA,2,2, = V*DlV(l4,IAH-VA*DIV<8,IA) rc^i IA : 3,2,=V¥DIV(15 * IA ' +VA *0^(9,IA) IF(KB .EQ. 1 ) GO TO 490 GAM(2> = GAM(2I ♦ C OEF ( I A, 1 , 2 )*FN( 1 , I A) I ♦ C0EF(IA,3,2)*FN(3,IAI CONTINUE COEF(IA,2,2l*FN(2,IAI DRV42800 DRV42900 DRV43000 DRV43100 DRV43200 DRV43300 DRV43400 DRV43500 ORV43600 DRV43700 DRV43800 DRV43900 DRV44000 0RV44100 DRV44200 DRV44300 DRV44400 DRV44500 DRV44600 DRV44700 DRV44800 DRV44900 DRV45000 DRV45100 DRV45200 DRV45300 0RV45400 DRV45500 DRV45600 DRV45700 DRV45800 0RV45900 0RV46000 DRV46100 DRV46200 DRV46300 DRV46400 DRV46500 DRV46600 DRV46700 DRV46800 DRV46900 DRV47000 DRV47100 DRV47200 DRV47300 DRV47400 DRV47500 DRV47600 DRV47700 DRV47800 DRV47900 DRV48000 0RV48100 DRV48200 l.OQ INTERPOLATES THE FUNCTIONAL VALUE USING TRIMITE 5CO GAM(1I»0. DO 510 IA»1,N1 V=0IV(16,IA>*DIV(16,IA>*DIV(16,IA) COEF{ IA,l,ll=V*DIV(13,IAI COEF( IA,2,1I = V*DIV( 14.IAI COEF{ IA,3 t l )=V*DIV(15,IA» IFCKB .EQ. 1) GO TO 510 GAM(l) = GAM(l) + COEF(IA,i,l)*FN(l,IAI 1 ♦ COEF< IAt3tll*FN(3t IAI 510 CONTINUE RETURN 520 WRITE(6,1C?01) 10001 FORMAT!* INVAL ID ORDERS • I RETURN END ♦ C0EF(IA,2,1H"FN<2, I A) DRV4830C DRV4840C DRV4850C DRV4860C DRV4870C DRV4880C DRV4890C DRV4900C DRV4910C DRV^920C DRV4930C DRV49400 DRV4950C DRV49600 DRV49700 DRV49800 DRV49900 C c c c c c c c c c c c c c c SUBROUTINE DECOMP *UL(L,K)-UL(L1,UI*JLCJ,KI IF(K2.EQ.l) RETURN DO 20 I»2,K2 11 » N ♦ I K5 » U - K3 12 * K2 - I DO 20 J«I,K2 K * U - J IF(DABS(UL(K3,K5M.LT.l.D-20» GO TO 15 UL(J,K>=UL(J f K)/UL(K3,K5) J1»J+1 J2 » 12 ♦ Jl LI * K3 DO 20 L»J1,J2 LI » LI ♦ 1 UL(L,K»=UL(L,K»-UL(L1,K5I*UL(J,KI J5=-l RETURN END DCP13300 DCP13400 DCP13500 DCP13600 DCP13700 DCP13800 DCP13900 DCP14000 DCP14100 DCP14200 DCP14300 DCP14400 DCP14500 0CP146C0 DCP14700 DCP14800 DCP14900 DCP15000 DCP15100 DCP15200 DCP15300 DCP15400 DCP15500 DCP15600 OCP15700 DCP15800 DCP15900 DCP16000 DCP16100 DCP16200 DCP16300 DCP16400 DCP16500 SUBROUTINE SOLVEJ N.UL, B, X, K7, ID,K> IMPLICIT REAL *8(A-H,0-ZI THIS PROGRAM USES THE DECOMPOSED MATRIX FROM DECQMP Tn PERFORM THE GAUSSIAN ELIMINATION T ° THE VARIABLES UL N K7 K3 K2 B X ID K IN THIS SUBPROGRAM ARE THE ARRAY THE NUMBER OF VARIABLES THE NUMBER OF DIAGONALS IN THE MATRIX iS S £5?JIKs THE MA,N 0UG0N4L thI "suLT N J 6 c V rm° R ° F ™ E E0UAr,0N I0 BE S0L ' E0 THE NUMBER OF ROWS IN THE RESULT VECTOR : T uL?? W IN RESULT VECT0R WHERE AN SWER IS PLACED THE REST OF THE VARIABLES ACT a"s COUNTERS OR A S PINTERS SLV10000 SLV10100 SLV10200 SLV10300 SLV10400 SLV10500 SLV10600 SLV10700 SLV10800 SLV10900 SLV11000 SLV11100 SLV11200 SLV11300 SLV11400 SLV11500 SLV11600 102. . TO POSITIONS IN UL DIMENSION UL(K7,i),X< ID, II, Bill COMMON /CNT/ ISTEP,NUMDE,NUMSO NN = N ♦ 1 NUMSO*NUMSO*l FIRST THE RESULT VECTOR IS INITILIZED TO THE CONSTANT VECTOR D01I=*1,N X(K,I »=B( I) K2=K7/2 K3=K2+1 K4=K2-1 IF(K2.EQ. 1> G0T09 THIS SECTION PERFORMS THE OPERATIONS OF THE FIRST PART OF THE GAUSSIAN ELIMINATION ON THE CONSTANT VECTOR D02 I = 2,K2 Il» 1-1 J3 = K3 - I DO 2 J=1,I1 J2 = J3 ♦ J 2 X(K,I >=X(K,I)-X(K,J)*UL(J2 ,11 9 DO 3 I=K3,N D03 J =1,K2 J2=I-J J3=K3-J ) X(K,I >=X(K, I >-X=X(K,NJ>-UL(KJ,NJ)*X(K,NJ1I 6 X(K,NJ) =X(K,NJI/UL(K3,NJ» 5 DO 10 IBK=K3,N I = NN - IBK DO 8 J = 1,K2 JK=K3+J IJ= I+J 8 X{K,n = X(K,I)-UL( JK, I»*X(K, IJ) 10 X(K,I >=X(K,I)/UL(K3,II RETURN END ! SLV1170Q SLV1180( SLV1190I SLV1200< SLV1210( SLV12200; SLV12300 SLV12400 SLV12500 SLV12600 SLV12700 SLV12800! SLV12900 SLV13000 SLV13100 SLV13200 SLV13300 SLV13400 SLV13500 SLV13600 SLV13700 SLV13800 SLV13900 SLV14000 SLV14100 SLV14200 SLV14300 SLV14400 SLV14500 SLV14600! SLV14700 SLV14800 SLV14900 SLV15000 SLV15100 SLV15200! SLV15300 SLV15400' SLV15500 SLV15600 SLV15700 SLV15800 SLV15900 SLV16000 SLV16100 SLV16200 SLV16300 SLV16400 SLV16500 SLV16600 SLV16700 SLV16800 SLV16900 SLV17000 SLV17100 103 C c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c 2. 3. SUBROUTINE OIFFUN( X,Y,T,D IVER.MO, INT ,N, IDER,ERR,OY,L. 12, L3 I THE PURPOSE OF THIS PROGRAM IS TO PROVIDE 1. THE TIME DERIVATIVE AS THE RESULT OF THE DIFFERENTIAL EQUATION DURING INTEGRATION THE TIME DERIVATIVE WHEN STARTING THE TIME DERIVATIVE FOR THE CALCULATION OF THE JACOB I AN IN THE FIRST CASE L * L2 « IN THE SECOND CASE L*1,L2»0 AND THE PROGRAM ESTIMATES THE TIMF POINTsI IVE ° F ™ E B0UN ° ARY P0INTS 8V EXTRAPOLATION FROM INTERIOR ARE T BYpIsSED. CASE LC ° fl2 ' 1 AN ° ™ E ERR0R ESTIMATION CALCULATIONS THE VARIABLES IN THIS SUBPROGRAM ARE X Y T DIVER MO INT N IDER ERR DY L L2 DV DVP THE REST HOLDS THE HOLDS THE TIME CONTAINS INTERPOLA ORDER OF PROVIDE P STARTING NUMBER OF CONTAINS CONTAINS CONTAINS INDICATES INDICATES TEMPORARI TEMPORARI OF THE VARIABLE PARTITION POINTS FUNCTIONAL VALUES THE COEFFICIENTS OF ERROR ESTIMATE TING POLYNOMIALS THE METHOD OINTS OF ANALOG f INTC 1*1 II AND THE POINT OF ERROR ESTIMATE ANALOG (INT(2tIII EQUATIONS NAMES OF DERIVATIVES FOR ALL ORDERS ERROR ESTIMATE VECTOR THE TIME DERIVATIVE IF FINDING TIME DERIVATIVE AT BOUNDARY JACOBIAN CALCULATION LY STORES SPACE DERIVATIVES LY STORES ERROR ESTIMATE FOR DV S ACT AS COUNTERS OR POINTERS. IMPLICIT REAL * 8 (A-H,Q-Z» DIMENSION SV(6» 1 o^^«:!Um:!:o?;::^j;^:i: ,dwer,9,mo ' i, ' ,nt,2 ' i, ' ider « 4 '- DIMENSION XT(5),XS(5I,YT<3,5I,YS<3,5I COMMON /DVR/ DVY( 4 » ,COEF ( 4, 3 f 31 , ALPHA! 2,4, 3 >, DI V( 16, 5 ) THIS SECTION DETERMINES THE TIME DERIVATIVES AT THE BOUNDARIES IF(L3.EQ.0I GO TO 4 DO 3 IA*l,MO SV(IAJ*Y(1,IA,1I SV(IA*3I=Y(1,IA,NI CONTINUE CALL BDRY(X f Y,T,MO,N,DY) DO 5 IA = l,MO DFNIOOOO DFN10100 DFN10200 DFN10300 DFN10400 DFN10500 DFN10600 DFN10700 DFN1Q800 DFN10900 DFNllOOO DFNlllOO DFN11200 DFN1130O DFN11400 DFN11500 DFN11600 DFN11700 DFN11800 DFN11900 DFN12000 DFN12100 DFN1220O DFN12300 DFN12400 DFN12500 DFN12600 DFN12700 DFN12800 DFN1290C DFN13000 DFN1310O DFN13200 DFN13300 DFN13400 DFN13500 DFN13600 DFN13700 DFN13800 DFN13900 DFN14000 DFN14100 DFN14200 DFN143C0 DFN14400 DFN14500 DFN14600 DFN14700 DFN14800 DFN14900 ioU c c c 10 20 C C c c 25 30 40 50 60 ERR(IA,1 I » DY( IA,1 I ERR(IA,N> » DY( IA,NI CONTINUE THIS SECTION DETERMINES WHICH POINTS ARE USED IN THE DERIVATIVE NM1*N-1 DO 100 IA * 2.NM1 IG»IA IF( INTd, IAI.LT.O) GO TO 10 IB=IA+INT(1,IA» IC=IA-1 DETERMINES STARTING POINT FOR ERR AND THE POWERS OF SPACING GO TO 20 IB=IA+1 IC=IA*INT(1,IAI ID=INT(2,IAI XI=X(IAI-X{ IC) IF (XI . EQ. XITMP .AND, MO . EQ. MOLD! GO TO 25 MOLD = MO XITMP * XI XIS=XI*XI IF (MO .LT. 2) GO TO 25 XIC=XI*XIS IF (MO .LT. 31 GO TO 25 XIF=XI*XIC GO TO (30,70,80 I, MO SETS PARAMETERS AND FORMS DERIVATIVES FOR LAGRANGIAN IF(IDER(1>.EQ.1I DV(1I=(Y(1,1,IBI-Y(1,1,IC)I/(2.*XII IFC I0ERC2l.EQ.il DV( 2 1 = ( Y( 1 , 1 ,1 B 1-2. *Y( 1, 1 , IA >+Y( 1 , 1, IC ) I /XIS DVP(1I=0. DVP(2I=0. USES COEFFICIENTS OF INTERPOLATING POLYNOMIAL AND CALCULATES ERROR DERIVATIVES. IF (L2 .EQ. 1) GO TO 65 DO 50 IE=1,M0 DO 5D IF=1,4 IH=ID+IF-1 IFl IDER(M0|.EQ.1J DVP( MOI =D I VER ( I F, I E , I A I *Y( 1 , IE, IH )+DVP( M3 t IF(IDER(MO+ll.EQ.il DVP ( M0*1I =D IVER ( I F+4, IE, I Al *Y( 1 , I E, IH l+DVP( MO 1+1 I CONTINUE IF(MO.NE.l) GO TO 60 IFC IDERC2l.EQ.il DVP(2I=D IVER (9, 1 , I A)*Y( 1, 1, ID+4 KDVPl 21 USES DIFFERENTIAL EQUATION TO OBTAIN TIME DERIVATIVES AND ERRORS IF( IDER(M0I.EQ.1I DVP( MOI =DVP(MO I *DV( MO I IF( IDER(M0*1I.EQ. II DVP ( MO* 1 I =DVP (MOU I +DVC MO+1 I DFN15000 OFN15100 DFN15200 DFN15330 DFN15400 DFN15500 DFN15600 DFN15700 DFN15800 DFN15900 DFN16000 DFN16100 DFN16200 DFN163Q0 DFN16400 DFN16500 DFN16600 DFN16700 DFN16800 DFN16900 DFN17000 0FN17100 DFN17200 DFN17300 DFN17400 DFN17500 DFN17600 DFN17700 DFN17800 DFN17900 DFN18000 DFN18100 DFN18200 DFN18300 DFN18400 DFN18500 DFN18600 DFN18700 DFN18800 DFN18900 DFN19000 DFN19100 DFN19200 DFN19300 DFN19400 DFN19500 DFN19600 DFN19700 DFN198D0 DFN19900 DFN20000 DFN20100 DFN20200 DFN20300 DFN20400 105 65 100 103 105 110 120 CALL UFUN(M0,XUA),Y CONTINUE ) RETURN IF XI XN DO II 12 C L .EQ. • xui ' X(N) 110 I * 1,5 ' I ♦ 1 » N - I XT(II » X(Il» xsm * xii2» 00 110 J * 1,M0 YT(J,II » DYU.Il) YS( J,I) « DY( J, 12) CONTINUE CALL DERIVUT.YT, I0ER CALL DERIV(XS,YS, IDER DO 120 I ■ 1,M0 DY( I» II * DC II DY( I,NI * ER(I) CONTINUE RETURN ,M0, 0, 1,D,M0,X1,D, XII »M0,0, l,0,MO t XN,ER,XII SETS PARAMETERS AND FORMS DERIVATIVES FOR HERMITE 70 DV(1I=Y(1,2,IGI DVP(1I*DV(1I IF( IDER(2).EQ.CI GO TO 75 DVP(2I=0. ic^i/(2!*xn (1,1 ' IB, " 2 '* Y{1,1,IG,4Y(1,1,ICn/,(IS ' ,Y,1,2,IB, " Y(U2 ' 75 IF( IDER(3).EQ.0I GO TO UO DVP(3»>=0. DV(3l=7.5000C0C*(Y(l t l,IBI-Y(l f l,IC)l/XIC-1.500000D*(Y(l,2.IBI 1 ♦8.*Y(1,2,IG) ♦ Y(1,2,ICII/XIS GO TO 40 SETS PARAMETERS AND FORMS DERIVATIVES FOR TRIMITE 80 DVIll=Y(l,2 t IG» DVP(1)*DV(1J DV(2I=Y(1,3,IG» DVP(2I=DV(2I IF( IDER(3). EQ.OI GO TO 85 DVP(3I»0. 1 +Y(l t 2, IC»l*18.000000»Y(l f 2,IGII/XIS+.37500000*(Y(l,3,IB» - OFN205OO DFN20600 DFN20700 DFN20800 DFN20900 DFN21000 DFN21100 DFN21200 DFN21300 DFN21400 DFN21500 0FN21600 DFN21700 DFN21800 DFN21900 DFN22000 DFN22100 DFN22200 OFN22300 DFN22400 DFN22500 DFN22 600 DFN22700 DFN22800 DFN22900 0FN23000 DFN23100 DFN23200 DFN23300 DFN23400 DFN23500 DFN23600 DFN23700 DFN23800 DFN23900 DFN24000 DFN24100 IDFN24200 0FN24300 DFN24400 DFN24500 DFN24600 DFN24700 DFN24800 DFN24900 DFN25000 0FN25100 DFN25200 DFN25300 DFN25400 DFN25500 DFN25600 DFN25700 OFN25800 DFN25900 106 c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c 2 Y(lf3tICII/XI DFN26O00 85 IF{ IDERm.EQ.O) GO TO 40 DFN26100 DVP(4)«0. DFN26200 DV(4)«72.000000*(Y(1,1,IBI-2.*Y(1,1,IG)«-Y(1,1,IC) )/XI F-19. 500000* DFN26300 l(Y{l,2t IBJ-Y(l,2,IC)>/XIC*1.5COOOOa*(Y(l,3,IB)-24.00000C*Y(l,3,IG)DFN26400 2 ♦Y(1,3,ICM/XIS DFN26500 GO TO 40 DFN26600 END DFN26700 SUBROUTINE 01 FSUB (N ,T ,Y .SAVE, H, HMIN, HMAX, EPSt YMAX, ERROR ,KFLAG , USTART,MAXDER,PSAVE,K7,HMY,X,DIVER,INT,IDER,ERR , MO , I OOUB, I QQQC, 2HSAV,PEQ,ER0RA,ITPSC) IMPLICIT REAL * 8(A-H,P-Z» REAL * 4 PEPSH,PRl,PR2tPR3,PERTST THE PARAMETERS TO THE SUBROUTINE DIFSUB HAVE THE FOLLOWING MEANINGS*. N SAVE HMIN HMAX THE NUMBER OF FIRST ORDER DIFFERENTIAL EQUATIONS. N MAY BE DECREASED ON LATER CALLS IF THE NUMBER OF ACTIVE EQUATIONS REDUCES. BUT IT MUST NOT BE INCREASED WITHOUT CALLING WITH JSTART » 0. THE INDEPENDENT VARIABLE. AN 8 BY N ARRAY CONTAINING THE DEPENDENT VARIABLES AND THEIR SCALED DERIVATIVES. Y(J+1,II CONTAINS THE J-TH DERIVATIVE OF Y(I) SCALED BY H**J/FACTORIAL(J) WHERE H IS THE CURRENT STEP SIZE. ONLY Yd, II NEED BE PROVIDED BY THE CALLING PROGRAM ON THE FIRST ENTRY. IF IT IS DESIRED TO INTERPOLATE TO NON MESH POINTS THESE VALUES CAN BE USED. IF THE CURRENT STEP SIZE IS H AND THE VALUE AT T + E IS NEEDED, FORM S = E/H, AND THEN COMPUTE NQ YUMT+EI » SUM YU + 1,1 )*S**J J=0 A BLOCK OF AT LEAST 12*N FLOATING POINT LOCATIONS USED BY THE SUBROUTINES. THE STEP SIZE TO BE ATTEMPTED ON THE NEXT STEP. H MAY BE ADJUSTED UP OR DOWN BY THE PROGRAM IN ORDER TO ACHEIVE AN ECONOMICAL INTEGRATION. HOWEVER, IF THE H PROVIDED BY THE USER DOES NOT CAUSE A LARGER ERROR THAN REQUESTED, IT WILL BE USED. TO SAVE COMPUTER TIME, THE USER IS ADVISED TO USE A FAIRLY SMALL STEP FOR THE FIRST CALL. IT WILL BE AUTOMATICALLY INCREASED LATER. THE MINIMUM STEP SIZE THAT WILL BE USED FOR THE INTEGRATION. NOTE THAT ON STARTING THIS MUST BE MUCH SMALLER THAN THE AVERAGE H EXPECTED SINCE A FIRST ORDER METHOD IS USED INITIALLY. THE MAXIMUM SIZE TO WHICH THE STEP WILL BE INCREASED DSB10000 DSB10100 DSB10200 DSB10300 DSB10400 DSB10500 0SB10600 DSB 10700 DSB10800 DSB10900 DSB11000 DSBIHOO DSB11200 DSB11300 DSB11400 DSB11500 DSB11600 DSB11700 DSB11800 OSB11900 DSB12000 DSB12100 0SBI2200 DSB12300 DSB12400 DSB12500 DSB12600 DSB12700 DSB12830 DSB12900 DSB13000 DSB13100 DSB13200 DSB13300 DSB13400 0SB13500 DSB13600 DSB13700 DSB13800 DSB13900 DSBKOOO DSB14100 107 EPS MF YMAX ERROR KFLAG JSTART MAXDER AN ARRAY OF OF EACH Y THE ERROR TEST CONSTANT, SINGLE STEP ERROR FSTimatpc DIVIDED BY YMAXUI MUST BE LESS THAN *?*, ESTlMATES IwUSTEO^^^ErSllJ* STEP AND/ ° R ° RDER IS THE METHOD INDICATOR. THE FOLLOWING ARE ALLOWED.. AN ADAMS PREDICTOR CORRECTOR IS USED. 1 A MULTI-STEP METHOD SUITABLE FOR STIFF l Y nt T lrLl S r° SBD ' IT WILL ALS0 WORK FOR NON STIFF SYSTEMS. HOWEVER THE USER MUST PROVIDE A SUBROUTINE PEDERV WHICH EVALUATES THE PARTIAL DERIVATIVES OF THE DIFFERENTIAL EQUATIONS WITH RESPECT TO THE Y«S. THIS IS DONE BY CALL PEDERV{T,Y,PW,M>. PW IS AN N BY N ARRAY WHICH MUST BE SET TO THE PARTIAL OF THE I-TH EQUATION WITH RESPECT TO THE J DEPENDENT VARIABLE IN PW(I.J). PW IS ACTUALLY STORED IN AN M BY M ARRAY WHERE M IS THE VALUE OF N USED ON THE FIRST CALL TO THIS PROGRAM. 2 THE SAME AS CASE 1, EXCEPT THAT THIS SUBROUTINE COMPUTES THE PARTIAL DERIVATIVES BY NUMERICAL DIFFERENCING OF THE DERIVATIVES. HENCE PEDERV IS N LOCATIONS WHICH CONTAINS THE MAXIMUM SEEN SO FAR. IT SHOULD NORMALLY BE SET TO AN ON 4 " R ^ E ? F E ^0 E 1 N M ^Ic S H H C H ^^ C ^ T . MNS THE """"» i COMPLETION CODE WITH THE FOLLOWING MEANINGS.. THE STEP WAS SUCCESSFUL THE STEP WAS TAKEN WITH H = HMIN, BUT THE REQUESTED ERROR WAS NOT ACHIEVED. THE MAXIMUM ORDER SPECIFIED WAS FOUND TO BE TOO LARGE. CORRECTOR CONVERGENCE COULD NOT BE ACHIEVED FOR H .GT. HMIN. THE REQUESTED ERROR IS SMALLER THAN CAN BE HANDLED FOR THIS PROBLEM. AN INPUT INDICATOR WITH THE FOLLOWING MEANINGS.. -1 REPEAT THE LAST STEP WITH A NEW H C PERFORM THE FIRST STEP. THE FIRST STEP MUST BE DONE WITH THIS VALUE OF JSTART SO THAT THE SUBROUTINE CAN INITIALIZE ITSELF. .c-r.oT +1 TAKE A NEW STEP CONTINUING FROM THE LAST. ir T IS ^\™ NQ ' THE CURR ENT ORDER OF THE METHOD AT EXIT. NQ IS ALSO THE ORDER OF THE MAXIMUM DERIVATIVE AVAILABLE. "AXIMUM THE MAXIMUM DERIVATIVE THAT SHOULD BE USED IN THE METHOD. SINCE THE ORDER IS EQUAL TO THE HIGHEST SI ltl l T E " SE0 ' THIS RESTR ICTS THE ORDER. IT mIsT BE LESS THAN 8 FOR ADAMS AND 7 FOR STIFF METHODS. ♦ 1 -1 -2 -3 -4 DSB14200 DSB143Q0 DSB14400 DSB14500 DSB14600 DSB14700 DSB14800 DSB14900 DSB15000 DSB15100 DSB15200 0SB15300 DSB15400 DSB15500 DSB15600 DSB15700 DSB15800 DSB15900 DSB16000 0SB16100 DSB16200 DSB16300 0SB16400 DSB16500 DSB16600 DSB16700 0SB16800 DSB16900 DSB17000 DSB17100 DSB17200 DSB17300 DSB17400 DSB17500 DSB17600 0SB17700 DSB17800 DSB17900 DSB18000 DSB18100 DSB18200 DSB18300 0SB18400 DSB18500 DSB18600 DSB18700 DSBI8800 DSB18900 DSB19000 DSB19100 DSB19200 DSB19300 DSB19400 DSB19500 DSB19600 PSAVE HMY HSAV PEG ERR ERORA COMMON DIMENS DIMENS DIMENS DIMENS I A(8I, DIMENS DIMENS DATA X A BLOCK OF AT LEAST N**2 FLOATING POINT LOCATIONS, STORES THE ASYMPTOTIC ERROR VECTOR SAVES THE PREVIOUS VALUES OF THE ASYMPTOTIC ERROR AND IS USED TO STORE DERIVATIVES STORES THE JACOBIAN CONTAINS THE SPATIAL ERROR ESTIMATE CONTAINS THE ESTIMATED, INTRODUCED ERROR /DVR/ DVY(4) ,COEF(4,3,3),ALPHA(2,4,3l,DIV(16,5l ION XT(5t2l v XI(2l ION YT(3,5> ION V <1 I ,D ION Y(5,15 t PERTST(7,3J ION ERORAd IONHMY(2,l» ,HSAV( 3,1) ,PEQIK7,1» ITMP /O.Q/ IVERCDt INTClIt IDERdlfERRCll , XIII YMAX<1U SAVEI7 , 1 1 , ERROR ( 1 1 , PSAVE( K7 , U , I THE COEFFICIENTS ORDER, THEREFORE DATA PERTST IN ONL /12.G 12.0 1 1. , 6. , 2. , • 3 DATA MF,MTYP /2,1 DATA A(2I / -1.0/ IF( ITPSC .EO. 1 J ! I TPSC » N14=N/M0 JRET =2 K2=K7/2 K3=K7/2+l IF = SAVE(N11,1»*H ERROR( II = EPSN 144 CONTINUE CALL DIFFUN(X,Y,T,DIVER,M0,lNT,N14,IDER,ERR,SAVE(N2,n,0,l,ll GO TO 310 142 CONTINUE CALL DIVERR(PE0,N,ERR0R,ERR,HSAV(NY2,1»,HMY,PERTST(NQ,1I,K7,1,M0 1 , X I J RET* 2 DO 150 I =* 1,N ERROR (I) » 0. NY11=NY1+I HMY(2,I)»HSAV(NY11,1»*H 150 CONTINUE HNEW » H K * 2 GO TO 100 C C REPEAT LAST STEP BY RESTORING SAVED INFORMATION. C 160 IF (NQ.EO.NQOLDI JSTART « 1 T»TOLD NO * NQOLD DS825200 DSB25300 DSB25400 DSB25500 DSB25600 DSB25700 DSB25800 DSB25900 DSB26000 DSB26100 DSB26200 DSB26300 DSB26400 DSB26500 DSB26600 DSB26700 DSB26800 DS826900 DSB27000 DSB27100 DSB27200 DSB27300 DS827400 DSB27500 DSB27600 DSB27700 DSB27800 DSB27900 DSB28000 DSB28100 DSB28200 DSB28300 DSB28400 DSB28500 DSB28600 DSB28700 DSB28800 0SB28900 DSB29000 DSB29100 DSB2920C DSB29300 DSB29400 DSB29500 DSB29600 DSB29700 DS829800 DSB29900 DSB30000 DSB30100 DSB30 200 DSB30300 DSB30400 DSB30500 DSB30600 110 K GO ■ NQ ♦ 1 TO 120 SET THE COEFFICIENTS THAT DETERMINE THE ORDER AND THE METHOD TYPE. CHECK FOR EXCESSIVE ORDER. THE LAST TWO STATEMENTS OF THIS SECTION SET IWEVAL .GT.J IF PW IS TO BE RE-EVALUATED BECAUSE OF THE ORDER CHANGE, AND THEN REPEAT THE INTEGRATION STEP IF IT HAS NOT YET BEEN DONE (IRET * il OR SKIP TO A FINAL SCALING BEFORE EXIT IF IT HAS BEEN COMPLETED (IRET ■ 21. 170 IF(NQ .GT. 41 GO TO 190 GO TO (221, 222, 223, 224), NQ 190 KFLAG » -2 RETURN C C C C c c c c c THE FOLLOWING COEFFICIENTS SHOULD BE DEFINED TO THE MAXIMUM ACCURACY PERMITTED BY THE MACHINE. THEY ARE IN THE ORDER USED., -1/2 -1/2,-1/2 -1/2,-5/8,-1/8 -1/2,-4 7/72,-1/6,-1/72 221 A(l )*-. 5000000000 GO TO 230 222 A(l »=-. 5000000000 A(3) = A(1» GO TO 230 223 AC1)=-. 5000000000 A(3»=-.6250000C0C0 A(4)=-. 125000000000 GO TO 230 224 A(l »=-. 500C0000000 A(3 )=-. 6527777777777778 A(4)=-. 1666666666666667 A( 5 »=-. 01 38888888 8888889 230 K = NQ>1 IDOUB = K 231 CONTINUE ENQ2 = .5/FL0AT(NQ ♦ II ENQ3 = .5/FL0AT(NQ ♦ 2» ENQ1 = 0.5/FL0AT(NQI PEPSH = EPS EUP = (PERTST(NQ,2J*PEPSH)**2 E = (PERTST(N0,1»*PEPSHI**2 EDWN=(PERTST(NQ,3I*PEPSH»**2 IF (EDWN.EQ.OI GO TO 780 BND = EPS*ENQ3/DFL0AT(N> IF( ITPSC .EQ. 1 I GO TO 2 240 IWEVAL = MF GO TO ( 2 50 , 680 >,IRET THIS SECTION COMPUTES THE PREDICTED VALUES BY EFFECTIVELY MULTIPLYING THE SAVED INFORMATION BY THE PASCAL TRIANGLE MATRIX. OSB30 700 DSB30800 DSB30900 DSB31000 DSB31100 DSB31200 DSB31300 DSB31400 DSB31500 DSB31600 DSB31700 0SB31800 DSB31900 DSB32000 DSB32100 DSB32200 0SB32300 DSB32400 DS832500 DSB32600 DSB32700 DSB32800 DSB32900 DSB33000 DS833100 DSB33200 DSB33300 DSB33400 DSB33500 DSB33600 DSB33700 DSB33800 DSB33900 DSB34000 DSB34100 DSB34200 DSB34300 DSB34400 DSB34500 DSB34600 DSB347C0 DSB34800 DSB34900 DSB35000 DSB35100 DSB35200 DSB35300 DSB35400 DSB35500 DSB35600 DSB35700 DSB35800 DSB35900 DSB36000 DSB36100 c 111 c 250 T = T ♦ H DSB36200 DO 260 J » 2,K DSB36300 DO 260 Jl = J,K DSB36400 J2 » K - Jl ♦ j - 1 0S836500 00 260 I = 1,N 0SB36600 ^ ,f,j •%n*i\^ i £» t vA\»?: l \ un * HHv,l • ,, ♦ H " Y,^ • I, lilHsH C UP TO 3 CORRECTOR ITERATIONS ARE TAKEN. CONVERGENCE IS TESTED n1«l7?n? C tV^Or'tIst^OnItanT." L6 " ™ N "» """ « ™™™°» "»™ S ™5!!f,?. w ?! I^SM'Si ^^SkV^^,* i!?° ? rl T H !^iFACTORIAL(K-ll»A(KI», AND IS THEREFORE PROPORTIONAL DSB37600 C TO THE ACTUAL ERRORS TO THE LOWEST POWER OF H PRESENT. IH**K» DSB377^S DO 270 I = 1,N DSB3780O 270 ERRORUI * 0. G 0SB37900 Li = 3 DSB38000 DO 430 L = 1.L1 , DSB38100 CALL DIFFUNU,Y,T,DIVER,M0,INT,N14,IDER, ERR, SAVECN2,1», 0,0,01 ollllllo llrl^J, HAS BEEN A CHANGE 0F 0RDER 0R T HERE HAS BEEN TROUBLE DSB38^0 C WITH CONVERGENCE, PW IS RE-EVALUATED PRIOR TO STARTING THE DSR3fl^n C THF^ E , C p? R rn TE ? A II 0N IN ™ E CASE ° F STIFF METH3 °S. IWEvIl IS ollll^O C THEN SET TO -1 AS AN INDICATOR THAT IT HAS SEEN DONE. DSB38800 IFdWEVAL.LT.il GO TO 370 nfo!oo2° GO TO 310 DSB39000 C DSB39100 C SETS THE JACOBIAN FOR USE WITH ASYMPTOTIC ERROR EQUATION DSBsJsSo 271 00 272 1=1, K7 DSB39400 DO 272 J=1,N DSB39500 272 PEQ(I,J)=PSAVE(I,J) DSB3960O C DSB39700 C PREPARES TO DECOMPOSE THE MATRIX PSAVE DSbI'wo IFURET.EQ.il GO TO 142 n?^°? 00 IFJJRET.EQ.3I GO TO 721 Sco^JSS 273 R = A(1»*H DSB40200 DO 28C I =1,«7 DSB4D300 K8=K2+2-I DSB4G400 K8= MAX0{1,K8I DSB40500 K9=N*K2 + 1-I DS6A0600 K9=MINl'(N,K9| DSB40700 DO 28C J=K8,K9 DSR4C8C0 IF( J .LE. MO .OR. J .GT. NMO. GO TO 275 DSB^O PSAVEd.JI = PSAVE(I,JI>R n«2J?«n GO TO 280 DSB4U00 275 PSAVEd.JI = PSAVE(I,JI * A(ll S?2tJ!?2 280 CONTINUE DSB41330 29C DO 300 1=1, N DSB41400 IF( I .LE. MO .OR. I .GT. NMO I GO TO 300 SIbJuOO 112- PSAVE(K3,II » 1.0 ♦ PSAVE(K3,I) 300 CONTINUE IWEVAL » -1 DO 308 I » It MO J » K2-I+2 J3 = DO 303 J2 « JtKT IF(DABS(PSAVE(J2, I) ) .LT. 1.0D-20) J3 » J3 ♦ 1 303 CONTINUE J = K7-J+1 IFU3 .GE. J) PSAVE(K3,II = l.ODO J3 = C Jl = NMO ♦ I J = K3+M0-I DO 305 J2 = l t J IF(DABS(PSAVE*SAVE(6,I I SAVE(6,II 42 NT NT - 1 42 5 42 7 430 00 420 I Y(1,I) Y(2,I> ERRORUI » ERRORm ♦ SAVE(6,I) IF (DABS(SAVE<6,II >.LE. ( BND*YMAX < 1 1 I I CONTINUE CALL DIVERR(PEQ,N,ERROR,ERR,HSAV(NY2,1I,HMY,PERTST(NQ,1I,K7,0,MO lfXI DO 4251 = 1, N NY11=NY5+I NY12=NY1+I HSAV (NYll,l»=HMY(2,n-HSAV(NY12,ll*H IF( I ,LE. MO .OR. I .GT. NMO) HSAV(NY11,1| * - HSAV(NY12,1» CONTINUE SOLVE THE ASYMPTOTIC ERROR EQUATION CALL SOLVE(N,PSAVE,HSAV(NY5*l,ll,HSAV,K7,3,3) DO 427 I»1,N HMYt 1 , I » = HMY ( 1 , II +A ( 1 l*HSAV( 3, I I HMY(2,I)=HMY(2,II- HSAV<3,II IF (NT.LE.OI GO TO 490 CONTINUE THE CORRECTOR ITERATION FAILEO TO CONVERGE IN 3 POSSIBILITIES ARE CHECKED FOR, IF H IS ALREADY THIS IS EITHER ADAMS METHOD OR THE STIFF METHOD MATRIX PW HAS ALREADY BEEN RE-EVALUATED, A NO CONVERGENCE EXIT IS TAKEN. OTHERWISE THE MATRIX PW IS RE-EVALUATED AND/OR THE STEP IS REDUCED TO TRY AND GET CONVERGENCE, TRIES, VARIOUS HMIN AND IN WHICH THE 440 T=TOLD IF UH.LE. (HMIN*1. 00001 ) ) . AND. ( ( I WEVAL - MTYP I. LT.-l I I IF HMF.EQ.OI.OR. ( I WEVAL. NE ,0 I » RACUM = RACUM*0.25DO IWEVAL = MF IRET1 « 2 GO TO 750 KFLAG * -3 DO 480 I » HMY(lfl) * HMYC2.II = DO 480 J Y(J,II H = HOLD NQ = NQOLD J ST ART « NQ RETURN GO TO 460 460 470 480 i,N HSAV(1,II HSAV(2,I I = i,K = SAVEIJtll THE CORRECTOR CONVERGED ANO CONTROL IS IF THE ERROR TEST IS O.K., AND TO 540 PASSED TO STATEMENT OTHERWISE, 520 DSB47200 DSB47300 DSB47400 DSB47500 DSB47600 DSB47700 DSB47800 DSB47900 DSB48000 DSB48100 DSB48200 DSB48300 DSB48400 DSB48500 DSB48600 DSB48700 DSB48800 DSB48900 DSB49000 DSB49100 DSB49200 DSB49300 DSB49400 DSB49500 DSB49600 DSB49700 DSB49800 DSB49900 DSB50000 DSB50130 DSB502OO DSB50300 DSB50400 DSB50500 . DSB50600 DSB50700 DSB508D0 DSB50900 DSB51000 DSB51100 DSB51200 DSB51300 DSB51400 DSB515C0 DSB51600 DSB51700 DSB51800 DS851900 DSB52000 DSB52100 DSB52200 DSB52300 DSB52400 DSB52500 DSB52600 Ill* C IF THE STEP IS O.K. IT IS ACCEPTED. IF I DOUB HAS BEEN REDUCED C TO ONE, A TEST IS MADE TO SEE IF THE STEP CAN BE INCREASED C AT THE CURRENT ORDER OR BY GOING TO ONE HIGHER OR ONE LOWER. C SUCH A CHANGE IS ONLY MADE IF THE STEP CAN BE INCREASED BY AT C LEAST 1.1. IF NO CHANGE IS POSSIBLE I DOUB IS SET TO 10 TO C PREVENT FURTHER TESTING FOR 10 STEPS. C IF A CHANGE IS POSSIBLE, IT IS MADE AND IDOUB IS SET TO C NQ ♦ 1 TO PREVENT FURTHER TESING FOR THAT NUMBER OF STEPS. C IF THE ERROR WAS TOO LARGE, THE OPTIMUM STEP SIZE FOR THIS OR C LOWER ORDER IS COMPUTED, AND THE STEP RETRIED. IF IT SHOULD C FAIL TWICE MORE IT IS AN INDICATION THAT THE DERIVATIVES THAT C HAVE ACCUMULATED IN THE Y ARRAY HAVE ERRORS OF THE WRONG ORDER C SO THE FIRST DERIVATIVES ARE RECOMPUTED AND THE ORDER IS SET C TO 1. C 490 D = 0.0 DO 5C0 I » l t N 500 D = D ♦ (ERR0R(I)/YMAX( 1)1**2 IWEVAL * IF (D.GT.EI GO TO 540 IF (K.LT.3) GO TO 520 C C STEP SUCCEEDED, UPDATE REMAINING VECTORS 510 520 3,K * 1,N « Y(J,I) ♦ A(J)*ERROR(I) DO 510 J - DO 510 I Y( J, I I KFLAG * +1 HNEW = H IF ( IDOUB. LE. II GO TO 550 IDOUB = IDOUB - 1 IF (IDOUB.GT. 1) GO TO 700 DO 53C I = 1,N 530 SAVEI7 ,1) =» ERROR(I) GO TO 700 ERROR TOO LARGE, REDUCE FAILURE FLAG AND TRY AGAIN 540 KFLAG = KFLAG - 2 IF (H.LE. (HMINM.OOOCl) ) GO TO 740 T = TOLD IF (KFLAG. LE. -11) GO TO 720 CHECKS TO SEE IF CHANGES SHOULD BE MADE IN STEP SIZE OR ORDER 550 PR2 = (D/E)**ENQ2*1.2 PR3 = l.E+20 IF UNQ.GE.MAXDER). OR. (KFLAG. LE.-l) ) GO TO 570 D = O.C DO 560 I = 1,N 560 D = D ♦ ((ERROR(I) - SAVE(7 , I ) )/YMAX( I ) 1**2 PR3 = (D/EUP)**ENQ3*1.4 570 PR1 = l.E+20 IF (NQ.LE.l) GO TO 590 D = 0.0 DSB52700 DSB52800 DSB52900 DSB53000 DSB53100 DSB53200 DSB53300 DSB53400 DSB53500 DSB53600 DSB53700 DSB53800 DSB53900 DSB54000 DSB54100 DSB54200 DSB54300 DSB54430 DSB 54500 DSB54600 DSB 54700 DSB54800 DSB54900 DSB55000 DSB55100 DSB55200 DSB55300 DSB55400 DSB55500 DSB55600 DSB55700 DSB55800 DSB55900 DSB56000 DSB56100 DSB56200 DSB56300 DSB56400 DSB56500 DSB56600 DSB56700 DSB56830 DSB56900 DSB57000 DSB57100 DSB57230 DSB57300 DSB57400 DSB57500 DSB57600 DSB57700 DSB57800 DSB57900 DSB58000 DSB58100 115 00 580 I * i,N 580 - ♦ (Y(K,I l/YMAX(I>l**2 PR1 * (D/EDWN>»*ENQ1*1.3 590 CONTINUE IF (PR2.LE.PR3) GO TO 650 IF (PR3.1T.PRU GO TO 660 600 R » 1.0/AMAXKPR1 ,l.E-4> NEWQ « NQ - 1 610 ID0U8 » 10 IF ((KFLAG.EQ.ll. AND.(R.LT. (l.lil ) GO TO 700 IF (NEWQ.LE.NQ) GO TO 630 COMPUTE AN ADDITIONAL SCALED DERIVATIVE IF ORDER IS INCREASED DO 620 I = 1,N 620 Y(NEWQ*1,II « ERRORd »*A(K) /DFLOAT( Kl 630 K = NEWQ ♦ 1 IF (KFLAG.EQ.ll GO TO 670 RACUM = RACUM*R IRET1 » 3 GO TO 750 640 IF (NEWQ. EQ.NQ) GO TO 250 NQ * NEWQ GO TO 170 650 IF (PR2.GT.PR1) GO TO 600 NEWQ * NQ R = 1.0/AMAXKPR2,l.E-4) GO TO 610 660 R =* 1.0/AMAX1(PR3,1.E-4I NEWQ » NQ ♦ 1 GO TO 610 670 IRET » 2 R = DMIN1(R,HMAX/DABS(H)| H » H^R HNEW = H IF (NQ. EQ.NEWQ) GO TO 680 NQ = NEWQ GO TO 170 680 Ri = 1.0 DO 690 J = 2,K Rl = R1*R DO 690 I = 1,N IF(J .EQ. 2) HMY(J,IJ » HMY(J,I» *R1 690 Y(J,I» = Y(J,n»Rl IDOUB = K 700 DO 710 I * 1,N YMAX(I) = DMAX1(YMAX(II,DABS(Y(1,I))J 710 ER0RA(n=(HMY(l, I l-HSAVI 1 ,11) /YMAXdl JSTART » NQ RETURN 720 IF (NQ.EQ.1J GO TO 780 CALL DIFFUN(X,Y,T, DIVER, MO, I NT, N14, 1 DER, ERR , SAVE( N2, 1 » , 1, , 1 J K c 1*3 EPSN = .16700 ♦ EPS/OFLOAT(NI DO 722 1=1, N 0S858200 0SB58300 DSB58400 DSB58500 0SB58600 DSB58700 DSB58800 0SB58900 DSB59000 0SB59100 DSB59200 DSB59300 DSB59400 DSB59500 DSB59600 DSB59700 DSB59800 DSB59900 DSB60000 DSB60100 OSB60200 DSB60300 DSB63400 DSB60500 0SB60600 DSB60700 DSB60800 DSB6C900 DSB61000 DSB61100 0SB61200 DSB61300 DSB61400 DSB61500 DSB61600 0SB61700 DSB61800 DSB61900 DSB62000 DSB62100 DSB62200 DSB62300 DSB62400 DSB62500 DSB62600 DSB62700 DSB628G0 DSB62900 DSB63000 DSB63100 DSB63200 DSB63300 DSB63400 DSB63500 DSB63600 .116 722 721 730 740 C C TH C TO C 750 760 770 780 R .« H/HOLO Yd, II « SAVEd.I I Nil * Nl ♦ I SAVE(2,I) ■ HOLD*SAVE(NU,l» Y(2,I I = SAVE(2,I I * R ERRORd I « EPSN CONTINUE CALL DIFFUN(X,Y,T,DIVER,MQ, INT,N14,IDER,ERR,SAVE(N2,ll,0,i,il GO TO 310 CONTINUE JRET-2 CALL DIVERR(PEQ,N,ERR0R,ERR,HSAV(NY2,1),HMY,PERTST(NQ,1I,K7,1,M0 1,XI R = H/HOLD DO 73C I * 1,N ERRORd I » 0,0 HMY(l,II=HSAVd,I I NY11=NY1*I HSAV(2,I)*H0LD*HSAV(NY11,1> HMY(2,I l=HSAV(2,I l*R CONTINUE NQ * 1 KFLAG » 1 GO TO 170 KFLAG * -1 HNEW * H JSTART * NQ RETURN IS SECTION SCALES ALL VARIABLES CONNECTED WITH H AND RETURNS THE ENTERING SECTION. RACUM = DMAX1 ( DABS ( HMIN/HOLD I , RACUM ) RACUM = DMIN1 (RACUM, DABS( HMAX/HOLDII Rl = 1.0 DO 760 J = 2,K Rl = R1*RACUM DO 760 I = 1,N IF( J .EG. 21 HMY(J,II = HSAV(J,II * Rl Y(J, I I = SAVE(J,I I*R1 H = HOLD*RACUM DO 770 I = 1,N HMYd,II=HSAVd,I I Yd, II = SAVEd ,1 I IDOUB GO TO KFLAG GO TO END = K ( 130 = -4 470 , 2 50 , 640 »,IRET1 DSB63700 DSB63800 DSB63900 DSB 640001 DSB64100 DSB64200 DSB6430C DSB64400 DSB64500 DSB 64600 DSB64700 DSB64800 DSB64900! DSB65000 DSB651001 DSB65200 DSB65300 DSB65400 DSB65500 DSB65600 DSB65700 DSB65800 DSB65900 DSB66000 DSB66100 DSB66200 DSB66300! DSB66400J DSB66500 i DSB66600 DSB66700 DSB66800 DSB66900 DSB67000 DSB67100 DSB67200': DSB67300" DSB67400' 1 DSB67500 DSB67600 DSB67700 DSB67800 DSB67900 DSB68000 DSB68100 DSB68200 DS868300 DSB68400 DSB 68500 SUBROUTINE DI VERR (PEQ,N, ERROR ,ERR,DY, Y, PERTST ,K7,M,M0 ,X) DVR10000 117 c c c c c c c c c c c c c c c c THE PURPOSE OF THIS SUBPROGRAM IS TO CALCULATE THE TIME DERIVATIVE OF THE ASSYMPTOTIC ERROR EQUATION THE VARIABLES IN THIS SUBPROGRAM ARE PEQ DY N ERROR ERR Y PERT ST K7 PARTIALS OF THE DIFFERENTIAL EQUATIONS WRT THE VARIABLES THE TIME DERIVATIVE VECTOR OF THE ASYMPTOTIC ERROR EQUATION NUMBER OF EQUATIONS THE TIME BASED TRUNCATION ERROR VECTOR THE SPACE BASED TRUNCATION ERROR VECTOR VALUE VECTOR OF THE ASYMPTOTIC ERROR CONSTANTS ASSOCIATED WITH THE ORDER THE NUMBER OF DIAGONALS IN THE JACOBIAN IMPLICIT REAL * 8 IA-H.O-ZI REAL * 4 PEPTST COMMON /DVR/ DVYC4) ,COEF ( 4, 3 ,3 ) , ALPHA ( 2,4, 3 > , DI V( 16, 5 1 DIMENSION DY(1» ,P EQ < K7, 1 ) ,E RROR(l) ,ERR<1> ,Y(2,1I DIMENSION XT(5I ,XS< 5 1 ,YT< 3, 5 1 ,VSI 3, 5 1 DIMENSION D(4),ER<4> DIMENSION X(l» K8 * K7/2 +1 K2 » K8 ♦ I K3 * K8 ♦ N DO 10 IA » 1,N Ni = MAXO(K2-IA,l ) N2 = MIN0U3-IA.K7I Kl = IA - K8 c c c DY(IA)»- (ERR(IA)I ♦ ( ERROR! I A J/PERTSTI MULTIPLIES THE APPROPRIATE VARIABLES BY THE PARTIA DO 5 IB = Nl , N2 DY(IA) * PEQ(IB.IA»*YU,K1*IB> ♦ DY(IA» 5 CONTINUE 10 c c c CONTINUE THIS SECTION DETERMINES THE VALUES AT THE BOUNDARY IF(M .EQ. 0) RETURN L=N/MO NO * i xi = xm XN » X(LI DO 115 J * 1,M0 DO 110 I » 1,5 11 » I ♦ 1 12 * L - I XT( I) = Xdll XS( II ■ X ( I 2 1 11 = I*MO*J YT(1,II * DY( H 1 12 * M0*(I2-1) ♦ J DVR10100 DVR1D200 DVR10300 DVR10400 DVR10500 DVR10600 DVR13700 DVR10800 DVR10900 DVR11000 DVR11100 DVR11200 DVR11300 DVR11400 DVR11500 DVRH600 DVR11700 DVR11800 DVR11900 DVR12000 DVR12100 DVR12200 DVR12300 DVR12400 DVR12500 DVR1260O DVR12700 DVR128C0 DVR12900 DVR13000 DVR13100 DVR13200 DVR1330O DVR13400 DVR1350O DVR13600 DVR13700 DVR13800 DVR13900 DVR14000 DVR14100 DVR14200 DVR14300 DVR144Q0 DVR14500 DVR14600 DVR147C0 DVR14800 DVR14900 DVR15000 DVR15100 DVR15200 DVR15300 DVR15400 DVR15500 118 YSCl.II - DYU2I . UB1C1AA 110 CONTINUE DVR15600 CALL DERIVCXT,YT,IDER,NO,0,l,C,Na,Xl,D,XII Dwlslnn CALL DERIV(XS,YS,IDER,NO,OU.O,NO.XN,ER,XII dvrisq™ dy(jj » Ddi n»„:;; 00 n * mo*(l-ii*j RXSJ!? 00 DYCIll » ERC1I 2X?} 6100 115 CONTINUE DVR16200 RETURN DVR16300 EN0 DVR16400 0VR16500 119 APPENDIX B EXAMPLE 1 120 903 30 50 IMPLICIT REAL *81A-H,Q-Z» DIMENSION Y(2,513),YTST<513»,M{8),E<513»,F{513I,TTST(8I DATA M/3, 7, 15, 31, 63, 127, 255, 511/ DATA TTST/.125,.25,.5,1.,2. ,4, ,8. ,16. / RP I =3. 141 592653 589793 QPI»RPI+RPl IP1 IP2 IP3 IP4 IP5 IP6 IP7 100 MAIN LOOP FOR DIFFERENT STEP SIZES DO 500 IA « 2,7 RMAX * 0,0 ID = 1 KA»1 KB»2 T=O.D NA=M( IA I WRITE(6,903I NA F0RMATP1SI30I NB=NA+1 NUMM=«0 NUMD=0 NC=NA+2 HX=1./DFL0AT(NB» HT=HX ZPI = HX * RPI INITILIZE DO 30 IB = 1,NC Y(KA, IB>=DSIN(ZPI*DFLOAT< IB-ll) A=.5*HT/(HX*HX> B=1.+2.*A V=1.-2.*A C=A NUMM=NUMM+4 NUMD = NUMD ♦ 1 E(ll=C.O F(1I = 0.0 ONE STEP OF THE INTEGRATION - SOLVE SYSTEM OF EQUATIONS Y(KA,1) * 0.0 DO 100 IB = 2,NB CPT=B-C*E( IB-ll E(IB»= A/CPT D= A*Y(KA,IB-1)+V*Y(KA,IB)*A*Y(KA,IB*1I F(IB) = (D+C*F( IB-1 H/CPT EXA10000 EXAniOC EXA10200 EXA10300 EXA10400 EXA10500 EXA10600 EXA10700 EXA13800 EXA10900 EXA11000 EXAiliOO EXA11200 EXA11300 EXA11400 EXA11500 EXA11600 EXA11700 EXA11800 EXA11900 EXA120CO EXA12100 EXA12200 EXA12300 EXA12400 EXA12500 EXA12600 EXA12700 EXA12800 EXA12900 EXA13000 EXA13100 EXA13200 EXA13300 EXA13400 EXA13500 EXA13600 EXA13700 EXA13800 EXA13900 EXA14000 EXA14100 EXA14200 EXA14300 EXAU400 EXA14500 EXA14600 EXA14700 EXA14800 EXA14900 EXA15000 EXA15100 EXA15200 EXA15300 121 Y(KB,NC»»0.0 00 2DC IB«1,NA IC*NC-IB 200 Y(KB,ICI»E< IC»»Y(KB, IC + ll+FUCI NUMM»NUMM*6*NA NUMD»NUMD*2*NA KA=3-KA KB* 3-KB T=T*HT IF( T .LT. TTSTCIOII GO TO 50 OUTPUT SECTION 00 300 IB=1,NA VAL=DEXP(-QPI*T»*DSIN HT = HX ZPI = HX * RPI INITILIZE DO 33 IB «= 1,NC D1 K = 'hT/HX OEXP( " rtX * OFLOAT{IB - 1,, * DSIN «^ p I*DFLOAT(IB-in A = . 5D0*D1/HX B = 1.00 ♦ 2.D0 *A C = A Al * A ♦ Dl A2 = l.D0-2.D0*A+HT A3 = A-Dl NUMM=NUMM+4 NUMD = NUMD ♦ 1 E(1» = 0.0 F(l »=0.0 ONE STEP OF THE INTEGRATION - SOLVE SYSTEM OF EQUATIONS Y(KA,1> = 0.0 DO IOC IB * 2,NB CPT*B-C*E( IB-ll EXB10000 EXB101DO EXB10200 EXB10300 EXB10400 EXB10500 EXB106OO EXB10700 EXB10800 EXB1C9G0 EXB11000 EXB11100 EX811200 EXB11300 EXB11400 EXB11500 EXB11630 EXB11700 EXB11800 EXB11900 EXB12000 EXB12100 EXB12200 EXB12300 EXB12400 EXB12500 EXB12600 EXB12700 EXB12800 EXB12900 EXB13000 EXB13100 EXB13200 EXB13300 EXB13400 EX813500 EXB13600 EXB13700 EXB13800 EXB13900 EX814000 EXB14100 EXB14200 EXB14300 EXB14400 EXB14500 EXB14600 EXB147C0 EXB14800 EXB14900 EXB15000 EXB15100 EXB15200 EXB15300 12U 130 200 ♦ A1*Y(KA,IB*1I 900 901 902 904 500 E f I B)« A/CPT D » A3 * Y(KA,IB-1) ♦ A2*Y(KA,IB» F( IB»*(D+C*F(IB-ll»/CPT Y(KB,NC)*0.0 DO 230 IB«1,NA IC=NC-IB Y(KB f IC»»E{ IC»*Y(KB, IC+1KFUCI NUMM=NUMM*6*NA NUMD=NUMD*2*NA KA=3-KA KB=3-KB T=T+HT IF( T .LT. TTST(ID) I GO TO 50 OUTPUT SECTION WRITE (6»900» T,HX,A,8 WRITE (6, 904> Y(KA,IP1+1»,Y(KA,IP2*1I,Y(KA,IP3+1»,Y(KA, IP4+1J, 1 Y(KA,IP5 + 1» ,Y(KA,IP6+1>,Y(KA,IP7+1I F0RMAT(1X,5D26.15) WRITE(6,901) NUMM,NUMD F0RMAT(1X,3I30> WRITE(6,902» FORMAT!* •» F0RMAT(1X,7D18.9I ID = ID ♦ 1 GO TO 50 IF(T I PI = IP2 = IP3 = IP4 = IP5 = IP6 '- IP7 = LT, IP1 IP2 IP3 IP4 IP5 IP6 IP7 161 * 2 CONTINUE CALL EXIT END EXB15400 EXB15500 EXB15600 EXB15700 EXB15800 EXB15900 EXB16000 EXB16100 EXB16200 EXB16300 EXB16400 EXB16500 EXB16600 EXB16700 EXB16800 EXB16900 EXB17000 EXB17100 EXB17200 EXB17300 EXB17400 EXB17500 EXB17600 EXB17700 EXB17800 EXB17900 EXB18000 EX818100 EXB18200 EXB18300 EXB18400 EXB18500 EXB18600 EXB18700 EXB18800 EXB18900 EXB19000 125 APPENDIX D EXAMPLE 3 126 IMPLICIT REAL *8(A-H,0-ZI EXC10000 DIMENSION TTSTI8I EXCIOIOO DIMENSION Y(2,102 8>,M(11),YTST(1028» , RMATI 3, 1028 i ,B { 10281 EXC 10200 DATA M/3, 7, 15, 31, 63, 127, 256, 511, 1023, 2045, 4095/ EXC10300 DATA TTST/. 125, .25, .5,1. ,2. .4,, 8., 16./ EXC10400 CALL UNOERZI»OFF»l EXC10500 RK1 » 1. DO EXC10600 RK2 » l.DO EXC10700 RPI*3. 141592653589793 EXC10800 IP1 * 1 EXC10900 IP2 * 2 EXC11000 IP3 » 3 EXC11100 IP4 = 4 EXC11200 I p 5 * 5 EXC11300 IP6 * 6 EXC1140C IP7 = 7 EXC11500 C EXC11600 C MAIN LOOP FOR DIFFERENT STEP SIZES EXC11700 C EXC11800 DO 530 IA * 2,7 EXC11900 RMAX « 0.0 EXC12000 ID * 1 EXC12100 K **l EXC12200 *B»2 EXC12300 T*0.0 EXC12400 NA=M(IA) EXC12500 WRITE(6,903> NA EXC12600 903 F0RMAT(«1«,I30I EXC12700 NB=NA*1 EXC12800 NUMM=0 EXC12900 NUMD=0 EXC13000 NC=NA+2 EXC13100 HX=1./DFL0AT(NBI EXC13200 QPI = RPI*HX EXC13300 HT=HX EXC13400 c EXC13500 C INITILIZE EXC13600 c EXC13700 DO 3D IB = i,NC EXC13800 V = DSIN(QPI*DFL0AT(IB-1)) EXC13900 Y(KA,IB» = V*V EXC14000 EXC14100 SET UP VALUES FOR MATRIX EXC14200 EXC 14300 R = HT/(HX*HX) EXC14400 Bl = l.ODO ♦ HX EXC14500 B2 = l.ODO ♦ HX EXC14600 TR = 2. DO * R EXC14700 TRB1 = TR*B1 EXC14800 TRB2 = TR*B2 EXC14900 RNATflfll » -R EXC15000 RMAT(2,1» = 2.00 ♦ TRB1 EXC15100 RMAT<3,1» * -TR EXC15200 RMAT(1,NC» » -TR EXC15300 127 RMAT(2,NC» » 2. DO ♦ TRB2 RMATOtNCI- -R RMAIN * 2.D0 ♦ TR RMANE » 2.00 -TR 10 DO ID I « RMAT(1,II RMAT(2,I) RMAT(3,I) CONTINUE 2,N8 » -R » RHAIN « -R DECOMPOSE THE MATRIX CALL DEC0MP B(NC) * TR*YIKA,NBI ♦ ( 2. D0-TRB2 I *Y( KA, NC ) C C SET UP THE RIGHT SIDE OF EQUATION C DO 23 I = 2,NB B(I) = R*Y(KA,I-1I ♦ RMANE*Y(KA,I I ♦ R*Y) 3C0 CONTINUE WRITE (6,9001 T,HX WRITE (6 ,904) Y( KA , I Pl + 1 ) , Y( KA, I P2*l ) , Y( KA , IP3 + 1 ) , Y( KA , I P4 + 1 I . 1 Y(KA,IP54l),Y(KA,IP6*l),Y(KA,IP7U) i??STU ( p6nY!sT!lJ7! Pn,YTST(IP2,,YTST(IP3,,YTST(IPM,YTST(IP5, » 900 F0RMAT(1X,5D26. 151 WRITE(6,90l) NUMM,NUMO 901 F0RMAT(1X,3I30) WRI TE(6,902» 902 FORMAT! • •) 9C4 F0RMAT(1X,7018.9) ID = 10 ♦ 1 IF(T .LT. 161 GO TO 50 IF(RMAX .LT, .001) CALL EXIT EXC 15400 EXC15500 EXC15600 EXC15700 EXC15800 EXC15900 EXC16000 EXC 16100 EXC16200 EXC 16300 EXC 16400 EXC16500 EXC16600 EXC16700 EXC16800 EXC16900 EXC17000 EXC17100 EXC17200 EXC17300 EXC17400 EXC1750C EXC17600 EXC17700 EXC17800 EXC17900 EXC 18000 EXC18100 EXC1820O EXC18300 EXC18400 EXC18500 EXC18600 EXC18700 EXC18800 EXC18900 EXC19000 EXC19100 EXC19200 EXC19300 EXC19400 EXC19500 EXC19600 EXC19700 EXC19800 EXC 19900 EXC20000 EXC20100 EXC20200 EXC20300 EXC20400 EXC20500 EXC20600 EXC23700 EXC20800 IP1 » IP1 * 2 IP2 « IP2 * 2 IP3 * IP3 * 2 IP* . ip* * 2 IP5 * IP5 * 2 IP6 « IP6 * 2 IP7 » IP7 * 2 500 CONTINUE CALL UNDERZ(»ON»> CALL EXIT 200 WRITE(6,2000) 2000 FORMAT! 1 BAD MATRIX 1 » CALL EXIT END 128 EXC 20900 EXC21000 EXC21100 EXC21200 EXC21300 EXC21400 EXC21500 EXC21600 EXC21700 EXC21800 EXC21900 EXC22000 EXC22100 EXC22200 129 VITA Leonard Andrew Larsen was born on April 22, 19UO. He attended elementary and secondary schools in Waukesha, Wisconsin. After graduating from Waukesha High School in 1958, he attended the University of Wisconsin, Madison, where the Bachelor of Science degree in Physics and Mathematics Education was conferred in 1962. Prior to attending the University of Illinois in 1967, he taught mathematics and physics at De Forest High School at De Forest, Wisconsin, for four years and at Conant High School in Hoffman Estates, Illinois, for one year. He received a Master of Science degree in Mathematics from the University of Illinois in August 1968. While continuing his education at the University of Illinois he taught at Illinois State University during the 1969 academic year. Beginning with the 1971 academic year, he has held the position of assistant professor of computer science at the University of Wisconsin, Eau Claire. He is currently an associate member of the Association for Computing Machinery and Sigma Xi . orm AEC-427 (6/68) AECM 3201 U.S. ATOMIC ENERGY COMMISSION UNIVERSITY-TYPE CONTRACTOR'S RECOMMENDATION FOR DISPOSITION OF SCIENTIFIC AND TECHNICAL DOCUMENT ( See Instructions on Reverse Side ) AEC REPORT NO. COO -ll*69 -0212 2. TITLE AUTOMATIC SOLUTIONS OF PARTIAL DIFFERENTIAL EQUATIONS TYPE OF DOCUMENT (Check one): f~1 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 |X] c. Other (Specify) Thesis RECOMMENDED ANNOUNCEMENT AND DISTRIBUTION (Check one): a. AEC's normal announcement and distribution procedures may be followed. □ b. Make available only within AEC and to AEC contractors and other U.S. Government agencies and their contractors. 1 I c. Make no announcement or distrubution. REASON FOR RECOMMENDED RESTRICTIONS: SUBMITTED BY: NAME AND POSITION (Please print or type) C. W. Gear, Professor and Principal Investigator Organization Department of Computer Science University of Illinois Urbana, Illinois 6l801 Signature ^KM^d^f^-- Date October 1972 FOR AEC USE ONLY AEC CONTRACT ADMINISTRATOR'S COMMENTS, IF ANY, ON ABOVE ANNOUNCEMENT AND DISTRIBUTION RECOMMENDATION: PATENT CLEARANCE: □ a. AEC patent clearance has been granted by responsible AEC patent group. □ b. Report has been sent to responsible AEC patent group for clearance. I I c. Patent clearance not required. BIBLIOGRAPHIC DATA SHEET 4. Title and Subtitle 1. Report No. UIUCDCS-R-72-5U6 AUTOMATIC SOLUTIONS OF PARTIAL DIFFERENTIAL EQUATIONS 7. Author(s) Leonard Andrew Lars en 9. Performing Organization Name and Address Department of Computer Science University of Illinois Urbana, Illinois 6l801 12. Sponsoring Organization Name and Address US AEC Chicago Operations Office 9800 South Cass Avenue Argonne, Illinois 15. Supplementary Notes 3. Recipient's Accession No. 5. Report Date October 1972 8. Performing Organization Rept. No. 10. Project/Task/Work Unit No. Ili69-Gear 11. Contract /Grant No. US AEC AT(ll-l)lU69 13. Type of Report & Period Covered Thesis (Ph.D.) 14. 16. Abstracts The problem of the present paper is to look for an automatic method that can be applied to some subset of partial differential equations. The primary value of a program involving an automatic method would be to provide the person who needs only a few runs with a particular type of equation the opportunity to find a solution, within acceptable tolerances, without having to write and debug a specialized program. 17. Key Words and Document Analysis. 17a. Descriptors partial differential equations automatic method spatial discretization 7b. Identifiers/Open-Ended Terms 7e. COSAT1 Field/Group 8. Availability Statement unlimited distribution ORM N TIS-3B ( 10-70) 19. Security Class (This Report) VflCLASSIFIED 20. Security Class (This Page UNCLASSIFIED 21- No. of Pages 129 22. Price USCOMM-DC 40329-P7 1