Advanced Math Improvements in Maple 2023
Maple 2023 includes many improvements to the math engine. In addition to those described on this page, Maple 2023 includes a new package for Quantifier Elimination and updates to the Graph Theory package.
Better answers for indefinite integrals of algebraic functions
Enhanced 'method' option for indefinite integrals
Better assumption support of MeijerG integrals
Support in solve for element-wise equations of rtables
Polyhedral sets
New support for Z-polyhedral sets
Improvements to parametric limits on branch cuts
New certified univariate complex solver
New method option for fsolve, for a univariate polynomial
New method option for fsolve, for a univariate non-polynomial equation
New support for lazy multivariate Puiseux series
The int command now uses Blake's heuristic to look for substitutions to find elementary forms for certain algebraic integrands. Blake's method works by considering algebraic function integrands of the form pxrxnmqx and then searching for substitutions (by solving systems of linear equations) that can put the integral into a simpler form, e.g. ∫pxrxnmqxⅆx=∫vuau+bmwuⅆu, for which compact elementary answers can often be easily computed.
Blake's heuristic has been available since Maple 2022 via the option method=pseudoelliptic but it is now tied directly into the default methods for computing integrals of algebraic functions.
Examples
The following examples produced a very complicated expression in earlier versions of Maple, but now produce a simple, elementary answer.
This is an example where in Maple 2022 the default answer comes from the elliptic integration method, and the answer is in terms of elliptic functions and is quite long.
intx3−2−x3−1x3+1x3+x2+1,x,method='elliptic'
2I3Ix−12−I323x+132+I32−Ix−12+I323EllipticF3Ix−12−I3233,I332+I323−x3−1+I2∑_α=RootOf_Z3+_Z2+1_αI2−1+2x−I3x+13+I3−I2−1+2x+I3−_α2+3+I3−_α2−2_α−1EllipticPi3Ix−12−I3233,3_α22+3_α−I_α232+32+3I32,I332+I32−x3−1
In Maple 2023 the answer is elementary and compact.
∫x3−2−x3−1x3+1x3+x2+1ⅆx
2arctanh−x3−1x
In Maple 2022 the following integral produced an answer in terms of elliptic functions and RootOfs of length 790119, but now the answer is elementary and very compact.
∫x4+x3+1−x4−2x3+1x4+2x3−1x4+2x3+x2−1ⅆx
arctanh−x4−2x3+1x
Finally, here is a fairly long integral with a compact new answer in Maple 2023. method=elliptic shows the answer previously returned in Maple 2022.
∫x3+2112x6−4x5+3x4−224x3+4x2+1122x3+x2−212x416x6−8x5−x4−32x3+8x2+16ⅆx
−95−2−12495x36+22arctanh22x3+x2−2x6−22+x32+12495arctanh22x3+x2−2x6+22−1122x3+x2−2x3+5328x2−16+222856−2286−226+22x3
intx3+2112x6−4x5+3x4−224x3+4x2+1122x3+x2−212x416x6−8x5−x4−32x3+8x2+16,x,method='elliptic'
−142x3+x2−23x3+532x3+x2−26x+142x3+x2−23+31I3107+6318136−16107+631813−Ix+107+63181312+112107+631813+16+I3107+6318136−16107+63181323107+6318136−16107+631813x−107+6318136−16107+631813+16−107+6318134−14107+631813−I3107+6318136−16107+6318132Ix+107+63181312+112107+631813+16−I3107+6318136−16107+63181323107+6318136−16107+631813EllipticF3−Ix+107+63181312+112107+631813+16+I3107+6318136−16107+63181323107+6318136−16107+6318133,−I3107+6318136−16107+631813−107+6318134−14107+631813−I3107+6318136−16107+631813262x3+x2−2+I2∑_α=RootOf16_Z6−8_Z5−_Z4−32_Z3+8_Z2+16_α1520_α4−760_α3−95_α2−1520_α+876107+631813−1107+631813−I22+12x+107+631813+1107+631813+I3107+631813−1107+631813107+631813−1107+6318131+6x−107+631813−1107+631813−3107+631813−3107+631813−I3107+631813−1107+631813I22+12x+107+631813+1107+631813−I3107+631813−1107+631813107+631813−1107+631813−320+61063107+6310613−1861063107+6310623+3316107+6310623−76107+6310613+64_α5+352_α4−68_α3−216_α2−528I106_α3107+6310623+3137I3_α3107+6310623−5730I3107+6310623_α2+66I3107+6310623_α+963I106107+6310623_α2−99I106107+6310613_α2+642I3107+6310613_α2−592I3_α5107+6310613+96I106_α5107+6310613−384I106_α5107+6310623+2288I3_α5107+6310623−12I106107+6310623_α−96I106_α107+6310613+570I3_α107+6310613+8I3_α4107+6310613+144I106_α4107+6310623−856I3_α4107+6310623+149I3_α3107+6310613−30I106_α3107+6310613−10_α31063107+6310613+176_α31063107+6310623+32_α51063107+6310613+128_α51063107+6310623−48_α43107+6310623106−321_α2107+63106231063−32106_α3107+6310613+4_α1063107+6310623−33_α21063107+6310613−76I3107+6310613−3316I3107+6310623+558I106107+6310623+18I106107+6310613+570_α107+6310613−66107+6310623_α+5730107+6310623_α2+642107+6310613_α2+149_α3107+6310613−3137_α3107+6310623+856_α4107+6310623+8_α4107+6310613−592_α5107+6310613−2288_α5107+6310623−360_αEllipticPi3−Ix+107+63181312+112107+631813+16+I3107+6318136−16107+63181323107+6318136−16107+6318133,−311063107+631061356+1063107+631062356−19107+631062384+829107+631061384+74_α521−_α421−149_α3168+I10628−107_α228+1942−107I_α2107+6310613106112+991I_α23107+6310613168−7I_α53107+63106133−I107+6310613_α41067−_α41063107+63106137+8I_α5107+631061310621−26I_α5107+631062310621+155I_α53107+631062321+137I_α107+631062310621+I_α107+631061310684+7I_α3107+631061324+I_α43107+63106132−44I_α4107+63106231067+523I_α43107+631062314−341I_α33107+6310613112+11I_α3107+631061310621+209I_α3107+6310623106168−825I_α33107+6310623112−3959I_α23107+6310623168−6515I_α3107+6310623168+443I_α2107+6310623106112+11_α31063107+631061321−5_α31063107+6310623168+8_α51063107+631061321+2_α51063107+631062321−11_α2107+63106231063112+106_α3107+631061384−2_α1063107+631062321−107_α21063107+6310613112−11_α107+631061356+95107+6310623_α56+107107+6310623_α256+955107+6310613_α256−3137_α3107+6310613336+149_α3107+6310623336+_α4107+631062342+107_α4107+631061342−143_α5107+631061321−37_α5107+631062321−5I_α310684−4I_α10621+4I_α510621−11I_α210656+101I3107+631061328−949I107+6310623328+319I107+631062310656−31I107+631061310656−95_α28,−I3107+6318136−16107+631813−107+6318134−14107+631813−I3107+6318136−16107+63181322x3+x2−264512
The int command now exposes more internal routines for indefinite integration through the method option. It is now possible to directly call the integration by parts routine Parts and a parallel version of the Risch algorithm ParallelRisch.
int(ln(exp(x)),x,method=parallelrisch);
−x22+lnⅇxx
int(4*x^2*arcsinh(x),x,method=parts);
4x3arcsinhx3−4x2x2+19+8x2+19
The definite integration method that works via MeijerG convolutions now does a better job of checking conditions on parameters so that they are only applied under proper assumptions. As well, the methods will prompt via warning when the method could have produced an answer using additional assumptions.
int( x*exp(-a*x^3 - b*x), x=0..infinity );
Warning, not enough assumptions for a closed form integral: a MeijerG transformation may be possible with additional assumptions: a <> 0, 0 <= Re(a), b <> 0, 0 <= Re(b)
intxexp−ax3−bx,x=0..∞
assume(a::real, b::real);
Warning, not enough assumptions for a closed form integral: a MeijerG transformation may be possible with additional assumptions: 0 < a, 0 < b
intxexp−a~x3−b~x,x=0..∞
assume(c>0, d>0);
int( x*exp(-c*x^3 - d*x), x=0..infinity );
d~d~32c~13−1+LommelS21,23,2d~3239c~3c~56
restart:
The solve command can now directly accept element-wise relations between vectors and matrices. Inputs of the form M, M=c, or M=N are transformed into a list of scalar equations (where M and N are rtables of the same dimensions and c is any scalar, and instead of = any binary relation of type,relation can be used).
solve( (<<x+1|y+1>,<z+2|w+2>>) = 1 );
w=−1,x=0,y=0,z=−1
With no relations, a right-hand side of 0 is assumed
solve( (<<x+1|y+1>,<z+2|w+2>>) );
w=−2,x=−1,y=−1,z=−2
solve( (<<x+1|y+1>,<z+2|w+2>>) = (<<1|1>,<2|2>>) );
w=0,x=0,y=0,z=0
solve( (<<x+1|y+1>,<z+2|w+2>>) < (<<1|1>,<2|2>>) );
w<0,x<0,y<0,z<0
While it is now possible to solve linear equations this way, it is advisable to use LinearAlgebra:-LinearSolve instead.
n := 200:
M := LinearAlgebra:-RandomMatrix(n):
b := LinearAlgebra:-RandomVector(n):
v := Vector( [ seq(cat(`v__`,i), i=1..n)]):
CodeTools:-Usage( eval(v, solve( M.v = b )) ):
memory used=110.74MiB, alloc change=40.17MiB, cpu time=2.10s, real time=1.90s, gc time=436.46ms
CodeTools:-Usage( LinearAlgebra:-LinearSolve(M, b) ):
memory used=54.78MiB, alloc change=1.01MiB, cpu time=1.28s, real time=1.16s, gc time=292.64ms
The IntegerHull command in the PolyhedralSets package has been extended to support higher dimensional polyhedral sets, both bounded and unbounded ones.
The following example is an unbounded three-dimensional set.
with(PolyhedralSets):
ineqs := [-x[1] - 132/205*x[2] - 62/205*x[3] <= -1358/123, -x[1] + x[2] + x[3] <= 1405/17, x[1] - 6/59*x[2] + 83/177*x[3] <= 3500/59];
ineqs≔−x1−132x2205−62x3205≤−1358123,−x1+x2+x3≤140517,x1−6x259+83x3177≤350059
poly := PolyhedralSet(ineqs, [x[1], x[2], x[3]]);
PolyhedralSets:-PolyhedralSet1358123−132205x2−62205x3≤x1,−140517+x2+x3≤x1,x1≤659x2−83177x3+350059,x1,x2,x3
IsBounded(poly);
false
IntegerHull(poly);
−20,36,26,−4,−25,103,−2,−30,107,−1,−36,117,0,−38,118,0,−36,118,1,−37,112,1,−34,117,2,−39,113,2,−38,114,10,−43,95,26,−51,60,399,−238,−776,403,−240,−785,453,−265,−897,1544,−811,−3342,101,260,−159,12072,−6041,−27054,−70,267,−337
A bounded four-dimensional example.
vertices := [[10, 10, 10, 10/3], [-140/8, -220/12, -10, -10/3], [60/8, 20, -100/12, -70/3], [-10/4, -100/12, 70/2, 35/3], [0, 0, 0, 50/3]];
vertices≔10,10,10,103,−352,−553,−10,−103,152,20,−253,−703,−52,−253,35,353,0,0,0,503
vars := [x[1], x[2], x[3], x[4]]:
poly := PolyhedralSet(vertices, [], vars);
PolyhedralSets:-PolyhedralSet−77753123+503694x2+85694x3+3112082x4≤x1,−67751117+27152234x2+6032234x3+8132234x4≤x1,x1≤−543+7586x2+59430x3+3430x4,x1≤−543+7586x2−33182x3+13413182x4,x1≤−145148+123185x2+153370x3+871480x4,x1,x2,x3,x4
true
−15,−16,−6,−2,−15,−15,−9,−4,−14,−15,−4,−1,−13,−13,−8,−1,−13,−12,−9,−5,−12,−13,−4,1,−12,−13,−4,2,−12,−12,−3,−3,−11,−12,−3,−1,−11,−11,−6,−3,−11,−11,−1,−3,−10,−8,−8,−5,−9,−6,−8,−8,−7,−7,−4,7,−7,−6,−5,3,−7,−3,−8,−10,−7,−3,−7,−10,−6,−4,−5,0,−5,−9,23,8,−5,−4,−4,5,−5,−4,−3,−2,−4,−5,3,10,−4,−5,3,11,−4,−3,−2,−3,−4,0,−6,−9,−3,−8,30,10,−3,−7,23,10,−3,−7,23,11,−3,−6,24,6,−3,3,−8,−12,−3,3,−7,−13,−2,−7,31,11,−2,−6,24,8,−2,−6,24,12,−2,−6,26,11,−2,−5,25,7,−2,−5,26,6,−2,−2,−1,13,−2,−1,−1,1,−2,3,−4,−8,−1,−5,25,12,−1,−5,29,9,−1,−4,19,13,−1,−4,27,7,−1,−3,12,13,−1,−3,12,14,−1,−1,1,14,−1,−1,2,0,−1,1,11,−4,0,−2,22,7,0,−1,8,14,0,−1,20,5,0,0,1,15,0,0,2,15,0,0,13,6,0,3,−3,4,0,8,−8,−16,1,−3,27,9,1,−2,23,9,1,0,8,3,1,0,22,3,1,2,0,12,1,2,5,9,1,2,9,6,1,3,−1,9,1,5,5,−5,1,8,−3,−13,2,1,9,3,2,1,21,3,2,2,2,14,2,2,4,1,2,5,7,0,2,5,13,−5,2,10,−7,−13,2,11,−8,−17,3,7,−2,3,3,7,10,−6,3,7,11,−7,3,8,−3,0,3,8,6,−7,3,11,−2,−15,4,2,17,6,4,8,8,−5,4,9,−2,−3,4,10,−1,−5,4,10,3,−8,4,11,−3,−13,4,11,−2,−13,4,12,−5,−15,4,12,0,−15,4,14,−7,−19,5,3,18,6,5,3,18,7,5,3,19,6,5,5,6,3,5,13,−5,−10,5,13,1,−14,5,14,−6,−13,5,15,−7,−18,5,15,−7,−17,6,9,3,−4,6,15,−4,−14,6,16,−6,−16,6,17,−7,−20,7,6,14,5,7,12,1,−7,7,13,0,−10,7,13,2,−9,7,13,3,−10,7,15,−3,−12,7,16,−4,−16,7,17,−4,−18,7,18,−6,−20,8,8,8,6,8,8,9,3,8,8,12,3,8,12,3,−4,9,10,8,1,9,10,8,2,9,10,9,1,9,13,4,−5,
The ZPolyhedralSets subpackage is new in Maple 2023. It is part of the PolyhedralSets package and deals with Z-polyhedral sets. A Z-polyhedral set is the intersection of a polyhedral set with an integer lattice.
with(PolyhedralSets): with(ZPolyhedralSets):
Here is a three-dimensional polyhedral set.
ineqs := [0 <= -16 + 2*y + z, 0 <= -72 + 4*x + 4*y + 3*z, 0 <= 2*y - z, 0 <= -24 + 4*x + 4*y - 3*z, 0 <= -4*x + 4*y + 3*z, 0 <= 48 - 4*x + 4*y - 3*z, 0 <= 48 - 4*x - 4*y + 3*z, 0 <= 8 - 2*y + z, 0 <= -24 + 4*x - 4*y + 3*z, 0 <= 24 - 2*y - z, 0 <= 24 + 4*x - 4*y - 3*z, 0 <= 96 - 4*x - 4*y - 3*z]:
ps := PolyhedralSet(ineqs):
Plot(ps);
Now we define the Z-polyhedral set of all integer lattice points in ps.
psz := ZPolyhedralSet(ineqs, [x, y, z]):
PlotIntegerPoints3d(psz);
The option parametric to the limit command has been extended to handle more cases where the expansion point may lie on the branch cut of a mathematical function for some real parameter values. Examples:
unassign('a'):
limit(ln(a+I*x), x=0, 'right', 'parametric');
ln−a+Iπa<0−∞a=0lna0<alimx→0+lna+Ixotherwise
limit(arccoth(a+I*x), x=0, 'right', 'parametric');
−arccoth−aa<−1−∞a=−1arctanha−Iπ2a<1∞a=1arccotha1<alimx→0+arccotha+Ixotherwise
limit(Li(a+I*x), x=0, 'right', 'parametric');
Eiln−a+Iπa<0Iπa=0Iπ+Eilnaa<1−∞a=1Eilna1<alimx→0+Lia+Ixotherwise
A new, faster method for isolating the complex roots of a univariate polynomial with complex(numeric) coefficients has been added to the RootFinding:-Isolate command. It is accessible through the option complex or method=HR and supports most of the existing options of RootFinding:-Isolate, as well as a new option domain for limiting the search space. See the corresponding help page for more details. Examples:
with(RootFinding):
Isolate(x^6-1, x, 'complex');
x=−1.,x=1.,x=0.5000000000+0.8660254038I,x=0.5000000000−0.8660254038I,x=−0.5000000000+0.8660254038I,x=−0.5000000000−0.8660254038I
Isolate(x^2-I, x, 'complex', 'output=interval');
x=8548396450010082566043051208925819614629174706176+854839645001008256604305I1208925819614629174706176,8548396450010102904031831208925819614629174706176+854839645001010290403183I1208925819614629174706176,x=−8548396450010102904031831208925819614629174706176−854839645001010290403183I1208925819614629174706176,−8548396450010082566043051208925819614629174706176−854839645001008256604305I1208925819614629174706176
Isolate(x^2-I, x, 'complex', 'digits'=100);
x=0.7071067811865475244008443621048490392848359376884740365883398689953662392310535194251937671638207864+0.7071067811865475244008443621048490392848359376884740365883398689953662392310535194251937671638207864I,x=−0.7071067811865475244008443621048490392848359376884740365883398689953662392310535194251937671638207864−0.7071067811865475244008443621048490392848359376884740365883398689953662392310535194251937671638207864I
Find 10 roots in the first quadrant.
Isolate(x^100+x^3-1, x, 'method=HR', 'domain'=[0,infinity+I*infinity], 'maxroots'=10);
x=0.9744076649,x=0.9781052241+0.1737357669I,x=0.9531760998+0.2980270829I,x=0.9680195472+0.2362266682I,x=0.9104543811+0.4181438477I,x=0.9827791401+0.1111417351I,x=0.9804852800+0.05036970385I,x=0.1198693017+0.9978851093I,x=0.4848690063+0.8826473106I,x=0.2461381235+0.9755862161I
The fsolve command has been updated with a new method option. When solving for roots of a univariate polynomial this option allows specification of the new solver within the RootFinding:-Isolate command, or the original modified Laguerre method of the NAG Library. The value for the option can be specified here as method=hefroots or method=NAG. The current default is method=NAG, unless fsolve determines that the problem is ill-conditioned.
[ fsolve(x^4-5.591265174*x^3+2.936595076*x^2+17.11912543*x-10.52786025, x=0..4, method=hefroots) ];
0.6180339885,2.718281830,3.872983344
When solving for roots of an equation in one variable that is not a polynomial, the new method option to fsolve lets you specify the method used. The value for the option can be specified here as method=subdivide or method=NextZero. This option must be accompanied by the maxsols option. The subdivide method repeatedly searches on new subintervals between roots already found, using a custom mix of approaches on each. The NextZero method utilizes the RootFinding:-NextZero command to repeatedly search past the last root found. The default is method=subdivide, which may be more suitable for equations that are not continuous or have non-isolated roots. For some problems the NextZero method may be significantly faster.
[ fsolve(BesselJ(0,x^3)-1/100*x^2-floor(x)+2, x=0..4, maxsols=10, method=subdivide) ];
2.040911370,2.290794873,2.446318006,2.642712321,2.748619065,2.920713827,2.995225008
[ fsolve(BesselJ(0,x^3)-1/100*x^2, x=0..4, maxsols=10, method=':-NextZero') ];
1.333416382,1.777268850,2.040911369,2.290794872,2.446318005,2.642712321,2.748619064,2.920713827,2.995225007,3.155131991
The MultivariatePowerSeries package was introduced in Maple 2021. It deals with multivariate power series in a lazy fashion, meaning that extra terms for a given result can be computed quickly.
In Maple 2023, this package can now deal with multivariate Puiseux series and univariate polynomials over multivariate Puiseux series.
with(MultivariatePowerSeries):
a := PowerSeries(exp(s + t)*ln(1+t));
ObjectPowerSeriesObject,Array(0..5, {(1) = 0, (2) = t, (3) = `+`(`*`(`/`(1, 2), `*`(`^`(t, 2))), `*`(s, `*`(t))), (4) = `+`(`*`(`/`(1, 3), `*`(`^`(t, 3))), `*`(`/`(1, 2), `*`(s, `*`(`^`(t, 2)))), `*`(`/`(1, 2), `*`(`^`(s, 2), `*`(t)))), (5) = `+`(`*`(`/`(1, 3), `*`(s, `*`(`^`(t, 3)))), `*`(`/`(1, 4), `*`(`^`(s, 2), `*`(`^`(t, 2)))), `*`(`/`(1, 6), `*`(`^`(s, 3), `*`(t))))}),5,mul_gen,s,t,ancestors
b := PuiseuxSeries(a, [s = sqrt(x), t = y^(1/3)*x]);
ObjectPuiseuxSeriesObject,modulePowerSeriesObjectoptionobject;localhpoly::Array,deg::nonnegint,gen::procedure,vars::setname,ancestors::record,algexpr::algebraic,dstyle::list:-identical:-maxterms,:-precision=:-nonnegint,:-identical:-∞;end module,x,y,s,t,12,0,1,13,
Truncate(1/b, 3, mode=absolute);
−12+1xy13−1xy13+12y13−x6y13+x2+x24y13−x4−x32120y13−y13x12+x3212+x2720y13+x32y1312−x248−x525040y13−x2y1324+x52240+5y23x224+x52y1372−x31440
The package can also compute a Puiseux factorization of a univariate polynomial over univariate Puiseux series. We set this up by defining a univariate Puiseux series.
c := PowerSeries(exp(t) * ln(1+t));
ObjectPowerSeriesObject,Array(0..5, {(1) = 0, (2) = t, (3) = `+`(`*`(`/`(1, 2), `*`(`^`(t, 2)))), (4) = `+`(`*`(`/`(1, 3), `*`(`^`(t, 3)))), (5) = 0}),5,algebraic_expression_gen_series,t,ancestors
d := PuiseuxSeries(c, [t = y^(1/3)]);
ObjectPuiseuxSeriesObject,modulePowerSeriesObjectoptionobject;localhpoly::Array,deg::nonnegint,gen::procedure,vars::setname,ancestors::record,algexpr::algebraic,dstyle::list:-identical:-maxterms,:-precision=:-nonnegint,:-identical:-∞;end module,y,t,13,
We now define e as dz2+1dz+1 and compute its Puiseux factorization.
e := UnivariatePolynomialOverPuiseuxSeries([d, 1/d, PuiseuxSeries(1)], z);
ObjectUnivariatePolynomialOverPowerSeriesObject,Array(0..2, {(1) = module PuiseuxSeriesObject () local pso, rays::(list(list(rational))), e::table, ord::list, ordCV::list, dstyle::(list(:-identical(:-maxterms, :-precision) = {:-nonnegint, :-identical(:-infinity)})), nonzero_pso_bound::{nonnegint, identical(undefined)}, smallest_term_bound::{nonnegint, identical(undefined)}, SublistLessThanV; option object; end module, (2) = module PuiseuxSeriesObject () local pso, rays::(list(list(rational))), e::table, ord::list, ordCV::list, dstyle::(list(:-identical(:-maxterms, :-precision) = {:-nonnegint, :-identical(:-infinity)})), nonzero_pso_bound::{nonnegint, identical(undefined)}, smallest_term_bound::{nonnegint, identical(undefined)}, SublistLessThanV; option object; end module})
e_factors := PuiseuxFactorize(e, returnleadingcoefficient=false):
for f in e_factors do print(Truncate(f, 4)); end do:
1y23
y43+y+y13z
1−y132−y2312−19y24−829y43720+y13z
Truncate(e - mul(f, f in e_factors), 6);
0
We see that the factors are linear or constant in z and they multiply to e.
Download Help Document
