home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 22 gnu / 22-gnu.zip / fweb140x.zip / demos / series.web < prev    next >
Text File  |  1993-10-30  |  36KB  |  1,024 lines

  1. %   REVISION HISTORY        
  2. %
  3. %
  4. %
  5. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  6. %
  7. %   LIMBO MATERIAL
  8. %
  9. \font\ninerm=cmr9
  10. \let\mc=\ninerm % medium caps for names like PASCAL
  11. \def\WEB{{\tt WEB}}
  12. \def\PASCAL{{\mc PASCAL}}
  13. \def\<{$\langle\,$}
  14. \def\>{$\,\rangle$}
  15. \def\C{{\bf C}}
  16. \def\unix{{\it unix}}
  17. \def\Unix{{\it Unix}}
  18. \def\today{\ifcase\month\or
  19.    January\or February\or March\or April\or
  20.    May\or June\or July\or August\or
  21.    September\or October\or November\or December\fi
  22.    \space\number\day, \number\year}
  23. \def\title{{\bf Power Series Operators and Tests}}
  24. \def\topofcontents{\hsize 6in
  25.   \vglue -30pt plus 1fil minus 1.5in
  26.   \centerline{\title}
  27.   \vskip 15pt
  28.   \centerline{Bart Childs and Timothy McGuire}
  29.   \smallskip
  30.   \centerline{Department of Computer Science}
  31.   \centerline{Texas A{\rm\char'046}M University}
  32.   \centerline{College Station, TX ~ 77843-3112}
  33.   \smallskip
  34.   \centerline{bart@@cs.tamu.edu ~ ~ mcguire@@cs.tamu.edu}
  35.   \centerline{409-845-5470      ~ ~ 409-845-1298}
  36.   \vskip 15pt
  37.   \centerline{\today}
  38.   \vfil
  39.   \def\?##1]{\hbox to 1in{\hfil##1.\ }}}
  40. \def\botofcontents{\vskip 0pt plus 1fil minus 1.5in}
  41. %
  42. %   END OF LIMBO MATERIAL
  43. %
  44. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  45. %
  46. %   BEGINNING OF WEB
  47. %
  48.  
  49. \def\BSl{{\rm\char'134}}
  50. \font\sc=cmcsc10
  51. \def\Fortran{{\sc Fortran}}
  52. \def\PS#1{#1_0 + #1_1 t  + \dots + #1_i t^i}
  53. \def\DPS#1{#1_1 + 2 #1_2 t + \dots + i #1_i t^{i-1}}
  54. \font\ninett=cmtt9
  55. \def\WEB{{\ninett WEB}}\let\web=\WEB
  56. \def\FWEB{{\ninett FWEB}}\let\fweb=\FWEB
  57. \def\zeron{\ifmmode 0, \dots\ , n\else $0, \dots\ ,n$\fi}
  58. @n
  59.  
  60. @* Power Series.  The power series is one of the more common
  61. tools in applied mathematics.  It often comes from the truncated
  62. Taylor or Maclaurin series.  The solution of ordinary differential
  63. equations can be expressed in this form but it often makes for
  64. a lot of work and long expressions.
  65.  
  66. Power series are an important theoretical tool in mathematics.
  67. They can also be used to great advantage in numerical computations.
  68. In many cases the solution of a differential equation can be
  69. expressed in terms of their power series expansions.
  70. One disadvantage of this representation is that the formal
  71. manipulation of these power series can often be quite complex.
  72.  
  73. The Frobenius form is an economical means of calculating these
  74. coefficients rather than using the repeated differentiation.
  75. We will describe various routines which will expedite some of
  76. the formal manipulations necessary to solve a fairly broad subset
  77. of the common differential equations.
  78.  
  79. The first part of the code will be a {\bf driver} or {\bf main}
  80. program that will test these operations.  This is presented first
  81. because it is done in a form of \Fortran, which has that orientation.
  82.  
  83. @ The main program will be skeletal at this point.  The use of
  84. \WEB{} allows us to write in a {\it pseudo-code} form in both
  85. top-down and bottom-up modes.
  86.  
  87. This program is written to run under the Cray {\tt cft77} and
  88. the usual compilers for |SUN|s, |VAX| {VMS}, and Data General's |AOS|.
  89. This is controlled by the macros |AOS|, |CRAY|, |SUN|, and |VAX|.
  90. It should be straightforward to accomodate other environments. 
  91. @^32-bit 64-bit@>
  92.  
  93. @ The {\tt ftangle} command line to produce |CRAY| code is:\hfil\break
  94. \qquad {\tt ftangle pst -m"CRAY" -d}\hfil\break
  95. The {\tt -d} converts  |do|/|end do| constructions into labeled |do|
  96. loops. Cray, like everybody else, has extended their compiler beyond
  97. the standard, but they omitted this logical extension.  Our coding
  98. style assumes that the standards extended by MIL-STD 1752 are
  99. present.  This added the |do|/|end do| and |implicit none|
  100. statements to \Fortran.
  101. @^32-bit 64-bit@>
  102.  
  103. {\vfil\it
  104. This is a prerelease form.  It will appear shortly as a technical
  105. report from the Department of Computer Science, \TeX{}as A{\rm\char'046}M
  106. University.}
  107.  
  108. @ Charles Karney of Princeton University furnished most of these macros
  109. for machine dependencies.  They came as part of a larger body of code released
  110. with John Krommes' \FWEB{} system.  The default is for |SUN| compilation.
  111. The following use of \FWEB{} features are quite ``{\tt C}'' like.
  112. @^32-bit 64-bit@>
  113.  
  114. @#if defined(CRAY)
  115. @m AOS 0
  116. @m CRAY 1
  117. @m SUN 0
  118. @m VAX 0
  119. @#elif defined(VAX)
  120. @m AOS 0
  121. @m CRAY 0
  122. @m SUN 0 
  123. @m VAX 1 
  124. @#elif defined(AOS)
  125. @m AOS 1
  126. @m CRAY 0
  127. @m SUN 0 
  128. @m VAX 0
  129. @#else
  130. @m AOS 0
  131. @m CRAY 0
  132. @m SUN 1 
  133. @m VAX 0 
  134. @#endif
  135.  
  136.  
  137. @ We start each \Fortran{} module with a command that defeats
  138. the archaic default typing of variable names.  We don't suggest
  139. that it should be removed from the compiler because of the
  140. need for upward compatability from millions of lines of
  141. existing code.  We question it remaining the default mode.
  142.  
  143. The Sun uses |implicit undefined(a-z)| while the other target environments
  144. have a rational form, |implicit none|.  We declare |implicit_none|
  145. and let the magic of \fweb{} macros substitute the form the target machine needs.
  146. In this and the next section we have used the \web{} command,
  147. {\tt@@f}.
  148. This instructs the formatting of the words {\tt none}, {\tt undefined},
  149. \dots{} to be done like the reserved word {\bf implicit}.
  150.  
  151. @f none implicit
  152. @f undefined implicit
  153. @f implicit_none implicit
  154.  
  155. @#if SUN
  156. @m implicit_none implicit undefined(a-z)
  157. @#else
  158. @m implicit_none implicit none
  159. @#endif
  160.  
  161. @ We introduce another macro to make it easy to run the tests
  162. with 64~bit precision in all the environments.  The Cray uses
  163. 64~bit precision for |real| variables by default.  The other
  164. machines will use this when |real*8| is specified.  However,
  165. the actual storage will vary among most of the machines.
  166.  
  167. Another macro is used to convert floating point constants
  168. by appending the ``{\tt d0}'' exponent on 32-bit machines.
  169. @^32-bit 64-bit@>
  170.  
  171. @f floating real
  172. @#if CRAY
  173.   @m floating real
  174.   @m const(x) x
  175. @#else
  176.   @m floating real*8
  177.   @m const(x) x##d0
  178.   @#endif
  179.  
  180.  
  181. @ The {\tt main} program will be a top view with only three additional
  182. statements, namely |implicit_none|, |stop| and |end|.  The
  183. declarations will come in several parts and each set of declarations
  184. will be documented in a literate manner.
  185. The appearance of the {\it semicolons} might frighten some
  186. \Fortran{} traditionalists.  We don't apologize because
  187. we have become accustomed to programming statements having
  188. ending markers, like sentences.!?
  189.  
  190. @a
  191.       program test_ps
  192.       implicit_none
  193.       @<Main program declarations@>
  194.       @<Main initializations@>
  195.       @<Test trigonometric, multiply, and square@>
  196.       @<Test power and square root@>
  197.       @<Test division@>
  198.       @<Test exponential@>
  199.       stop
  200.       end
  201.             
  202. @ We will declare a |parameter| |n| to be the maximum index of
  203. our power series.  In most problems the power series will have an
  204. infinite range of convergence.  However, the concept of numeric
  205. continuation is important and is similar to the concept of analytic
  206. continuation.  The range of convergence for the usual power series
  207. is based upon the location of singularities.  The power series
  208. has a limited useful range due to the increasing cost
  209. of evaluating each coefficient and the nature of floating point
  210. roundoff.  Thus, each power series will have a maximum
  211. index.  Previous tests have shown that values of approximately 12
  212. give the maximum efficiency for reasonable accuracies.  We will
  213. name this parameter |n| because it will be used so frequently.
  214. @^analytic continuation@>
  215. @^numeric continuation@>
  216. @^continuation@>
  217.  
  218. @<Main program declarations@>=
  219.       integer n
  220.       parameter (n=12)
  221.  
  222. @ We will need several arrays for this testing.  The variables
  223. |x| and |y| will be generic (coefficients of) power series.
  224. We will calculate the power series that represent functions of
  225. these generic power series.  Then we need arrays
  226. |sn|, |cs|, |ex|, and |sum|
  227. to be the power series that represent $\sin(x)$,
  228. $\cos(x)$, $\exp(x)$, and $\sin^2(x)+\cos^2(x)$,
  229. respectively.  The argument may be |y| instead of |x|.
  230. We approximate
  231. $$x(t) \approx \sum^n_{i=0} x_i t^i$$
  232.  
  233. Our arrays will generally have a range of \zeron{} which
  234. means there will be $n+1$ elements.  We will refer to the range
  235. rather than the number of elements to avoid the confusion
  236. caused by some programmers considering the upper limit of
  237. the index to be the ``size'' of the array.  We also declare
  238. an array |unity| which is used for testing purposes.
  239.  
  240. @<Main program declarations@>=
  241.       floating x(0:n), y(0:n), unity(0:n)
  242.       floating sn(0:n), cs(0:n), ex(0:n)
  243.       floating sum(0:n)
  244.       
  245. @ We need a local |integer| variable and a |format| array.
  246. The `*' format yields different output on the different
  247. machines.  We use |PS| to get consistent output of the
  248. arrays of power series coefficients across these machines.
  249. The |write(*,*)| is a neat feature of \Fortran~77, but its
  250. results vary and we want the same (at least similar)
  251. results on several different machines.  Thus, we will use
  252. formatted output in several places.
  253.  
  254. @<Main program declarations@>=
  255.       integer i
  256.       character*80 PS
  257.  
  258. @ This module is a bottom up part of our code.  We have already
  259. written the set of functions (which start a few modules later) and
  260. we must make ready to incorporate them.
  261.  
  262. We also need to tell our main program the type of our user
  263. functions.  Some compilers also require that |subroutine|s be
  264. declared external if you use |implicit_none|.  The |subroutine|s
  265. |ps_trig| and |ps_diff| are also used in these tests.
  266.     
  267. @<Main program declarations@>=
  268.     floating ps_mult, ps_div, ps_exp, ps_sqr, ps_pwr, ps_sqrt
  269.     external ps_mult, ps_div, ps_exp, ps_sqr, ps_pwr, ps_sqrt
  270.         external ps_trig, ps_diff
  271.  
  272. @ The power series |unity| and |x| are really quite simple.
  273. The represent |unity|$(t)\equiv1$ and |x|$(t) \equiv t$.
  274. We use the series |x| because the sine,
  275. cosine, and exponential coefficients for this power series are
  276. obvious.  This ``visual check'' is a part of our testing.
  277. This is also the module where we initialize the |format|
  278. array |PS|.
  279.  
  280. @<Main initializations@>=
  281.       PS = '(1x,5E15.7/(3x5E15.7))'
  282.       do  i = 0, n
  283.         unity(i) = const(0.0)
  284.         x(i) = const(0.0)
  285.         end do
  286.       unity(0) = const(1.0)
  287.       x(1) = const(1.0)
  288.    
  289. @ The trigonometric identity of $\sin^2(x)+\cos^2(x)=1$ is a simple
  290. test for these routines.  We expect the results will be close,
  291. within roundoff error. Since the exponential
  292. function is related to the trigonometric functions, we will also
  293. calculate it for this specific form of $x(t)$.
  294. We will output the results of the calculations for visual
  295. inspection.  Notice that we called a |subroutine| to calculate the
  296. trigonometric functions and that we use |function|s for
  297. multiplication, squaring, and exponentiation.  The reasons for
  298. this are more clear when |ps_trig| is developed.
  299.  
  300. @<Test trigonometric, multiply, and square@>=
  301.       do i = 0, n
  302.         call ps_trig(x,sn,cs,i)
  303.         sum(i)= ps_sqr(sn,i) + ps_mult(cs,cs,i)
  304.         ex(i) = ps_exp(x,i,ex)
  305.         end do
  306.       write(*,*)
  307.       write(*,*) "Basic functions, sin, cos, exp:"
  308.       write(*,PS) sn
  309.       write(*,PS) cs
  310.       write(*,PS) ex
  311.  
  312. @ We could also output the series |sum|.  However, we know that
  313. it should be close to the series |unity|.  We have developed a
  314. |subroutine|, |ps_diff| for such comparisons.  We will call this
  315. rather than outputting the array.
  316.  
  317. @<Test trigonometric, multiply, and square@>=
  318.       write(*,*)
  319.       write(*,*) " Comparing the trig identity with unity"
  320.       call ps_diff( unity, sum, n )
  321.  
  322. @ We now are ready to test power and square root.  The generic series
  323. |x| has a null value for its $0^{\hbox{\sevenrm th}}$ coefficient.
  324. We will define the
  325. generic series |y| to be $y=\exp(x)$ to avoid problems that could
  326. arise from the null $0^{\hbox{\sevenrm th}}$ coefficient.
  327.  
  328. @<Test power and square root@>=
  329.       do i=0,n
  330.         y(i)=ps_exp(x,i,y)
  331.         end do
  332.  
  333. @ We need two more arrays for testing the power and square root
  334. functions.  The first will contain the series that is the square
  335. of the original function and the second will contain the square
  336. root of the square, resulting in a complete circle, we hope.
  337.  
  338. @<Main program declarations@>=
  339.       floating square(0:n), sqrroot(0:n)
  340.  
  341. @ Here we use the obvious identity $ y = \sqrt{(y^2)}$.
  342. We expect that we will be close, within roundoff error.
  343. @^32-bit 64-bit@>
  344.  
  345. @<Test power and square root@>=
  346.       do i = 0, n
  347.         square(i)=ps_pwr(y,const(2.0),i,square)
  348.         sqrroot(i)=ps_sqrt(square,i,sqrroot)
  349.         end do
  350.       write(*,*)
  351.       write(*,*) "Comparing the square root of y^2.0 to y"
  352.       call ps_diff( y, sqrroot, n )
  353.  
  354. @ We now are ready to test |ps_div|. The declaration is
  355. obvious and the test is designed to be as simple as possible
  356. without using a trivial array.
  357.  
  358. @<Main program declarations@>=
  359.       floating quotient(0:n)
  360.  
  361. @ Here we use the obvious identity $ y = y*y/y $.  We will
  362. multiply and divide as shown.
  363.  
  364. @<Test division@>=
  365.      do i = 0, n
  366.         square(i) = ps_sqr(y,i)
  367.         quotient(i)=ps_div(square,y,i,quotient)
  368.         end do
  369.       write(*,*)
  370.       write(*,*) "Comparing (y^2)/y to y"
  371.       call ps_diff( y, quotient, n )
  372.       
  373. @ Lastly, we want to test the exponential routine.  We will again make use
  374. of |x|, and we will need the additional arrays |minus_y|, |ex_inv|, 
  375. and |prod| to be the power series that represent $-y$, $e^{-y}$, and 
  376. $e^y e^{-y}$, respectively.  The array, |ex|, was declared in
  377. the first test.
  378.  
  379. @<Main program declarations@>=
  380.     floating minus_y(0:n)
  381.     floating ex_inv(0:n), prod(0:n)
  382.  
  383. @ We are using the obvious identity of $e^y e^{-y} = 1$.  This should be the
  384. result we get.  Our expected results would be this simple series with
  385. some or all terms affected slightly by roundoff error.  The results
  386. will be worse with 32~bit arithmetic.
  387.  
  388. @<Test exponential@>=
  389.       do i = 0, n
  390.     minus_y(i) = -y(i)
  391.         ex(i)=ps_exp(y,i,ex)
  392.     ex_inv(i)=ps_exp(minus_y,i,ex_inv)
  393.     prod(i)=ps_mult(ex,ex_inv,i)
  394.         end do
  395.       write(*,*)
  396.       write(*,*) "Comparing exp(x)*exp(-x) with unity"
  397.       call ps_diff( unity, prod, n )
  398.  
  399. @* Comparing test results.
  400. We mentioned the utility |ps_diff| which is used to compare power
  401. series coefficients.  This |subroutine| has three arguments.  The
  402. first is the series that is the answer, the second is the
  403. calculated series, and the third is the maximum index of both series.
  404.  
  405. @a
  406.       subroutine ps_diff( a, c, n)
  407.         implicit_none
  408.         @<Difference variables@>
  409.         @<Difference initialization@>
  410.         do i = 0, n
  411.           if( a(i) != const(0.0) ) then
  412.             @<Non zero checks@>
  413.           else
  414.             @<Zero checks@>
  415.             end if
  416.           end do
  417.         @<Output differences@>
  418.         return
  419.         end
  420.         
  421. @ We will calculate the relative difference of array elements
  422. that should be non-zero and absolute value of elements that
  423. should have been zero.  We only report the maximum values of
  424. each.
  425.  
  426. The variable |i| is a dummy loop variable.
  427. We use the variables |in| and |iz| to point to these elements.
  428. The value of |nz| will be the number of zero elements in
  429. the array.  We don't output a message about zero coefficients
  430. if there aren't any.
  431. The variable |rn| has the maximum relative difference of a
  432. non-zero element.
  433.  
  434.         @<Difference variables@>=
  435.         integer n
  436.         floating a(0:n), c(0:n)
  437.         integer i, iz, in, nz
  438.         floating rn
  439.         
  440. @ The initializations are rather straightforward.  We use
  441. the negative one to have a convenient signal that we may
  442. have had no problems.
  443.  
  444.         @<Difference initialization@>=
  445.         in = -1
  446.         iz = -1
  447.         rn = const(0.0)
  448.         nz = 0
  449.         
  450. @ If a value is non-zero, then we will calculate the relative
  451. difference by dividing the absolute value of the difference
  452. by the absolute value of the original.
  453.  
  454.         @<Non zero checks@>=
  455.         if( abs(a(i)-c(i))/abs(a(i)) > rn ) then
  456.           rn = abs(a(i)-c(i))/abs(a(i))
  457.           in = i
  458.           end if
  459.  
  460. @ The checking for those elements that should be zero is simple.
  461. The first time we remember it, the rest of the time we have to
  462. check to make sure it is larger.
  463.  
  464.         @<Zero checks@>=
  465.         nz = nz + 1
  466.         if( nz == 1) then
  467.           iz = i
  468.         else
  469.           if( abs(c(i)) > abs(c(iz)) ) iz = i
  470.           end if
  471.  
  472. @ The output will consist of one line for the non-zero coefficients
  473. and if there is at least one zero coefficient, one line for its output.
  474.  
  475.         @<Output differences@>=
  476.         if( in != -1 ) then
  477.           write(*,"("" Maximum relative change, index"",G15.7,I3)") rn,in
  478.         else
  479.           write(*,*) " Non zero coefficients are identical."
  480.           end if
  481.         if( nz > 0 ) then
  482.           if( iz != -1 ) then
  483.             write(*,"("" Worst from 0.0, index "",G15.7,I3)") c(iz), iz
  484.           else
  485.             write(*,"("" All "",I2,"" zeroes remained true."")") nz
  486.             end if
  487.           end if
  488.  
  489. @* Power Series Operations.  The basic arithmetic operations on power
  490. series (addition, subtraction, multiplication, and division) are 
  491. relatively straightforward.  Addition  and subtraction are trivial, 
  492. as they involve only adding and subtracting the appropriate coefficients. 
  493. Multiplication is likewise trivial, but involves enough computation 
  494. to make a function call worthwhile.  Division is quite a bit more
  495. complicated, but still can be computed by a means which is equivalent 
  496. to that of the longhand division technique for polynomials taught in 
  497. elementary algebra.  
  498.  
  499. Other operations to be performed are the exponential, sine and cosine, 
  500. square root, and arbitrary power of a power series.  We will explain in 
  501. detail how these functions are calculated.
  502.  
  503. It should be stressed here that in all the operations described, the 
  504. function or subroutine computes only one term of the series each time it 
  505. is called.  Although this may seem counter-intuitive at first, it is this 
  506. feature which allows us to be able to use these functions efficiently in 
  507. the complicated application in which we will see them.
  508.  
  509. Again, we will consistently approximate:
  510. $$x(t) \approx \sum^n_{i=0} x_i t^i$$
  511. This power series, and many others, will be developed by the use
  512. of Frobenius recurrence relations.  Of course, we may use $a(t)$,
  513. $b(t)$, and a wide array of symbols for these truncated power series.
  514.  
  515. @* Multiplication.  We first consider computing the product of two power 
  516. series.  The inputs to this function are two vectors |a| and |b|,
  517. with indices \zeron. The product of these is 
  518. $$\sum_{k = 0}^n c_k t^k    = \left ( \sum_{k = 0}^n a_k t^k \right ) 
  519. \left ( \sum_{k = 0}^n b_k t^k  \right )$$
  520. where $c_k$ expands to $$c_k = a_0 b_k + a_1 b_{k-1} + \dots + a_k b_0.$$
  521.  
  522. The values of elements $a_0, \dots , a_i$ and $b_0, \dots , b_i$ are known 
  523. at the time of the call.  The function result, |ps_mult|, is the term $c_i$
  524. defined above. 
  525.  
  526. @a
  527.    floating function ps_mult(a,b,i)
  528.     implicit_none
  529.     @<Declare |ps_mult| variables@>
  530.     @<Compute |ps_mult| result@>
  531.    return
  532.    end
  533.  
  534. @ We will declare the input parameters to the function.  
  535. The arrays |a| and |b| probably have an actual range of \zeron{}.
  536. We declare a subset range of $0, \dots\ , i$.
  537. We must also declare a single local loop index. 
  538. \itemitem{|i|} is the index of the product term to be computed.
  539. \itemitem{|a|} and |b| are as described above.
  540. \itemitem{|k|} is a local variable used as a loop index.
  541.  
  542.     @<Declare |ps_mult| variables@>=
  543.     integer i
  544.     floating a(0:i), b(0:i)
  545.     integer k
  546.  
  547.  
  548. @  The computation is a rather straightforward implementation of the above 
  549. formula for $c_k$, using a loop to sum the result.  We note that the
  550. saving of one addition will not work in \Fortran~66.  That does
  551. not count because it did not have zero indices anyway.
  552.   
  553.       @<Compute |ps_mult| result@>=
  554.       ps_mult = a(0)*b(i)
  555.     do k=1,i
  556.       ps_mult = ps_mult + a(k)*b(i-k)
  557.       end do
  558.  
  559.  
  560. @* Division. To compute the quotient of two power series, we use an adaptation
  561. of longhand division of polynomials (with some ``neat tricks'' to help
  562. simplify the computation.)  We compute the quotient $q$ as follows:
  563. $$q(t) = {{\PS{u}} \over {\PS{d}}}$$
  564.  
  565. If we attempt to perform the long division symbolically at this point, we
  566. might quickly despair at the successively increasing complexity of each
  567. coefficient.  However, we observe that each successive coefficient of the
  568. quotient series depends on the previous coefficients. Simplifying, we obtain 
  569. the following formula:  
  570. $$ q_i = {u_i - \sum_{k=0}^{i-1} q_k d_{i-k} \over d_0}$$
  571.  
  572. The following function will return the coefficient of the $i^{th}$ term of
  573. of the quotient $q$, where $u$ is the numerator and $d$ is the denominator.
  574. The quotient $q$ is also an argument because the $q_0, \dots , q_{i-1}$
  575. are needed to compute $q_i$.
  576.  
  577. @a
  578.    floating function ps_div( u, d, i, q )
  579.     implicit_none
  580.     @<Declare |ps_div| variables@>
  581.     @<Compute |ps_div| result@>
  582.     return
  583.     end
  584.  
  585.  
  586. @ We will declare the parameters to the function.  |u| and |d| are 
  587. vectors with indices of \zeron. The vector of quotient coefficients,
  588. |q|, must be also be sent, as described above.   
  589. We must declare a single local loop index. 
  590. \itemitem{|i|} is the index of the quotient term to be computed.
  591. \itemitem{|u|} is the numerator, |d| is the denominator, and |q| 
  592. is the quotient.  All of these have an upper limit of |i|.
  593. \itemitem{|k|} is a local variable used as a loop index.
  594.     @<Declare |ps_div| variables@>=
  595.     integer i   
  596.           floating u(0:i), d(0:i), q(0:i) 
  597.     integer k    
  598.  
  599. @ The computation is a straightforward application of the above formula for 
  600. $q_i$, using a loop to perform the required sum.  The function 
  601. result is summed in |ps_div| and copied into $q_i$ .
  602.   
  603.     @<Compute |ps_div| result@>=
  604.     ps_div = u(i)  
  605.     do  k = 0, i-1
  606.           ps_div = ps_div - q(k) * d(i-k)
  607.       end do
  608.     ps_div = ps_div/d(0)
  609.     q(i)   = ps_div        
  610.  
  611. @* Exponential.  To compute the exponential function of a power series it is  
  612. necessary to use our knowledge of elementary calculus. We derive the series 
  613. as follows:   If 
  614. $$e(t) = \exp(a(t))$$ 
  615. then by differentiating (and omitting the $(t)$ of each term) we obtain 
  616. $$\dot e = \exp(a) \dot a = e \dot a$$
  617. or
  618. $$ \DPS{e} = (\PS{e}) (\DPS{a}) $$
  619. If we equate the coefficients of $t^{i-1}$ on either side we 
  620. obtain the relation: 
  621. $$i e_i = i a_i e_0 + (i-1) a_{i-1} e_1 + \dots + a_1 e_{i-1}$$
  622. or
  623. $$e_i = {\sum_{k=1}^i k a_k e_{i-k} \over i}$$
  624. It is clear that $e_0 = \exp (a_0)$ from the definition of $e(t)$.
  625. The $k a_k$ term is a manifestation of $\dot a(t)$.
  626.  
  627. @a
  628.       floating function ps_exp( a, i, e)
  629.         implicit_none
  630.     @<Declare |ps_exp| variables@>
  631.     if( i == 0)  then
  632.       @<Perform the trivial case for |ps_exp|@>
  633.         else
  634.           @<Perform the general case for |ps_exp|@>
  635.           end if
  636.     e(i) = ps_exp
  637.     return
  638.     end
  639.  
  640. @  We declare the parameters for the |ps_exp|.  The argument of the
  641. the exponential function is |a|, |i| is the current 
  642. term that we are calculating, and |e| is the vector of the previous results 
  643. (|e| needs to be  passed because the previous coefficients are used in 
  644. computing the next one.)  We also declare a single local loop index.
  645. \itemitem{|i|} is the index of the exponential term to be computed.
  646. \itemitem{|a|} and |e| are the zero-indexed vectors described above.  Their
  647. upper limit is |i|. 
  648. \itemitem{|k|} is a local variable used as a loop index.
  649.     @<Declare |ps_exp| variables@>=
  650.     integer i
  651.         floating  a(0:i), e(0:i)  
  652.         integer k  
  653.  
  654. @ The $0^{th}$ term is the trivial case, where $e_0 = \exp (a_0)$.  
  655.  
  656.     @<Perform the trivial case for |ps_exp|@>=
  657.     ps_exp = exp(a(0))
  658.  
  659. @ The general case is a straightforward implementation of the formula for
  660. $e_i$ described above.
  661.  
  662.     @<Perform the general case for |ps_exp|@>=
  663.     ps_exp = a(1)*e(i-1)
  664.     do  k = 2, i
  665.       ps_exp = ps_exp + k * a(k) * e(i-k)
  666.       end do
  667.     ps_exp = ps_exp/i
  668.  
  669.  
  670. @* Trigonometric Functions. To compute the trigonometric functions of a 
  671. power series we will use a ``trick'' similar to that of the exponential 
  672. function. We derive the sine and  cosine of a power series as follows:
  673. If 
  674. $$s = \sin(a)$$ 
  675. and
  676. $$c = \cos(a)$$ 
  677. then by differentiating, we obtain 
  678. $$\dot s = \cos(a) \dot a = c \dot a$$
  679. and
  680. $$\dot c = - \sin(a) \dot a = - s \dot a$$
  681. Expanded, these become
  682. $$ \DPS{s} = (\PS{c})(\DPS{a})$$
  683. and
  684. $$ \DPS{c} = - (\PS{s})(\DPS{a})$$
  685.  
  686. If we equate the coefficients of $t^{i-1}$ on either side we 
  687. obtain the recurrences:
  688. $$i s_i = i a_i c_0 + (i-1) a_{i-1} c_1 + \dots + a_1 c_{i-1}$$
  689. and
  690. $$i c_i = i a_i s_0 + (i-1) a_{i-1} s_1 + \dots + a_1 s_{i-1}$$
  691. with
  692. $$s_0 = \sin(a_0)$$
  693. and 
  694. $$c_0 = \cos(a_0)$$
  695.  
  696. It is thus convenient to calculate the sine and cosine series 
  697. together, although only one of them might be present in the equation.
  698.  
  699.    
  700. @ Now for the \.{FORTRAN} code.  This will be a subroutine because we
  701. actually calculate (and in a sense return) two values with each call,
  702. namely $s(i)$ and $c(i)$.
  703.  
  704.  
  705. @a
  706.        subroutine ps_trig(a,s,c,i)
  707.        implicit_none
  708.        @<Declare |ps_trig| variables@>
  709.        if (i==0) then
  710.            @<Do the trivial case for |ps_trig|@>
  711.        else
  712.            @<Do the general case for |ps_trig|@>
  713.          end if
  714.        return
  715.        end
  716.  
  717. @  The arguments must be declared.  They are: |a|, the argument of the
  718. sine and cosine functions; |s| and |c|, the trig functions; and |k|, the
  719. index of the coefficients being calculated.  The first three of these are
  720. power series and therefore have a lowest subscript of {\bf zero}.  |k| is 
  721. used as a loop control variable.
  722.  
  723.        @<Declare |ps_trig| variables@>=
  724.        integer i
  725.        floating a(0:i), s(0:i), c(0:i)
  726.  
  727.        integer k
  728.  
  729. @  When calculating the zeroth coefficient, we use the library functions.
  730.  
  731.          @<Do the trivial case for |ps_trig|@>=
  732.            s(0) = sin(a(0))
  733.            c(0) = cos(a(0))
  734.  
  735. @  The general coefficient, {\it i.e.}  when $i \not = 0$, requires a 
  736. multiplication of power series operators.  Remember that we are calculating
  737. sin and cos both because they are dependent upon each other.
  738.  
  739.          @<Do the general case for |ps_trig|@>=
  740.      s(i) = c(i-1) * a(1)
  741.      c(i) = s(i-1) * a(1)
  742.      do k=2,i
  743.        s(i) = s(i) + c(i-k) * (k*a(k))
  744.        c(i) = c(i) + s(i-k) * (k*a(k))
  745.        end do
  746.      s(i) = s(i)/i
  747.       c(i) = -c(i)/i
  748.  
  749.  
  750.  
  751. @* Square.  The next function we will need for the power series operations
  752. is the function to square the vector.  It should be noted that the square
  753. function is simply a special case of multiplication.  We have optimized it 
  754. here for reasons of efficiency.
  755.  
  756. We pair the identical terms in the sum so that they need not be computed for a
  757. second time.  The modification of the recurrence formula can be expressed
  758. as follows:
  759. $$c_i = \cases {(a_0)^2 &if $i=0$\cr
  760.       2 \sum\limits_{k=0}^{i-1 \over 2} a_k a_{i-k} &if $i>0$ and odd\cr
  761. (a_{i \over 2})^2
  762.        + 2 \sum\limits_{k=0}^{\lfloor{i-1 \over 2} \rfloor}
  763.         a_k a_{i-k} &if $i>0$ and even \cr}
  764. $$
  765.  
  766. @a
  767.       floating function ps_sqr(a,i)
  768.     implicit_none
  769.     @<Declare |ps_sqr| variables@>
  770.     if (i == 0) then
  771.       ps_sqr = a(0)*a(0)
  772.     else    
  773.       @<Compute |ps_sqr| sum@>
  774.       @<Check for an odd number of terms and update |ps_sqr|@>
  775.       end if   
  776.         return
  777.         end
  778.  
  779. @ We will declare the input parameters to the function.  The vector
  780. |a| has indices \zeron{} and is input to the function.
  781. We must also declare a single local loop index.
  782. \itemitem{|i|} is the index of the square term to be computed.
  783. \itemitem{|a|} is the input vector of power series coefficients.
  784. \itemitem{|k|} is a local variable used as a loop index.
  785.  
  786.            @<Declare |ps_sqr| variables@>=
  787.     integer i  
  788.     floating a(0:i)
  789.     integer k 
  790.  
  791.  
  792. @ The computation is a rather straightforward adaptation of the |ps_mult| 
  793. computation.  It is made more efficient by pairing the identical terms in 
  794. the sum so that they are computed only once.  The |if| statement is 
  795. necessary around the |do| to prevent the execution of the |do| loop when
  796. computing the coefficient of the $0^{th}$ term.
  797. @^32-bit 64-bit@>
  798.  
  799.     @<Compute |ps_sqr| sum@>=
  800.     ps_sqr = a(0)*a(i)
  801.      do k=1,(i-1)/2
  802.       ps_sqr = ps_sqr + a(k)*a(i-k)
  803.       end do
  804.     ps_sqr = const(2.0)*ps_sqr
  805.  
  806. @ We now need to check for an odd number of terms.  If so the middle term 
  807. has not been added to |ps_sqr| and we need to add it to our sum.
  808.  
  809.     @<Check for an odd number of terms and update |ps_sqr|@>=
  810.     if (mod(i,2) == 0) then
  811.       ps_sqr = ps_sqr + a(i/2)**2
  812.       end if
  813.  
  814.  
  815. @* Power.  To compute the coefficients $r_i$ of the power series resulting
  816. from raising a power series $a$ to an arbitrary real power $p$, we write
  817. $$ r = a^p$$  
  818. Differentiation yields
  819. $$ \dot r =  p a^{p-1} \dot a$$
  820. After multiplying both sides by $a$ we obtain
  821. $$ {\dot r} a = {p r} {\dot a}$$
  822. Expanded out, this is
  823. $$\displaylines{\quad(\DPS{r})(\PS{a}) \hfill\cr
  824. \hfill= p (\PS{r})(\DPS{a})\quad\cr}$$ 
  825. Thus, by equating coefficients for $t^{i-1}$ we have
  826. $$ \sum_{k=1}^i k r_k a_{i-k} = p \sum_{k=1}^i k a_k r_{i-k}$$
  827. Solving this we can obtain the coefficients $r_i$ by the recurrence
  828. $$ r_i = {i p a_i r_0 + \sum_{k=1}^{i-1} k (p a_k r_{i-k} - r_k a_{i-k})
  829.            \over i a_0}$$
  830.  
  831. @a
  832.       floating function ps_pwr( a, p, i, r )
  833.     @<Declare variables for |ps_pwr|@>
  834.     if( i == 0 ) then
  835.       @<Perform the trivial case for |ps_pwr|@>
  836.     else
  837.       @<Perform the general case for |ps_pwr|@>
  838.       end if
  839.     r(i) = ps_pwr
  840.         return
  841.         end
  842.  
  843. @ We will declare the parameters to the function.  |a| is our input power
  844. series, |p| is the power to which we wish to raise |a|,  |i| is the current
  845. term that we are calculating, and |r| is the vector of previously computed
  846. coefficients (|r| needs to be passed because the previous coefficients are
  847. used in computing the next one.)  We also declare a single local loop index
  848. |k|.
  849.  
  850.     @<Declare variables for |ps_pwr|@>=
  851.     implicit_none
  852.     integer i
  853.     floating a(0:i), p, r(0:i)
  854.  
  855.     integer k
  856.  
  857. @ the $0^{th}$ term is the trivial case, where $r_0 = (a_0)^p$.
  858.  
  859.     @<Perform the trivial case for |ps_pwr|@>=
  860.       ps_pwr = a(0)**p
  861.  
  862. @ Now, we can perform the calculation for all other cases by the recurrence
  863. formula derived above.
  864.  
  865.     @<Perform the general case for |ps_pwr|@>=
  866.       ps_pwr = p * i * a(i) * r(0)
  867.       do k = 1, i-1
  868.         ps_pwr = ps_pwr + (k)*( p*a(k)*r(i-k) - r(k)*a(i-k) )
  869.             end do
  870.       ps_pwr = ps_pwr / ( i * a(0) )
  871.  
  872.  
  873.  
  874.  
  875. @* Square Root.  As a final example, we will compute the square root |r| of a 
  876. power series |a|.  (We could use the |ps_pwr| function with |p| = 0.5, but 
  877. this will be more efficient.) The derivation is reasonably straightforward:
  878. $$r = \sqrt {a}$$ is equivalent to $$r^2 = a$$ or
  879. $$(\PS{r})^2 = \PS{a}$$
  880. By equating coefficients on $t^i$  we obtain
  881. $$r_0 r_i + r_1 r_{i-1} + \dots + r_i r_0 = a_i$$
  882. Solving for $r_i$ yields the recurrence 
  883. $$r_i = {1 \over 2 r_0} \left [ a_i - \sum_{k=1}^{i-1} r_k r_{i-k} \right ]$$
  884. with $r_0 = \sqrt a_0$.  The implementation takes
  885. advantage of the symmetry of the sum, like it was done in the
  886. |ps_square| routine.
  887.  
  888. @a
  889.           floating function ps_sqrt( a, i, r)
  890.     @<Declare |ps_sqrt| variables@>
  891.     if( i == 0 ) then
  892.       @<Perform the {\it zero} case for |ps_sqrt|@>
  893.     else
  894.       @<Perform the general case for |ps_sqrt|@>
  895.           end if
  896.     r(i) = ps_sqrt
  897.     return
  898.     end
  899.  
  900. @ We will declare the parameters to the function. |a| is the zero-indexed 
  901. vector of coefficients for the power series we are taking the square root
  902. of, |i| is the current term that we are calculating, and |r| is the vector 
  903. of previous results (which needs to be passed because the previous 
  904. coefficients are used in computing the next one.)  We also declare a 
  905. single local loop index.
  906. \itemitem{|i|} is the index of the square root term to be computed.
  907. \itemitem{|a|} and |r| are the zero-indexed vectors described above and 
  908. their upper limit is |i|.
  909. \itemitem{|k|} is a local variable used as a loop index.
  910.  
  911.     @<Declare |ps_sqrt| variables@>=
  912.     implicit_none
  913.           integer i 
  914.     floating a(0:i),  r(0:i)
  915.     integer k 
  916.  
  917.  
  918. @ The $0^{th}$ term is a trivial case, where $r_0 = \sqrt a_0$.
  919.  
  920.       @<Perform the {\it zero} case for |ps_sqrt|@>=
  921.       ps_sqrt = sqrt( a(0) )
  922.  
  923. @ The general case is not quite straightforward. We have a special form
  924. when $i=1$ and have to check for the odd central term because we
  925. avoid summing both ends of the sum to take advantage of the symmetry.
  926.  
  927. @^32-bit 64-bit@>
  928.  
  929.       @<Perform the general case for |ps_sqrt|@>=
  930.           if(i==1) then
  931.             ps_sqrt = const(0.5)*a(1)/r(0)
  932.           else
  933.             ps_sqrt = const(0.0)
  934.             do k = 1, (i-1)/2
  935.           ps_sqrt = ps_sqrt + r(k)*r(i-k)
  936.               end do
  937.             if( mod(i,2)==0) then
  938.               ps_sqrt = ps_sqrt + const(0.5)*r(i/2)**2
  939.               end if
  940.             ps_sqrt = ( const(0.5)*a(i) - ps_sqrt) / r(0)
  941.             end if
  942.  
  943. @* Additional work.  There are a few more operations that are needed
  944. for some applications.
  945. They may not be added in the next month or so, but they are planned.
  946. These include |ps_shift| which would allow the multiplication or
  947. division of a power series by the independent variable to an
  948. arbitrary integral power, like in Bessel's equations.
  949. Of course, this can be done by defining that operand in terms of
  950. the |ps_mult| or |ps_div| operations and using a series (like |x| in
  951. the tests above).  That would not give maximum efficiency and can
  952. cause some more problems when considering singularities.
  953.  
  954. Another routine more in line with those already shown is
  955. integration, like in integro-differential equations.
  956. This has been done but we want to let it settle in out minds
  957. and use it in several samples before releasing its format and
  958. implied need of maintenance for a reasonable future.
  959.  
  960. Gibbons also included a |reciprocal| routine.  We have chosen
  961. not to implement it because our target is differential equations
  962. and such terms are rare in the ones we use.  Gibbons included
  963. a |log| routine which we probably need but have not implemented.
  964.  
  965. \noindent
  966. June 1, 1990.
  967.  
  968. @* References. The following references are available to some extent
  969. and show a number of uses and background work.  They are annotated
  970. to a slight extent and additional references are welcome.
  971.  
  972. {\parindent0sp \everypar{\hangindent2pc\hangafter1}
  973.  
  974. Doiron, H. H.  Numerical Integration via Power Series Expansions,
  975.     M.S. Thesis, University of Houston, Houston, Texas, August 1967.
  976. {\it The primary results in this thesis show that these integration
  977. methods are fast.  These and other tests usually show that the
  978. use of these methods saves approximately 80 percent of the CPU time.
  979. It uses some archaic data structures rather to do some of these
  980. functions (like sine and cosine).  The algorithms here are superior.}
  981.  
  982. Fehlberg, E.  Numerical Integration of Differential Equations by Power
  983.     Series Expansions, Illustrated by Physical Examples.  
  984.     NASA Technical Note No. TN D-2356, October 1964.
  985. {\it The two examples are a restricted three-body problem and the
  986. motion of an electron in the field of a magnetic dipole.  The results
  987. indicate a required CPU time of 15 to 20 percent of the
  988. Runge-Kutta-Nystr\"om method.}
  989.  
  990. Gibbons, A.  A Program for the Automatic Integration of Differential 
  991.     Equations using the Method of Taylor Series, Computer Journal, 
  992.     vol. 3, pp. 108-111  (1960).
  993. {\it A good source of algorithms and discussion of the procedures.
  994. The term ``practical radius of convergence'' is introduced which
  995. is like the ``numeric continuation'' term used in this paper.}
  996.  
  997. Henrici, P. Automatic Computations with Power Series,
  998.     {\sl Journal of the ACM}, Vol.~3, no.~1, Jan.~1956.
  999.     {\it An early reference of doing similar work in a symbolic
  1000.     fashion. His operations are not as efficient as those described in
  1001.     this paper. Applications included are combinatorial analysis,
  1002.     asymptotic expansions, computing Legendre polynomials, and
  1003.     solving differential equations.}
  1004.     
  1005. }
  1006.  
  1007. There are many other references to the use of power series as
  1008. integrators but they don't really affect the spirit of this
  1009. work.  We will include them in the next revision.  Two primary
  1010. authors of these works are listed but I did not have time
  1011. to get the complete references.  I just did not want the impression
  1012. to get out that we were not aware of their works.
  1013.  
  1014. {\parindent0sp \everypar{\hangindent2pc\hangafter1}
  1015.  
  1016. Change
  1017.  
  1018. Corliss
  1019.  
  1020. }
  1021.  
  1022. @* Index.
  1023.  
  1024.