| Up | Next | Prev | PrevTail | Tail |
This package implements formal Laurent power series expansions in one variable using the domain mechanism of REDUCE. This means that power series objects can be added, multiplied, differentiated etc. like other first class objects in the system. A lazy evaluation scheme is used in the package and thus terms of the series are not evaluated until they are required for printing or for use in calculating terms in other power series. The series are extendible giving the user the impression that the full infinite series is being manipulated. The errors that can sometimes occur using series that are truncated at some fixed depth (for example when a term in the required series depends on terms of an intermediate series beyond the truncation depth) are thus avoided.
The package was originally based on an earlier truncated power series package developed by Julian Padget in the 1980’s. The name of the original package was TPS and this was never changed. The alternative (more accurate) name EPS was perhaps rejected because of possible confusion with the acronym for encapsulated PostScript.
In the first subsection below a brief description of the main operators available for series expansion are given together with some examples of their use.
The most important operator is ps which is used as follows:
ps(EXP:algebraic, VAR:kernel, ABOUT:algebraic):algebraic.
The ps operator returns a Laurent power series object (a tagged domain element) representing the univariate formal Laurent power series expansion of EXP with respect to the dependent variable VAR about the expansion point ABOUT. EXP may itself contain power series objects. If the function has a pole at the expansion point then the correct Laurent series expansion will be produced.
The algebraic expression ABOUT should simplify to an expression which is independent of the dependent variable VAR, otherwise an error will result. If ABOUT is the identifier infinity then the power series expansion about \(\infty \) is obtained in ascending powers of 1/VAR.
Examples
a := ps(sin x, x, 0); ps(sin a, x, 0); ps(cos x/x^2, x, 0); ps(x/(1+x),x,infinity);
As power series objects are domain elements they may be added, subtracted, multiplied and divided in the normal way. For example if A and B are power series objects with the same expansion variable and expansion point:
a+b; a*b; 1/b; a/b;
will produce power series objects representing the sum, product, reciprocal, and quotient of the power series objects A and B respectively.
Differentiation
Similarly, if A is a power series object depending on X then the input df(a, x); will produce the power series expansion of the derivative of A with respect to X.
Integration
The power series expansion of an integral may also be obtained (even if REDUCE cannot evaluate the integral in closed form). An example of this is
ps(int(exp(exp x),x),x,0);
Note that if the integration variable is the same as the expansion variable, the integration package is not called. If on the other hand the two variables are different the integrator is called to integrate each of the coefficients in the power series expansion of the integrand. The constant of integration is zero by default.
Note that the Laurent series domain is not closed under integration with respect to the expansion variable; if the term of degree -1 is non-zero a logarithmic singularity error will occur on integration.
Exponentiation The Laurent series domain is closed under exponentiation by an integer power. Thus, with respect to integer exponentiation, power series are first class objects and for example the following results in automatic expansion of the final result:
a:= ps(cos x,x,0); b:= ps(sin x,x,0); a^2+b^(-2);
However, for more general exponents automatic expansion does not occur. For example given power series a and b defined as above, the following commands are necessary:
ps(a^(1/2),x,0); ps(a^pi,x,0); ps(a^b,x,0);
Note any power of a power series of order zero (that is with a non-zero term of degree zero) can be expanded as a power series (again of order zero) provided only that the power is non-singular at the expansion point. As the third example above shows the exponent may itself be a power series.
However in general the Laurent series domain is not closed under exponentiation. If the result is to be a Laurent series some restrictions on the allowed values of the exponent and order of the original series are necessary. Namely, if the order of the power series is non-zero (\(\sigma \) say) and the exponent is rational with denominator \(q\) say, then \(\sigma q\) must be integral. If the exponent is rational, but \(\sigma q\) is not an integer, a branch point error is generated. For other exponents a logarithmic singularity error is usually generated. For example,
a := ps(1-cos x,x,0); % series has order 2 ps(a^(1/2),x,0); % series has order 1 ps(a^(2/3),x,0); % branch point error ps(a^pi,x,0); % logarithmic singularity error
Power series of user defined functions
New user-defined functions may be expanded provided the user provides a rule or rule list defining the derivative of the function and optionally its value at the expansion point. For example
operator u; let df(u(~x),~x)= exp(e^x); let u(0) = e; ps(u(sin x),x,0);
Of course the rules defined must be such that the function actually has a Taylor series expansion about the specified point.
Restrictions and Known Bugs
Currently automatic expansion of quotients with an integer denominator does not normally occur. One must use:
a:=ps(sin x,x,0); ps(a/5,x,0); or on rational; % or on rounded; a/5;
Currently the following does not produce a power series object (although the result is formally valid):
a := ps(cos x, x, 0); ps(2^a,x,0); % instead use: ps(2^cos x,x,0);
If A is a power series object and X is a variable which evaluates to itself then expressions such as a*x or int(a, x); do not automtically expand to a single power series object (although the result returned is formally valid). Instead expressions such as ps(a*x,x,0) and ps(int(a,x),x,0 should be used.
Currently the handling of essential sigularities is rather erratic; sometimes an Essential Singularity or Logarithmic Singularity error message is output, but often the system fails rather ungracefully.
There is no simple way to write the results of power series calculation to a file and read them back into REDUCE at a later stage.
Taylor Series Expansion
The operator pstaylor may be used as follows:
pstaylor(EXP:algebraic, VAR:kernel, ABOUT:algebraic):algebraic.
which uses the classic Taylor series algorithm for expanding EXP and returning an extendible Taylor series object.
The pstaylor operator may be useful in contexts where the operator ps fails to build a suitable recurrence relation automatically and reports too deep a recursion in ps!:unknown!-crule. A typical example is the expansion of the \(\Gamma \) function about an expansion point which is not a non-positive integer 2.
Note, however, that pstaylor always returns a Taylor series whose order is non-negative. Attempting to use pstaylor to expand a function about a pole will fail with a zero divisor error message.
Also in many cases the use of an automatically generated recurrence relation built by ps is more efficient than using pstaylor, particularly if a large number of terms is required; expansion of tan is a typical example where the number of terms in the nth derivative grows exponentially.
If the command ps or pstaylor is terminated by a semi-colon, a power series object is compiled and then a number of terms of the power series expansion are evaluated and printed.
psexplim Operator
The expansion is carried out as far as the value specified by an internal variable (with a default value of 6). This variable can be accessed via the operator psexplim.
psexplim(UPTO:integer):integer. or psexplim():integer
If psexplim is called with an integer value, the internal variable is updated to the value of UPTO and its previous value is returned. If psexplim is called with no argument the current value is unaltered and that value is returned.
If psexplim is used to increase the expansion limit, sufficient information is stored in the power series object to enable the additional terms to be calculated without recalculating the terms already obtained.
If the command is terminated by a dollar symbol, a power series object is compiled and the first term is calculated, but no output is printed.
When the switch psprintorder is ON the trailing terms of power series beyond psexplim are represented in print by a big-O notation, otherwise, three dots are printed. This switch is ON by default. However, if expression being expanded is a polynomial in the expansion variable and all non-zero terms have been output then the big-O or trailing dots are omitted to indicate that the series is complete.
In this section a number of accessor functions which allow the user to extract information such as the dependent variable, expansion point, a particular term etc. of a power series object.
psdepvar(TPS:power series object):identifier.
The operator psdepvar returns the expansion variable of the power series object TPS. TPS should evaluate to a power series object or an integer, otherwise an error results. If TPS is an integer, the identifier undefined is returned.
psexpansionpt(TPS:power-series-object):algebraic.
The operator psexpansionpt returns the expansion point of the power series object TPS. TPS should evaluate to a power series object or an integer, otherwise an error results. If TPS is an integer, the identifier undefined is returned. If the expansion is about infinity, the identifier infinity is returned.
psfunction(TPS:power-series-object):algebraic.
The operator psfunction returns the function whose expansion gave rise to the power series object TPS. TPS should evaluate to a power series object or an integer, otherwise an error results.
psterm(TPS:power-series-object, NTH:integer):algebraic.
The operator psterm returns the NTH term of the existing power series object TPS. If NTH does not evaluate to an integer or TPS to a power series object an error results. It should be noted that an integer is treated as a power series.
psorder(TPS:power-series-object):integer.
The operator psorder returns the order, that is the degree of the first non-zero term, of the power series object TPS. TPS should evaluate to a power series object or an error results. If TPS is zero, the identifier undefined is returned.
pstruncate(TPS:power-series-object, POWER:integer):algebraic.
This procedure truncates the power series TPS discarding terms of order higher than POWER. The series is extended automatically if the value of POWER is greater than the order of last term calculated to date. For example
a := ps(sin x, x, 0); pstruncate(a, 11);
will output the eleventh order polynomial resulting in truncating the series for \(sin x\) after the term involving \(x^{11}\).
If POWER is less than the order of the series then \(0\) is returned. If POWER does not simplify to an integer or if TPS is not a power series object then a Reduce error result.
In order to functionally invert a power series the operator psreverse is used.
psreverse(TPS:power-series-object) :power-series-object
Four cases arise:
If TPS is not a power series object after evaluation an error results.
Some examples:
ps(sin x,x,0); psreverse(ws); % produces series for asin x about x=0. ps(exp x,x,0); psreverse ws; % produces series for log x about x=1. ps(sin(1/x),x,infinity); psreverse(ws); % series for 1/asin(x) about x=0.
In order to functionally compose two power series the operator pscompose is used.
pscompose(TPS1:power-series-object, TPS2:power-series-object) :power-series-object
The power series TPS1 and TPS2 are functionally composed; that is to say that TPS2 is substituted for the expansion variable in TPS1 and the result expressed as a power series. The dependent variable and expansion point of the result coincide with those of TPS2. The following conditions apply to power series composition:
If these conditions do not hold the series cannot be composed (with the current algorithm terms of the inverted series would involve infinite sums) and a REDUCE error occurs.
Some examples:
a:=ps(exp y,y,0); b:=ps(sin x,x,0); pscompose(a,b); % Produces the power series expansion of exp(sin x) % about x=0. a:=ps(exp z,z,1); b:=ps(cos x,x,0); pscompose(a,b); % Produces the power series expansion of exp(cos x) % about x=0. a:=ps(cos(1/x),x,infinity); b:=ps(1/sin x,x,0); pscompose(a,b); % Produces the power series expansion of cos(sin x) % about x=0.
If an expression is known for the nth term of a power series, an extendible power series object may be constructed by the operator pssum
pssum(J:kernel = LOWLIM:integer, COEFF:algebraic, X:kernel, ABOUT:algebraic, POWER:algebraic) :power-series-object
The formal power series sum for J from LOWLIM to infinity of
COEFF*(X-ABOUT)**POWER
when ABOUT is finite or zero, whereas if ABOUT is infinity
COEFF*(1/X)**POWER
is constructed and returned. This enables power series whose general term is known to be constructed and manipulated using the other procedures of the power series package.
J and X should be distinct simple kernels. The algebraics ABOUT, COEFF and POWER should not depend on the expansion variable X, similarly the algebraic ABOUT should not depend on the summation variable J. The algebraic POWER should be a strictly increasing integer-valued function of J for J in the range LOWLIM to infinity.
Some examples:
pssum(n=0,1,x,0,n*n); % Produces the power series summation for n=0 to % infinity of x**(n*n). pssum(n=1,n,x,0,n); % Produces the power series summation for n=1 to % infinity of n*x**n. pssum(m=1,(-1)**(m-1)/(2m-1),y,1,2m-1); % Produces a power series which is actually the % expansion of atan(y-1) about y=1. pssum(j=1,-1/j,x,infinity,j); % Produces a power series which is actually the % expansion of log(1-1/x) about the point at infinity. pssum(n=0,1,x,0,2n**2+3n) + pssum(n=1,1,x,0,2n**2-3n); % Produces the power series summation for n=-infinity % to +infinity of x**(2n**2+3n).
It should be noted that a formal power series is produced which may not have a non-zero radius of convergence; the second example above illustrates this. Nevertheless these formal series may be added, multiplied, differentiated etc. by the TPS package. Of course, in general the result may also have a zero radius of convergence.
pscopy(TPS:power-series-object):power-series-object
This procedure returns a copy of the power series TPS. The copy has no shared sub-structures in common with the original series. This enables substitutions to be performed on the series without side-effects on previously computed objects. For example:
clear a; b := ps(sin(a*x)), x, 0); b where a => 1;
will result in a being set to 1 in each of the terms of the power series and the resulting expressions being simplified. Owing to the way power series objects are implemented using Lisp vectors, this has the side-effect that the value of b is changed. This may be avoided by copying the series with pscopy before applying the substitution, thus:
b := ps(sin(a*x)), x, 0); pscopy b where a => 1;
pschangevar(TPS:power-series-object, X:kernel):power-series-object
The operator pschangevar changes the dependent variable of the power series object TPS to the variable X. TPS should evaluate to a power series object and X to a kernel, otherwise an error results. Also X should not appear as a parameter in TPS. The power series with the new dependent variable is returned.
psordlim(UPTO:integer):integer or psordlim():integer
An internal variable is set to the value of UPTO (which should evaluate to an integer). The value returned is the previous value of the variable. The default value is 100. If psordlim is called with no argument, the current value is returned.
The significance of this control is that the system attempts to find the order of the power series required, that is the order is the degree of the first non-zero term in the power series. If the order is greater than the value of this variable an error message is given and the computation aborts. This prevents infinite loops in certain cases, for example:
a:=ps(1-(cos x)^2,x,0); b :=ps((sin x)^2,x,0); b-a;
This will also occur in the rather unlikely situation where the expression being expanded is
| Up | Next | Prev | PrevTail | Front |
Hosted by