New Packages for Advanced Mathematics Research in Maple 8

Copyright 2002 Waterloo Maple Inc.

We showcase some of the new Maple 8 packages for advanced mathematics, including

• Calculus of variations
• Symbolic-numeric algorithms for polynomials
• Matrix polynomial algebra
• New conversion and evaluation routines for special functions

Calculus of Variations with the VariationalCalculus  package

Introduction

The new VariationalCalculus  package provides routines for solving problems in the calculus of variations, which studies nature's most "efficient" curves and surfaces.  Classical problems from the calculus of variations include:

• Find the shortest path between two points on a 3-D surface, such as a cone (a geodesic problem).
• Shape a ramp between two heights such that a ball rolling down it reaches the bottom in minimum time (the Brachistochrone problem).
• Find the shape of a soap film having minimum surface area spanning a given wire frame.

Such problems can often be solved with the Euler-Lagrange equation, which generalizes the Lagrange Multiplier Theorem for minimizing functions of real variables subject to constraints.  The Euler-Lagrange equation is easy to write down in general but notoriously difficult to write down and solve for most practical problems.  The VariationalCalculus  package automates the construction and analysis of the Euler-Lagrange equation.

>

Example: The Brachistochrone Problem

The Brachistochrone problem can be stated as follows: Given two endpoints in the plane, find the curve y(x) between them such that a ball of unit mass rolls along the curve under the influence of gravity in minimum time.

>

restart;

>	with(VariationalCalculus);

First we write down the falling time over an infintesimal distance dx in terms of y(x)  and yInit , assuming the gravitational constant is 1.  This is found in standard textbooks on classical mechanics.

>	fallTime := sqrt( (1+diff(y(x),x)^2)/(2*(yInit-y(x))) );

We then use the EulerLagrange  function to compute the Euler-Lagrange equations for this functional in terms of y(x)  and its derivatives.  This returns a set of ODEs.

>	EL := EulerLagrange( fallTime, x, y(x) ) ;

Extract the first ODE from the list

>	ODE1 := remove(has, EL, diff(y(x),x,x))[1];

Compute y(x)  as a parametric curve (x(s), y(s)) using the dsolve  command with the parametric  option.

>	sol := dsolve(ODE1, parametric, y(x));

>	sol := subs( {_T=s}, sol);

>	assign(sol);

Specify the endpoints of the curve: (0,1)  and (1, 1/2)

>	xInit := 0; yInit := 1; xEnd := 1; yEnd := 1/2;

The parameter s  is not equivalent to time t , so we must solve for the value of s  at which

>	solve(y(s)=yEnd, s);

>	sEnd := %[1];

Use the endpoints to solve for the unknown constants K[1]  and _C1

>	_C1 := solve(subs({s=sEnd}, x(s)) = xEnd, _C1);

The finite-length curve (x(s), y(s))  is parametrized over , so

>	limit(x(s), s=-infinity);

>	K[1] := fsolve( % = xInit, K[1]);

The path along which a ball would roll in minimum time under the influence of gravity.  The result is surprising because the curve actually dips below the right end point.

>	plot([x(s),y(s), s=-10..10], xInit..xEnd, 0..yInit);

SNAP: Symbolic-Numeric Algorithms for Polynomials

Introduction

Many polynomials used in practice have inexactly known coefficients. You may have a polynomial that describes your data, fitted using least-squares. You may have a formula with parameters for each coefficient, but the values of the parameters may be known only through experimental measurements. Such polynomials arise in many applied problems, from direct mathematical models through normal forms and bifurcation theory to geometric problems described in floating-point arithmetic for Computer-Aided Design.

The idea that we don't know the coefficients precisely moves us from the domain of algebra to the domain of analysis. In particular, this allows us to bring to bear modern tools of numerical analysis on problems of computational algebra and geometry.

Experience with numerical linear algebra suggests that it will be useful to replace classical discrete algorithms for polynomial problems with continuous, iterative algorithms for approximate polynomial problems.

The new SNAP  package in Maple 8 provides tools for the algebraic manipulation of numeric polynomials. It includes support for quotient, remainder and numeric gcd, numerical computations that give the last numerically stable euclidean reduction, an approximation to the distance of a closest common root of two relatively prime numeric polynomials and similar operations.

Examples using the SNAP Package

>

restart;

>	with(SNAP);

The AreCoprime  function determines whether two numeric polynomials are relatively prime up to a given error bound . For example,

>	a := 0.1*z^2+1.5*z-0.2;

>	b := 0.2*z^3+0.15;

>	AreCoprime(a,b,z,0.5);

>	AreCoprime(a,b,z,0.1);

The EpsilonGCD(a,b,z)  function returns a univariate numeric polynomial g  together with a positive float epsilon such that g  is an epsilon-GCD for the input polynomials ( a , b ).

>	a := -.2313432836*z^4+.3500000000e-2*z^3-.1753694030*z^2-.3397276119*z-.3395522388e-3;

>	b := -.2313432836*z^3+.3731343284e-2*z^2-.1753731343*z-.3395522388;

>	EpsilonGCD(a,b,z);

