OS/2 Shareware BBS: 22 gnu/22-gnu.zip/fweb140x.zip/demos/series.web

Text File | 1993-10-30 | 36KB | 1,024 lines

%   REVISION HISTORY
%
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   LIMBO MATERIAL
%
\font\ninerm=cmr9
\let\mc=\ninerm % medium caps for names like PASCAL
\def\WEB{{\tt WEB}}
\def\PASCAL{{\mc PASCAL}}
\def\<{$\langle\,$}
\def\>{$\,\rangle$}
\def\C{{\bf C}}
\def\unix{{\it unix}}
\def\Unix{{\it Unix}}
\def\today{\ifcase\month\or
   January\or February\or March\or April\or
   May\or June\or July\or August\or
   September\or October\or November\or December\fi
   \space\number\day, \number\year}
\def\title{{\bf Power Series Operators and Tests}}
\def\topofcontents{\hsize 6in
  \vglue -30pt plus 1fil minus 1.5in
  \centerline{\title}
  \vskip 15pt
  \centerline{Bart Childs and Timothy McGuire}
  \smallskip
  \centerline{Department of Computer Science}
  \centerline{Texas A{\rm\char'046}M University}
  \centerline{College Station, TX ~ 77843-3112}
  \smallskip
  \centerline{bart@@cs.tamu.edu ~ ~ mcguire@@cs.tamu.edu}
  \centerline{409-845-5470      ~ ~ 409-845-1298}
  \vskip 15pt
  \centerline{\today}
  \vfil
  \def\?##1]{\hbox to 1in{\hfil##1.\ }}}
\def\botofcontents{\vskip 0pt plus 1fil minus 1.5in}
%
%   END OF LIMBO MATERIAL
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   BEGINNING OF WEB
%

\def\BSl{{\rm\char'134}}
\font\sc=cmcsc10
\def\Fortran{{\sc Fortran}}
\def\PS#1{#1_0 + #1_1 t  + \dots + #1_i t^i}
\def\DPS#1{#1_1 + 2 #1_2 t + \dots + i #1_i t^{i-1}}
\font\ninett=cmtt9
\def\WEB{{\ninett WEB}}\let\web=\WEB
\def\FWEB{{\ninett FWEB}}\let\fweb=\FWEB
\def\zeron{\ifmmode 0, \dots\ , n\else $0, \dots\ ,n$\fi}
@n

@* Power Series.  The power series is one of the more common
tools in applied mathematics.  It often comes from the truncated
Taylor or Maclaurin series.  The solution of ordinary differential
equations can be expressed in this form but it often makes for
a lot of work and long expressions.

Power series are an important theoretical tool in mathematics.
They can also be used to great advantage in numerical computations.
In many cases the solution of a differential equation can be
expressed in terms of their power series expansions.
One disadvantage of this representation is that the formal
manipulation of these power series can often be quite complex.

The Frobenius form is an economical means of calculating these
coefficients rather than using the repeated differentiation.
We will describe various routines which will expedite some of
the formal manipulations necessary to solve a fairly broad subset
of the common differential equations.

The first part of the code will be a {\bf driver} or {\bf main}
program that will test these operations.  This is presented first
because it is done in a form of \Fortran, which has that orientation.

@ The main program will be skeletal at this point.  The use of
\WEB{} allows us to write in a {\it pseudo-code} form in both
top-down and bottom-up modes.

This program is written to run under the Cray {\tt cft77} and
the usual compilers for |SUN|s, |VAX| {VMS}, and Data General's |AOS|.
This is controlled by the macros |AOS|, |CRAY|, |SUN|, and |VAX|.
It should be straightforward to accomodate other environments.
@^32-bit 64-bit@>

@ The {\tt ftangle} command line to produce |CRAY| code is:\hfil\break
\qquad {\tt ftangle pst -m"CRAY" -d}\hfil\break
The {\tt -d} converts  |do|/|end do| constructions into labeled |do|
loops. Cray, like everybody else, has extended their compiler beyond
the standard, but they omitted this logical extension.  Our coding
style assumes that the standards extended by MIL-STD 1752 are
present.  This added the |do|/|end do| and |implicit none|
statements to \Fortran.
@^32-bit 64-bit@>

{\vfil\it
This is a prerelease form.  It will appear shortly as a technical
report from the Department of Computer Science, \TeX{}as A{\rm\char'046}M
University.}

@ Charles Karney of Princeton University furnished most of these macros
for machine dependencies.  They came as part of a larger body of code released
with John Krommes' \FWEB{} system.  The default is for |SUN| compilation.
The following use of \FWEB{} features are quite ``{\tt C}'' like.
@^32-bit 64-bit@>

@#if defined(CRAY)
@m AOS 0
@m CRAY 1
@m SUN 0
@m VAX 0
@#elif defined(VAX)
@m AOS 0
@m CRAY 0
@m SUN 0
@m VAX 1
@#elif defined(AOS)
@m AOS 1
@m CRAY 0
@m SUN 0
@m VAX 0
@#else
@m AOS 0
@m CRAY 0
@m SUN 1
@m VAX 0
@#endif


@ We start each \Fortran{} module with a command that defeats
the archaic default typing of variable names.  We don't suggest
that it should be removed from the compiler because of the
need for upward compatability from millions of lines of
existing code.  We question it remaining the default mode.

The Sun uses |implicit undefined(a-z)| while the other target environments
have a rational form, |implicit none|.  We declare |implicit_none|
and let the magic of \fweb{} macros substitute the form the target machine needs.
In this and the next section we have used the \web{} command,
{\tt@@f}.
This instructs the formatting of the words {\tt none}, {\tt undefined},
\dots{} to be done like the reserved word {\bf implicit}.

@f none implicit
@f undefined implicit
@f implicit_none implicit

@#if SUN
@m implicit_none implicit undefined(a-z)
@#else
@m implicit_none implicit none
@#endif

@ We introduce another macro to make it easy to run the tests
with 64~bit precision in all the environments.  The Cray uses
64~bit precision for |real| variables by default.  The other
machines will use this when |real*8| is specified.  However,
the actual storage will vary among most of the machines.

Another macro is used to convert floating point constants
by appending the ``{\tt d0}'' exponent on 32-bit machines.
@^32-bit 64-bit@>

@f floating real
@#if CRAY
  @m floating real
  @m const(x) x
@#else
  @m floating real*8
  @m const(x) x##d0
  @#endif


@ The {\tt main} program will be a top view with only three additional
statements, namely |implicit_none|, |stop| and |end|.  The
declarations will come in several parts and each set of declarations
will be documented in a literate manner.
The appearance of the {\it semicolons} might frighten some
\Fortran{} traditionalists.  We don't apologize because
we have become accustomed to programming statements having
ending markers, like sentences.!?

@a
      program test_ps
      implicit_none
      @<Main program declarations@>
      @<Main initializations@>
      @<Test trigonometric, multiply, and square@>
      @<Test power and square root@>
      @<Test division@>
      @<Test exponential@>
      stop
      end

@ We will declare a |parameter| |n| to be the maximum index of
our power series.  In most problems the power series will have an
infinite range of convergence.  However, the concept of numeric
continuation is important and is similar to the concept of analytic
continuation.  The range of convergence for the usual power series
is based upon the location of singularities.  The power series
has a limited useful range due to the increasing cost
of evaluating each coefficient and the nature of floating point
roundoff.  Thus, each power series will have a maximum
index.  Previous tests have shown that values of approximately 12
give the maximum efficiency for reasonable accuracies.  We will
name this parameter |n| because it will be used so frequently.
@^analytic continuation@>
@^numeric continuation@>
@^continuation@>

@<Main program declarations@>=
      integer n
      parameter (n=12)

@ We will need several arrays for this testing.  The variables
|x| and |y| will be generic (coefficients of) power series.
We will calculate the power series that represent functions of
these generic power series.  Then we need arrays
|sn|, |cs|, |ex|, and |sum|
to be the power series that represent $\sin(x)$,
$\cos(x)$, $\exp(x)$, and $\sin^2(x)+\cos^2(x)$,
respectively.  The argument may be |y| instead of |x|.
We approximate
$$x(t) \approx \sum^n_{i=0} x_i t^i$$

Our arrays will generally have a range of \zeron{} which
means there will be $n+1$ elements.  We will refer to the range
rather than the number of elements to avoid the confusion
caused by some programmers considering the upper limit of
the index to be the ``size'' of the array.  We also declare
an array |unity| which is used for testing purposes.

@<Main program declarations@>=
      floating x(0:n), y(0:n), unity(0:n)
      floating sn(0:n), cs(0:n), ex(0:n)
      floating sum(0:n)

@ We need a local |integer| variable and a |format| array.
The `*' format yields different output on the different
machines.  We use |PS| to get consistent output of the
arrays of power series coefficients across these machines.
The |write(*,*)| is a neat feature of \Fortran~77, but its
results vary and we want the same (at least similar)
results on several different machines.  Thus, we will use
formatted output in several places.

@<Main program declarations@>=
      integer i
      character*80 PS

@ This module is a bottom up part of our code.  We have already
written the set of functions (which start a few modules later) and
we must make ready to incorporate them.

We also need to tell our main program the type of our user
functions.  Some compilers also require that |subroutine|s be
declared external if you use |implicit_none|.  The |subroutine|s
|ps_trig| and |ps_diff| are also used in these tests.

@<Main program declarations@>=
    floating ps_mult, ps_div, ps_exp, ps_sqr, ps_pwr, ps_sqrt
    external ps_mult, ps_div, ps_exp, ps_sqr, ps_pwr, ps_sqrt
        external ps_trig, ps_diff

@ The power series |unity| and |x| are really quite simple.
The represent |unity|$(t)\equiv1$ and |x|$(t) \equiv t$.
We use the series |x| because the sine,
cosine, and exponential coefficients for this power series are
obvious.  This ``visual check'' is a part of our testing.
This is also the module where we initialize the |format|
array |PS|.

@<Main initializations@>=
      PS = '(1x,5E15.7/(3x5E15.7))'
      do  i = 0, n
        unity(i) = const(0.0)
        x(i) = const(0.0)
        end do
      unity(0) = const(1.0)
      x(1) = const(1.0)

@ The trigonometric identity of $\sin^2(x)+\cos^2(x)=1$ is a simple
test for these routines.  We expect the results will be close,
within roundoff error. Since the exponential
function is related to the trigonometric functions, we will also
calculate it for this specific form of $x(t)$.
We will output the results of the calculations for visual
inspection.  Notice that we called a |subroutine| to calculate the
trigonometric functions and that we use |function|s for
multiplication, squaring, and exponentiation.  The reasons for
this are more clear when |ps_trig| is developed.

@<Test trigonometric, multiply, and square@>=
      do i = 0, n
        call ps_trig(x,sn,cs,i)
        sum(i)= ps_sqr(sn,i) + ps_mult(cs,cs,i)
        ex(i) = ps_exp(x,i,ex)
        end do
      write(*,*)
      write(*,*) "Basic functions, sin, cos, exp:"
      write(*,PS) sn
      write(*,PS) cs
      write(*,PS) ex

@ We could also output the series |sum|.  However, we know that
it should be close to the series |unity|.  We have developed a
|subroutine|, |ps_diff| for such comparisons.  We will call this
rather than outputting the array.

@<Test trigonometric, multiply, and square@>=
      write(*,*)
      write(*,*) " Comparing the trig identity with unity"
      call ps_diff( unity, sum, n )

@ We now are ready to test power and square root.  The generic series
|x| has a null value for its $0^{\hbox{\sevenrm th}}$ coefficient.
We will define the
generic series |y| to be $y=\exp(x)$ to avoid problems that could
arise from the null $0^{\hbox{\sevenrm th}}$ coefficient.

@<Test power and square root@>=
      do i=0,n
        y(i)=ps_exp(x,i,y)
        end do

@ We need two more arrays for testing the power and square root
functions.  The first will contain the series that is the square
of the original function and the second will contain the square
root of the square, resulting in a complete circle, we hope.

@<Main program declarations@>=
      floating square(0:n), sqrroot(0:n)

@ Here we use the obvious identity $ y = \sqrt{(y^2)}$.
We expect that we will be close, within roundoff error.
@^32-bit 64-bit@>

@<Test power and square root@>=
      do i = 0, n
        square(i)=ps_pwr(y,const(2.0),i,square)
        sqrroot(i)=ps_sqrt(square,i,sqrroot)
        end do
      write(*,*)
      write(*,*) "Comparing the square root of y^2.0 to y"
      call ps_diff( y, sqrroot, n )

@ We now are ready to test |ps_div|. The declaration is
obvious and the test is designed to be as simple as possible
without using a trivial array.

@<Main program declarations@>=
      floating quotient(0:n)

@ Here we use the obvious identity $ y = y*y/y $.  We will
multiply and divide as shown.

@<Test division@>=
     do i = 0, n
        square(i) = ps_sqr(y,i)
        quotient(i)=ps_div(square,y,i,quotient)
        end do
      write(*,*)
      write(*,*) "Comparing (y^2)/y to y"
      call ps_diff( y, quotient, n )

@ Lastly, we want to test the exponential routine.  We will again make use
of |x|, and we will need the additional arrays |minus_y|, |ex_inv|,
and |prod| to be the power series that represent $-y$, $e^{-y}$, and
$e^y e^{-y}$, respectively.  The array, |ex|, was declared in
the first test.

@<Main program declarations@>=
    floating minus_y(0:n)
    floating ex_inv(0:n), prod(0:n)

@ We are using the obvious identity of $e^y e^{-y} = 1$.  This should be the
result we get.  Our expected results would be this simple series with
some or all terms affected slightly by roundoff error.  The results
will be worse with 32~bit arithmetic.

@<Test exponential@>=
      do i = 0, n
    minus_y(i) = -y(i)
        ex(i)=ps_exp(y,i,ex)
    ex_inv(i)=ps_exp(minus_y,i,ex_inv)
    prod(i)=ps_mult(ex,ex_inv,i)
        end do
      write(*,*)
      write(*,*) "Comparing exp(x)*exp(-x) with unity"
      call ps_diff( unity, prod, n )

@* Comparing test results.
We mentioned the utility |ps_diff| which is used to compare power
series coefficients.  This |subroutine| has three arguments.  The
first is the series that is the answer, the second is the
calculated series, and the third is the maximum index of both series.

@a
      subroutine ps_diff( a, c, n)
        implicit_none
        @<Difference variables@>
        @<Difference initialization@>
        do i = 0, n
          if( a(i) != const(0.0) ) then
            @<Non zero checks@>
          else
            @<Zero checks@>
            end if
          end do
        @<Output differences@>
        return
        end

@ We will calculate the relative difference of array elements
that should be non-zero and absolute value of elements that
should have been zero.  We only report the maximum values of
each.

The variable |i| is a dummy loop variable.
We use the variables |in| and |iz| to point to these elements.
The value of |nz| will be the number of zero elements in
the array.  We don't output a message about zero coefficients
if there aren't any.
The variable |rn| has the maximum relative difference of a
non-zero element.

        @<Difference variables@>=
        integer n
        floating a(0:n), c(0:n)
        integer i, iz, in, nz
        floating rn

@ The initializations are rather straightforward.  We use
the negative one to have a convenient signal that we may
have had no problems.

        @<Difference initialization@>=
        in = -1
        iz = -1
        rn = const(0.0)
        nz = 0

@ If a value is non-zero, then we will calculate the relative
difference by dividing the absolute value of the difference
by the absolute value of the original.

        @<Non zero checks@>=
        if( abs(a(i)-c(i))/abs(a(i)) > rn ) then
          rn = abs(a(i)-c(i))/abs(a(i))
          in = i
          end if

@ The checking for those elements that should be zero is simple.
The first time we remember it, the rest of the time we have to
check to make sure it is larger.

        @<Zero checks@>=
        nz = nz + 1
        if( nz == 1) then
          iz = i
        else
          if( abs(c(i)) > abs(c(iz)) ) iz = i
          end if

@ The output will consist of one line for the non-zero coefficients
and if there is at least one zero coefficient, one line for its output.

        @<Output differences@>=
        if( in != -1 ) then
          write(*,"("" Maximum relative change, index"",G15.7,I3)") rn,in
        else
          write(*,*) " Non zero coefficients are identical."
          end if
        if( nz > 0 ) then
          if( iz != -1 ) then
            write(*,"("" Worst from 0.0, index "",G15.7,I3)") c(iz), iz
          else
            write(*,"("" All "",I2,"" zeroes remained true."")") nz
            end if
          end if

@* Power Series Operations.  The basic arithmetic operations on power
series (addition, subtraction, multiplication, and division) are
relatively straightforward.  Addition  and subtraction are trivial,
as they involve only adding and subtracting the appropriate coefficients.
Multiplication is likewise trivial, but involves enough computation
to make a function call worthwhile.  Division is quite a bit more
complicated, but still can be computed by a means which is equivalent
to that of the longhand division technique for polynomials taught in
elementary algebra.

Other operations to be performed are the exponential, sine and cosine,
square root, and arbitrary power of a power series.  We will explain in
detail how these functions are calculated.

It should be stressed here that in all the operations described, the
function or subroutine computes only one term of the series each time it
is called.  Although this may seem counter-intuitive at first, it is this
feature which allows us to be able to use these functions efficiently in
the complicated application in which we will see them.

Again, we will consistently approximate:
$$x(t) \approx \sum^n_{i=0} x_i t^i$$
This power series, and many others, will be developed by the use
of Frobenius recurrence relations.  Of course, we may use $a(t)$,
$b(t)$, and a wide array of symbols for these truncated power series.

@* Multiplication.  We first consider computing the product of two power
series.  The inputs to this function are two vectors |a| and |b|,
with indices \zeron. The product of these is
$$\sum_{k = 0}^n c_k t^k    = \left ( \sum_{k = 0}^n a_k t^k \right )
\left ( \sum_{k = 0}^n b_k t^k  \right )$$
where $c_k$ expands to $$c_k = a_0 b_k + a_1 b_{k-1} + \dots + a_k b_0.$$

The values of elements $a_0, \dots , a_i$ and $b_0, \dots , b_i$ are known
at the time of the call.  The function result, |ps_mult|, is the term $c_i$
defined above.

@a
   floating function ps_mult(a,b,i)
    implicit_none
    @<Declare |ps_mult| variables@>
    @<Compute |ps_mult| result@>
   return
   end

@ We will declare the input parameters to the function.
The arrays |a| and |b| probably have an actual range of \zeron{}.
We declare a subset range of $0, \dots\ , i$.
We must also declare a single local loop index.
\itemitem{|i|} is the index of the product term to be computed.
\itemitem{|a|} and |b| are as described above.
\itemitem{|k|} is a local variable used as a loop index.

    @<Declare |ps_mult| variables@>=
    integer i
    floating a(0:i), b(0:i)
    integer k


@  The computation is a rather straightforward implementation of the above
formula for $c_k$, using a loop to sum the result.  We note that the
saving of one addition will not work in \Fortran~66.  That does
not count because it did not have zero indices anyway.

      @<Compute |ps_mult| result@>=
      ps_mult = a(0)*b(i)
    do k=1,i
      ps_mult = ps_mult + a(k)*b(i-k)
      end do


@* Division. To compute the quotient of two power series, we use an adaptation
of longhand division of polynomials (with some ``neat tricks'' to help
simplify the computation.)  We compute the quotient $q$ as follows:
$$q(t) = {{\PS{u}} \over {\PS{d}}}$$

If we attempt to perform the long division symbolically at this point, we
might quickly despair at the successively increasing complexity of each
coefficient.  However, we observe that each successive coefficient of the
quotient series depends on the previous coefficients. Simplifying, we obtain
the following formula:
$$ q_i = {u_i - \sum_{k=0}^{i-1} q_k d_{i-k} \over d_0}$$

The following function will return the coefficient of the $i^{th}$ term of
of the quotient $q$, where $u$ is the numerator and $d$ is the denominator.
The quotient $q$ is also an argument because the $q_0, \dots , q_{i-1}$
are needed to compute $q_i$.

@a
   floating function ps_div( u, d, i, q )
    implicit_none
    @<Declare |ps_div| variables@>
    @<Compute |ps_div| result@>
    return
    end


@ We will declare the parameters to the function.  |u| and |d| are
vectors with indices of \zeron. The vector of quotient coefficients,
|q|, must be also be sent, as described above.
We must declare a single local loop index.
\itemitem{|i|} is the index of the quotient term to be computed.
\itemitem{|u|} is the numerator, |d| is the denominator, and |q|
is the quotient.  All of these have an upper limit of |i|.
\itemitem{|k|} is a local variable used as a loop index.
    @<Declare |ps_div| variables@>=
    integer i
          floating u(0:i), d(0:i), q(0:i)
    integer k

@ The computation is a straightforward application of the above formula for
$q_i$, using a loop to perform the required sum.  The function
result is summed in |ps_div| and copied into $q_i$ .

    @<Compute |ps_div| result@>=
    ps_div = u(i)
    do  k = 0, i-1
          ps_div = ps_div - q(k) * d(i-k)
      end do
    ps_div = ps_div/d(0)
    q(i)   = ps_div

@* Exponential.  To compute the exponential function of a power series it is
necessary to use our knowledge of elementary calculus. We derive the series
as follows:   If
$$e(t) = \exp(a(t))$$
then by differentiating (and omitting the $(t)$ of each term) we obtain
$$\dot e = \exp(a) \dot a = e \dot a$$
or
$$ \DPS{e} = (\PS{e}) (\DPS{a}) $$
If we equate the coefficients of $t^{i-1}$ on either side we
obtain the relation:
$$i e_i = i a_i e_0 + (i-1) a_{i-1} e_1 + \dots + a_1 e_{i-1}$$
or
$$e_i = {\sum_{k=1}^i k a_k e_{i-k} \over i}$$
It is clear that $e_0 = \exp (a_0)$ from the definition of $e(t)$.
The $k a_k$ term is a manifestation of $\dot a(t)$.

@a
      floating function ps_exp( a, i, e)
        implicit_none
    @<Declare |ps_exp| variables@>
    if( i == 0)  then
      @<Perform the trivial case for |ps_exp|@>
        else
          @<Perform the general case for |ps_exp|@>
          end if
    e(i) = ps_exp
    return
    end

@  We declare the parameters for the |ps_exp|.  The argument of the
the exponential function is |a|, |i| is the current
term that we are calculating, and |e| is the vector of the previous results
(|e| needs to be  passed because the previous coefficients are used in
computing the next one.)  We also declare a single local loop index.
\itemitem{|i|} is the index of the exponential term to be computed.
\itemitem{|a|} and |e| are the zero-indexed vectors described above.  Their
upper limit is |i|.
\itemitem{|k|} is a local variable used as a loop index.
    @<Declare |ps_exp| variables@>=
    integer i
        floating  a(0:i), e(0:i)
        integer k

@ The $0^{th}$ term is the trivial case, where $e_0 = \exp (a_0)$.

    @<Perform the trivial case for |ps_exp|@>=
    ps_exp = exp(a(0))

@ The general case is a straightforward implementation of the formula for
$e_i$ described above.

    @<Perform the general case for |ps_exp|@>=
    ps_exp = a(1)*e(i-1)
    do  k = 2, i
      ps_exp = ps_exp + k * a(k) * e(i-k)
      end do
    ps_exp = ps_exp/i


@* Trigonometric Functions. To compute the trigonometric functions of a
power series we will use a ``trick'' similar to that of the exponential
function. We derive the sine and  cosine of a power series as follows:
If
$$s = \sin(a)$$
and
$$c = \cos(a)$$
then by differentiating, we obtain
$$\dot s = \cos(a) \dot a = c \dot a$$
and
$$\dot c = - \sin(a) \dot a = - s \dot a$$
Expanded, these become
$$ \DPS{s} = (\PS{c})(\DPS{a})$$
and
$$ \DPS{c} = - (\PS{s})(\DPS{a})$$

If we equate the coefficients of $t^{i-1}$ on either side we
obtain the recurrences:
$$i s_i = i a_i c_0 + (i-1) a_{i-1} c_1 + \dots + a_1 c_{i-1}$$
and
$$i c_i = i a_i s_0 + (i-1) a_{i-1} s_1 + \dots + a_1 s_{i-1}$$
with
$$s_0 = \sin(a_0)$$
and
$$c_0 = \cos(a_0)$$

It is thus convenient to calculate the sine and cosine series
together, although only one of them might be present in the equation.


@ Now for the \.{FORTRAN} code.  This will be a subroutine because we
actually calculate (and in a sense return) two values with each call,
namely $s(i)$ and $c(i)$.


@a
       subroutine ps_trig(a,s,c,i)
       implicit_none
       @<Declare |ps_trig| variables@>
       if (i==0) then
           @<Do the trivial case for |ps_trig|@>
       else
           @<Do the general case for |ps_trig|@>
         end if
       return
       end

@  The arguments must be declared.  They are: |a|, the argument of the
sine and cosine functions; |s| and |c|, the trig functions; and |k|, the
index of the coefficients being calculated.  The first three of these are
power series and therefore have a lowest subscript of {\bf zero}.  |k| is
used as a loop control variable.

       @<Declare |ps_trig| variables@>=
       integer i
       floating a(0:i), s(0:i), c(0:i)

       integer k

@  When calculating the zeroth coefficient, we use the library functions.

         @<Do the trivial case for |ps_trig|@>=
           s(0) = sin(a(0))
           c(0) = cos(a(0))

@  The general coefficient, {\it i.e.}  when $i \not = 0$, requires a
multiplication of power series operators.  Remember that we are calculating
sin and cos both because they are dependent upon each other.

         @<Do the general case for |ps_trig|@>=
     s(i) = c(i-1) * a(1)
     c(i) = s(i-1) * a(1)
     do k=2,i
       s(i) = s(i) + c(i-k) * (k*a(k))
       c(i) = c(i) + s(i-k) * (k*a(k))
       end do
     s(i) = s(i)/i
      c(i) = -c(i)/i



@* Square.  The next function we will need for the power series operations
is the function to square the vector.  It should be noted that the square
function is simply a special case of multiplication.  We have optimized it
here for reasons of efficiency.

We pair the identical terms in the sum so that they need not be computed for a
second time.  The modification of the recurrence formula can be expressed
as follows:
$$c_i = \cases {(a_0)^2 &if $i=0$\cr
      2 \sum\limits_{k=0}^{i-1 \over 2} a_k a_{i-k} &if $i>0$ and odd\cr
(a_{i \over 2})^2
       + 2 \sum\limits_{k=0}^{\lfloor{i-1 \over 2} \rfloor}
        a_k a_{i-k} &if $i>0$ and even \cr}
$$

@a
      floating function ps_sqr(a,i)
    implicit_none
    @<Declare |ps_sqr| variables@>
    if (i == 0) then
      ps_sqr = a(0)*a(0)
    else
      @<Compute |ps_sqr| sum@>
      @<Check for an odd number of terms and update |ps_sqr|@>
      end if
        return
        end

@ We will declare the input parameters to the function.  The vector
|a| has indices \zeron{} and is input to the function.
We must also declare a single local loop index.
\itemitem{|i|} is the index of the square term to be computed.
\itemitem{|a|} is the input vector of power series coefficients.
\itemitem{|k|} is a local variable used as a loop index.

           @<Declare |ps_sqr| variables@>=
    integer i
    floating a(0:i)
    integer k


@ The computation is a rather straightforward adaptation of the |ps_mult|
computation.  It is made more efficient by pairing the identical terms in
the sum so that they are computed only once.  The |if| statement is
necessary around the |do| to prevent the execution of the |do| loop when
computing the coefficient of the $0^{th}$ term.
@^32-bit 64-bit@>

    @<Compute |ps_sqr| sum@>=
    ps_sqr = a(0)*a(i)
     do k=1,(i-1)/2
      ps_sqr = ps_sqr + a(k)*a(i-k)
      end do
    ps_sqr = const(2.0)*ps_sqr

@ We now need to check for an odd number of terms.  If so the middle term
has not been added to |ps_sqr| and we need to add it to our sum.

    @<Check for an odd number of terms and update |ps_sqr|@>=
    if (mod(i,2) == 0) then
      ps_sqr = ps_sqr + a(i/2)**2
      end if


@* Power.  To compute the coefficients $r_i$ of the power series resulting
from raising a power series $a$ to an arbitrary real power $p$, we write
$$ r = a^p$$
Differentiation yields
$$ \dot r =  p a^{p-1} \dot a$$
After multiplying both sides by $a$ we obtain
$$ {\dot r} a = {p r} {\dot a}$$
Expanded out, this is
$$\displaylines{\quad(\DPS{r})(\PS{a}) \hfill\cr
\hfill= p (\PS{r})(\DPS{a})\quad\cr}$$
Thus, by equating coefficients for $t^{i-1}$ we have
$$ \sum_{k=1}^i k r_k a_{i-k} = p \sum_{k=1}^i k a_k r_{i-k}$$
Solving this we can obtain the coefficients $r_i$ by the recurrence
$$ 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})
           \over i a_0}$$

@a
      floating function ps_pwr( a, p, i, r )
    @<Declare variables for |ps_pwr|@>
    if( i == 0 ) then
      @<Perform the trivial case for |ps_pwr|@>
    else
      @<Perform the general case for |ps_pwr|@>
      end if
    r(i) = ps_pwr
        return
        end

@ We will declare the parameters to the function.  |a| is our input power
series, |p| is the power to which we wish to raise |a|,  |i| is the current
term that we are calculating, and |r| is the vector of previously computed
coefficients (|r| needs to be passed because the previous coefficients are
used in computing the next one.)  We also declare a single local loop index
|k|.

    @<Declare variables for |ps_pwr|@>=
    implicit_none
    integer i
    floating a(0:i), p, r(0:i)

    integer k

@ the $0^{th}$ term is the trivial case, where $r_0 = (a_0)^p$.

    @<Perform the trivial case for |ps_pwr|@>=
      ps_pwr = a(0)**p

@ Now, we can perform the calculation for all other cases by the recurrence
formula derived above.

    @<Perform the general case for |ps_pwr|@>=
      ps_pwr = p * i * a(i) * r(0)
      do k = 1, i-1
        ps_pwr = ps_pwr + (k)*( p*a(k)*r(i-k) - r(k)*a(i-k) )
            end do
      ps_pwr = ps_pwr / ( i * a(0) )




@* Square Root.  As a final example, we will compute the square root |r| of a
power series |a|.  (We could use the |ps_pwr| function with |p| = 0.5, but
this will be more efficient.) The derivation is reasonably straightforward:
$$r = \sqrt {a}$$ is equivalent to $$r^2 = a$$ or
$$(\PS{r})^2 = \PS{a}$$
By equating coefficients on $t^i$  we obtain
$$r_0 r_i + r_1 r_{i-1} + \dots + r_i r_0 = a_i$$
Solving for $r_i$ yields the recurrence
$$r_i = {1 \over 2 r_0} \left [ a_i - \sum_{k=1}^{i-1} r_k r_{i-k} \right ]$$
with $r_0 = \sqrt a_0$.  The implementation takes
advantage of the symmetry of the sum, like it was done in the
|ps_square| routine.

@a
          floating function ps_sqrt( a, i, r)
    @<Declare |ps_sqrt| variables@>
    if( i == 0 ) then
      @<Perform the {\it zero} case for |ps_sqrt|@>
    else
      @<Perform the general case for |ps_sqrt|@>
          end if
    r(i) = ps_sqrt
    return
    end

@ We will declare the parameters to the function. |a| is the zero-indexed
vector of coefficients for the power series we are taking the square root
of, |i| is the current term that we are calculating, and |r| is the vector
of previous results (which needs to be passed because the previous
coefficients are used in computing the next one.)  We also declare a
single local loop index.
\itemitem{|i|} is the index of the square root term to be computed.
\itemitem{|a|} and |r| are the zero-indexed vectors described above and
their upper limit is |i|.
\itemitem{|k|} is a local variable used as a loop index.

    @<Declare |ps_sqrt| variables@>=
    implicit_none
          integer i
    floating a(0:i),  r(0:i)
    integer k


@ The $0^{th}$ term is a trivial case, where $r_0 = \sqrt a_0$.

      @<Perform the {\it zero} case for |ps_sqrt|@>=
      ps_sqrt = sqrt( a(0) )

@ The general case is not quite straightforward. We have a special form
when $i=1$ and have to check for the odd central term because we
avoid summing both ends of the sum to take advantage of the symmetry.

@^32-bit 64-bit@>

      @<Perform the general case for |ps_sqrt|@>=
          if(i==1) then
            ps_sqrt = const(0.5)*a(1)/r(0)
          else
            ps_sqrt = const(0.0)
            do k = 1, (i-1)/2
          ps_sqrt = ps_sqrt + r(k)*r(i-k)
              end do
            if( mod(i,2)==0) then
              ps_sqrt = ps_sqrt + const(0.5)*r(i/2)**2
              end if
            ps_sqrt = ( const(0.5)*a(i) - ps_sqrt) / r(0)
            end if

@* Additional work.  There are a few more operations that are needed
for some applications.
They may not be added in the next month or so, but they are planned.
These include |ps_shift| which would allow the multiplication or
division of a power series by the independent variable to an
arbitrary integral power, like in Bessel's equations.
Of course, this can be done by defining that operand in terms of
the |ps_mult| or |ps_div| operations and using a series (like |x| in
the tests above).  That would not give maximum efficiency and can
cause some more problems when considering singularities.

Another routine more in line with those already shown is
integration, like in integro-differential equations.
This has been done but we want to let it settle in out minds
and use it in several samples before releasing its format and
implied need of maintenance for a reasonable future.

Gibbons also included a |reciprocal| routine.  We have chosen
not to implement it because our target is differential equations
and such terms are rare in the ones we use.  Gibbons included
a |log| routine which we probably need but have not implemented.

\noindent
June 1, 1990.

@* References. The following references are available to some extent
and show a number of uses and background work.  They are annotated
to a slight extent and additional references are welcome.

{\parindent0sp \everypar{\hangindent2pc\hangafter1}

Doiron, H. H.  Numerical Integration via Power Series Expansions,
    M.S. Thesis, University of Houston, Houston, Texas, August 1967.
{\it The primary results in this thesis show that these integration
methods are fast.  These and other tests usually show that the
use of these methods saves approximately 80 percent of the CPU time.
It uses some archaic data structures rather to do some of these
functions (like sine and cosine).  The algorithms here are superior.}

Fehlberg, E.  Numerical Integration of Differential Equations by Power
    Series Expansions, Illustrated by Physical Examples.
    NASA Technical Note No. TN D-2356, October 1964.
{\it The two examples are a restricted three-body problem and the
motion of an electron in the field of a magnetic dipole.  The results
indicate a required CPU time of 15 to 20 percent of the
Runge-Kutta-Nystr\"om method.}

Gibbons, A.  A Program for the Automatic Integration of Differential
    Equations using the Method of Taylor Series, Computer Journal,
    vol. 3, pp. 108-111  (1960).
{\it A good source of algorithms and discussion of the procedures.
The term ``practical radius of convergence'' is introduced which
is like the ``numeric continuation'' term used in this paper.}

Henrici, P. Automatic Computations with Power Series,
    {\sl Journal of the ACM}, Vol.~3, no.~1, Jan.~1956.
    {\it An early reference of doing similar work in a symbolic
    fashion. His operations are not as efficient as those described in
    this paper. Applications included are combinatorial analysis,
    asymptotic expansions, computing Legendre polynomials, and
    solving differential equations.}

}

There are many other references to the use of power series as
integrators but they don't really affect the spirit of this
work.  We will include them in the next revision.  Two primary
authors of these works are listed but I did not have time
to get the complete references.  I just did not want the impression
to get out that we were not aware of their works.

{\parindent0sp \everypar{\hangindent2pc\hangafter1}

Change

Corliss

}

@* Index.