The distance to common divisors is near 0, meaning that an epsilon-GCD exists.

>	DistanceToCommonDivisors(a, b, z);

>	AreCoprime(a, b, z, .7609811447e-11);

>	c := z^2+3.1*z-2;

>	d := 2*z^3+1.5;

>	EpsilonGCD(c, d, z);

The distance to common divisors is relatively large, meaning an epsilon-GCD could not be found.

>	DistanceToCommonDivisors(c, d, z);

>	AreCoprime(c, d, z, .8);

The EuclideanReduction(a,b,z)  function computes in a stable manner a polynomial pair that is equivalent to the input polynomials ( a , b ) but  has smaller degrees.

>	a := z^6-12.4*z^5+50.18112+62.53*z^4-163.542*z^3+232.9776*z^2-170.69184*z;

>	b := z^5-17.6*z^4+118.26*z^3-372.992*z^2-274.09272+538.3333*z;

>	EuclideanReduction(a,b,z);

Matrix Polynomial Algebra

Introduction

Matrices of polynomials have many algebraic properties that are similar to ordinary polynomials. They may be added, multiplied, have degrees and may be placed in normal form. They  may even have right and left factors and right or left greatest common factors.   The MatrixPolynomialAlgebra  package provides tools for the algebraic manipulation of matrix polynomials. In particular, there are tools that perform basic polynomial operations , such as determining degrees and coefficients; perform algebraic operations , such as matrix polynomial left and right division, and determining greatest common divisors and least common multiples; convert to special forms of matrix polynomials , such as normal forms and reduced forms; and determine bases for the kernel of a matrix polynomial .

Techniques based on matrix polynomials are  used in  many applications, for example in linear system theory, control theory  and signal processing. For example, multi-input/multi-output (MIMO) linear systems of the form  can be studied (after the taking of Laplace transforms) through the transfer function .  Similarly, when designing feedback control systems, matrix polynomial equations of the form  may arise, which yield the controller transfer matrix .  In both applications given above, one needs to be able to work with polynomial matrix fractions. But polynomial matrix fractions  do not  have a unique representation, and so a  fundamental question is how to get the "best" matrix fraction  representation. For example, "best" means that  is  irreducible  or that the denominator  is in canonical form . Hence a need for coprimeness tests and matrix polynomial GCDs and for algorithms computing various reduced and canonical forms such as the Hermite, Popov and Smith forms.

Examples using the Matrix Polynomial Algebra package

>

restart;

>	with( MatrixPolynomialAlgebra );

>	A := <<3+x,4,x^2-1>|<1,x,4>|<-4*x,2,-x^3>>;

The Coeff  command extracts the coefficients of the powers of a variable in a polynomial matrix.

>	seq( Coeff(A, x, i), i=0..3);

The LeftDivision(A, B, x)  function divides B into A on the left, producing  a  quotient Q  and a remainder R  such that A = B . Q + R .  The RightDivision () function works analogously. For example, if

>	A := Matrix(2,2,[[-9*z^2 - 3*z + 1, 12*z^2 + 10*z],                  [-3*z^3+2*z^2-z,4*z^3+2*z - 2*z^2]]); B := Matrix(2,2,[[-3*z^3 + 6*z^2 + 5*z + 1, -12*z^2-13*z],                  [z^4 + z^3 + z^2, -4*z^3 - 3*z + 3*z^2]]);

then a quotient Q  and remainder   R  can be determined via

>	Q, R := op( LeftDivision(A, B, z) );

 The result can then be verified via

>	map(expand, A - (B . Q + R));

The MatrixGCRD(A,B,z)  function computes the greatest common right divisor of A  and B  along with the associated matrix diophantine equation

>	C := MatrixGCRD( A, B, z, 'U' );

with the information about the diophantine equation given in the list U

>	map( expand, U[1] . A + U[2] . B );

The HermiteForm[row](A)  function computes the row reduced echelon form of an m x n  rectangular matrix of univariate polynomials in x  over the field of rational numbers Q , or rational expressions over Q

>	A := <<3+x,4,x^2-1>|<1,x,4>|<-4,2,-1>>;

>	H := HermiteForm(A, x);

>

Modular Linear Algebra

The LinearAlgebra:-Modular  sub-package is a highly efficient suite of tools for performing dense linear algebra computations in , the integers modulo m over the positive range.

As a programmer package, the key focus is on efficient computation mod m, so for all routines that accept modular data in Matrix/Vector form, no checking of the input data (to verify that all data is in the range 0..m-1) is performed.

>

restart;

>	with(LinearAlgebra:-Modular);

>	Mat := Mod(13, [[1, x, x^2], [x-2, 3, x/(x-1)]], x=4, float[8]);

>	RowReduce(13, Mat, 2 ,3, 3, 'det', 0, 0, 0, 0, true): Mat;

>	dtype := integer[kernelopts(wordsize)/8]: V1 := Mod(5, [5, 2, 3, 4, 1], dtype);

>	V2 := Mod(5,[3,1,1,9,1],'transpose',dtype);

>	Multiply(5,V1,V2);

>	Multiply(13,V2,V1);

>

Conversions between Special Functions

A flexible conversion facility in the mathematical language is as important as a dictionary in spoken languages. Maple 8 introduces such a tool in a computer algebra system for the first time, implemented as a net of conversion routines permitting the expression of any mathematical function in terms of another one, whenever that is possible, as a finite sum of terms. When the parameters of the functions being converted depend on symbols in a rational manner, any assumptions on these symbols made by the user are taken into account at the time of performing the conversions (see assuming ).

• Example:  

>

restart;

>	ee := exp(1/2*z)* WhittakerM(0,1/2,z);

>	convert( %, MeijerG );

>	convert( %, exp );

>	convert( %, hypergeom, include=exp );

Converting the original expression ee  directly to hypergeom  we arrive at an identity for hypergeometric functions

>	% = convert(ee, hypergeom);

>	simplify( (lhs-rhs)(%) );

The main idea behind the Maple 8 conversion network is to split the set of mathematical functions into three main subsets, according to their hypergeometric representation as 2F1, 1F1 and 0F1 hypergeometric functions. Then, generally speaking, conversions - when possible - can be performed among functions of the same class. The Maple 8 routines allow you to convert to any of 58 mathematical functions, including all the elementary transcendental ones and most of the special functions of mathematical physics:

arccos ,        arccosh ,     arccot ,      arccoth ,      arccsc ,     
arccsch ,       arcsec ,      arcsech ,     arcsin ,       arcsinh ,    
arctan ,        arctanh ,     cos ,         cosh ,         cot ,        
coth ,          csc ,         csch ,        exp ,          ln ,         
sec ,           sech ,        sin ,         sinh ,         tan ,        
tanh ,          GAMMA ,       dilog ,       polylog ,      AiryAi ,     
AiryBi ,        HankelH1 ,    HankelH2 ,    BesselI ,      BesselJ ,    
BesselK ,       BesselY ,     MeijerG ,     hypergeom ,    Ei ,         
erf ,           KummerM ,     KummerU ,     WhittakerM ,   WhittakerW ,
CylinderU ,     CylinderV ,   CylinderD ,   LaguerreL ,    HermiteH ,   
GegenbauerC ,   LegendreP ,   LegendreQ ,   ChebyshevT ,   ChebyshevU ,
EllipticE ,     EllipticK ,   JacobiP                                

• Example:  Conversions between 2F1 functions

>	JacobiP(1/2,0,-1,z);

>	convert(%,EllipticE);

>	convert(%,LegendreP);

>	convert(%,GegenbauerC);

>	convert(%,hypergeom);

• Example:  Conversions between 1F1 functions

>

erf(z);

>	convert(%,HermiteH);

>	convert(%,KummerU);

>	normal(convert(%,WhittakerW));

Taking into account the possible split of the set of mathematical functions into hypergeometric classes as well as their classification found in typical handbooks like Abramowitz and Stegun, the Maple 8 routines also permits converting to a function class ; that is, attempting to rewrite a given expression by using any of the functions of a specified class.

• Example:  

>	1/2*LegendreP(-1/2,-1/2,-1-2*z^2+4*z)*(-2*z^2+4*z)^(1/4)/(-2-2*z^2+4*z)^(1/4)*(z-1);

>	convert(%,arctrig);

The function classes understood by the Maple 8 mathematical function conversion network are:

trig ,               trigh ,        arctrig ,                
arctrigh ,           elementary ,   GAMMA_related ,          
Kelvin ,             Airy ,         Hankel ,                 
Bessel_related ,     0F1 ,          Ei_related ,             
erf_related ,        Kummer ,       Whittaker ,              
Cylinder ,           1F1 ,          Orthogonal_polynomials ,
Elliptic_related ,   Legendre ,     Chebyshev                 

The new conversion routines also introduce the concept of "rule conversion", that is, a conversion where the input and output are expressed using the same function but applying identity rules to it. For example, the HermiteH  function

>	H1 := HermiteH(a,z);

can have its (first) parameter raised or lowered by one. Below, "raise a" stands for "raise the first parameter":

>	H2 := convert(H1, HermiteH, "raise a");

For functions depending on more parameters, for example, LegendreP  or JacobiP , you can raise the second or third parameters using "raise b" or "raise c". The rules known by the Maple 8 conversion routines are:

"raise a" ,       "lower a" ,      "normalize a" ,
"raise b" ,       "lower b" ,      "normalize b" ,
"raise c" ,       "lower c" ,      "normalize c" ,
"mix a and b" ,   "1F1 to 0F1" ,   "0F1 to 1F1"     

The mathematical identities which can be derived using the new conversion routines can be checked visually in Maple 8 using the new plotcompare  routine, which draws four 3-D plots, comparing the real and imaginary parts of two expressions when evaluated over a complex variable z = x + I y .

All the Maple 8 conversion routines listed in the tables above understand a uniform set of optional arguments for restricting the conversion process in varied manners, for example, performing (or excluding) the conversion process only for some functions (see convert,to_special_function ), therefore providing the necessary flexibility for optimal symbolic manipulation.