Symbolic Computation with Roots of Transcendental Functions

Adam Strzebonski

Wolfram Research Inc.

Polynomial equations

Linear equations can be solved without introducing any additional functions.

Reduce[2 x+(3/7)==0,x]

To express solutions of polynomial equations with degrees ≤ 4 we need two additional functions

(Sqrt[_]),((_)^(1/3)):ℂ ⟶ ℂ

Reduce[(x^2)+2 x+3==0,x]

Reduce[(x^4)+2 x+3==0,x,Quartics->True]//First

To solve arbitrary polynomial equations we need to introduce an infinite family of functions

Root:(Underoverscript[⊔, (d=3), ∞])(Subscript[ℂ[x], d])×{1,…,d} ⟶ ℂ

or equivalently

Root:(Underoverscript[⊔, (d=3), ∞])(ℂ^(d+1))×{1,…,d} ⟶ ℂ

Reduce[(x^5)+2 x+5==0,x]

Transcendental (locally holomorphic) equations

A locally holomorphic fuction has a countable number of roots, so in principle we could use a representation for roots of transcendental function which would be very similar to what we have in the polynomial case

Root:𝒪[x]×ℕ ⟶ ℂ

The problem with this approach is that we have no natural canonical way of assigning indices to roots of transcendental functions, and in most cases we can only find roots in bounded regions. For this reason our representation uses approximations to identify roots. An approximate number represents a circle in the complex plane (or a real interval). Such representation of an exact root of f is valid if we can prove there is exactly one root of f in the region corresponding to the given approximate number.

Root[{f, a}] represents the root of the univariate pure function f in the region represented by the approximate number a. f has been proven to have exactly one simple root in the region represented by a.

Example: Kepler's equation

Reduce[x-((Sin[x])/2)==1&&0≤x≤2,x]

r=%[[2]]

NumericQ[r]

Precision[r]

N[r,100]

Since a root is identified by its numeric approximation, transcendental Root objects cannot contain symbolic parameters. We can still graph the solution of Kepler's equation by defining a function which solves the equation at each step. (There are of course more efficient way of graphing solutions to Kepler's equation, because this is a well known and well investigated equation. However, the method involving transcendental Root objects could be applied to graphing roots of any analytic equation in a bounded interval.)

f[M_?NumberQ,e_?NumberQ]:=x/.ToRules[Reduce[Rationalize[x-e Sin[x]==M,0]&&0≤x≤20,x]]

Plot[Evaluate[Table[f[M,e],{e,0,0.8,0.2}]],{M,0,4Pi}]

Example: The n-dimensional unit ball with maximal volume

The volume of an n-dimensional ball is equal to (((π^(n/2)) (R^n) )/(Γ((n/2)+1))). This illustrates how the volume of the n-dimensional unit ball changes with n.

vol[n_]:=((π^(n/2))/(Gamma[(n/2)+1]))
ListPlot[vol/@Range[20]]

Treating the volume as a function of a real parameter n, we can compute the maximum.

Maximize[{vol[n],1≤n≤20},n]

Example: Planck's radiation function

Planck' s radiation function is the function

f[x_]:=(15/((π^4) (x^5) ((E^(1/x))-1)))

p=Plot[f[x],{x,0,2}]

Find the maximum point.

m=Reduce[f'[x]==0,x,Reals]

Find the inflection points.

i=Reduce[f''[x]==0,x,Reals]

Show[{p,Graphics[{{Red,PointSize[0.02],Point[{x,f[x]}/.{ToRules[N[m]]}]},{Green,PointSize[0.02],Point[{x,f[x]}/.{ToRules[N[i]]}]}}]}]

Wavelength with the greatest radiance as a function of temperature.

f[λ_,T_]:=((2 h (c^2))/((λ^5) ((E^((h c )/(λ k T)))-1)));
k=1.3806504 (10^(-23));
h=6.62606896 (10^(-34));
c=299792458;

λ[T_?NumberQ]:=x/.ToRules[Reduce[Rationalize[D[f[x,T],x],0]==0,x,Reals]]

Plot[λ[T],{T,1000,10000},AxesOrigin->{0,0}]

Multiple roots

Unlike in the polynomial case, there are no algorithms for square free factorization of arbitrary transcendental functions. Hence a transcendental equation we are trying to solve may have multiple roots and we do not have a way to distinguish between a multiple roots and a closter of close but distinct roots. What we may be able to prove is that a function has a specific number of roots, counted with multiplicities, in a specified region.

Root[{f, a, mult}] represents a root of the univariate pure function f in the region represented by the approximate number a. f has been proven to have exactly mult roots, counted with multiplicities, in the region represented by a.

This equation has three double roots.

Reduce[409+120 x+400 (x^2)-400 ((Cos[x])^2)-120 Sin[x]-800 x Sin[x]==0&&Abs[x]<1,x]

r=%[[3,2]]

NumericQ[r]

N[r,100]

Root clusters

A Root object with multiplicity higher than one may in fact represent a cluster of close roots.

f=TrigReduce[(Sin[x]-x-1-I 10^-100)(Sin[x]-x-1+I 10^-100)];
r=Reduce[f==0&&Abs[x-1]<2,x][[2,2]]

When a higher precision approximation of the Root object is requested, it turns out that the Root object represents two distinct roots. In this case Mathematica gives a warning message, picks one of the roots and gives its approximation.

N[r,100]

The choice of root is remembered and the same root is used to compute higher order approximations.

N[r,200]

Examples: Equations solvable in terms of transcendental Root objects

Exp-log equations over the reals.

Reduce[((E^x)/(x^3))-(x+1) Log[x]-Log[Log[x]]==1,x,Reals]

Exp-log inequalities over the reals.

Reduce[x Log[x]-(E^x)>-2,x,Reals]

Elementary function equations over bounded intervals.

Reduce[Sin[Cos[x]]-4 (x^2)+(E^x)==1&&-1<x<1,x]

Elementary function equations over bounded complex regions.

Reduce[(((x+1))^(Sin[x]))==Cos[(x^2)-Tan[x]]+(I/2)&&0≤Re[x]≤1&&0≤Im[x]≤1,x]

Holomorphic function equations over bounded intervals. For nonelementary functions we have to rely on NIntegrate to find the number of roots.

Reduce[Cos[x]-BesselJ[5,x]==(1/2)&&0≤x≤10,x]

Holomorphic function equations over bounded complex regions.

Reduce[Gamma[x]-Log[x]==1&&Abs[x-2]<(3/2),x]

Here Reduce found some solutions but was not able to prove there are no other solutions.

Reduce[Zeta[x]==(E^(I x))-1&&Abs[x]<5,x]

Examples: Periodic equations

Periodic analytic function equations over the reals.

Reduce[(E^(Sin[x]))-Cos[3 Sin[x]]==0,x,Reals]

An equation involving a holomorphic function with a purely imaginary period over a vertical stripe in the complex plane.

Reduce[AiryAi[(E^x)]==3 (E^(-x))+1&&0≤Re[x]≤2,x]

An equation involving a doubly-periodic locally holomorphic function over the complex plane. Reduce cannot prove completeness here, because the function has poles within each period region.

Reduce[JacobiCS[x,3]==1,x]

Examples: High degree sparse algebraic equations

Polynomial Root objects are restricted to degrees not exceeding

$MaxRootDegree

Creation of a polynomial Root object requires factorization of its defining polynomial. For high degrees computing factorization would be too slow.

Creation of a transcendental Root object does not require factorization of its defining polynomial and its real roots can be isolated using the exp-log root isolation algorithm.

Reduce[(x^1000001)-5 (x^2)+1==0,x,Reals]

Reduce[(Sqrt[x])-3 (x^(1/300000))+5 (x^2)==0,x,Reals]

Examples: Exact optimization

Minimize an exp-log function.

Minimize[{(E^(x (E^x)))+(x^(Log[x])),x>0},x]

N[%,30]

Maximize an elementary function over a bounded interval.

Maximize[{Sin[Cos[x]]-(E^(((x-1))^2)),0<x<2},x]

Maximize a periodic elementary function.

Maximize[Sin[Cos[x]]-Cos[Sin[x]],x]

Maximize an analytic function over a bounded interval.

Maximize[{Erf[x]-((E^x)/(x^2)),1≤x≤3},x]

Example: Piecewise definite integration

Integrate a piecewise function with an exp-log inequality condition.

(∫_0^∞)x Boole[(E^x)-(x^x)>2]\[DifferentialD]x

Direct specification of transcendental Root objects

Transcendental Root objects can be entered directly.

FindRoot[Sin[x]-Gamma[x],{x, 1}]

In this case we are able to prove there is exactly one root equal to the approximation.

Root[{Sin[#]-Gamma[#]&,%[[1,2]]}]

The result is an exact numeric object. We can compute its approximation to an arbitrary precision.

N[%,100]

The multiplicity is computed automatically.

Root[{Cos[#]+1&,N[Pi,20]}]

If it can prove there are no roots equal to the specified approximation, Root gives an error message and returns unevaluated.

Root[{Sin[#1]-Gamma[#1]&,1`20}]

The result is not numeric.

NumericQ[%]

Here we cannot prove or disprove the existence of a root.

r=Root[{Abs[#1]-(2^(1/3))&,N[(((-2))^(1/3))]}]

The result is not numeric.

NumericQ[%]

Algorithm: Real root isolation for exp-log functions

For exp-log functions, i.e. arbitrary compositions of Plus, Times, Power and Log we have a real root isolation algorithm. It allows to find the real roots and the intervals on which the function is identically zero, assuming that we have a zero test for constants, including constants involving transcendental Root objects. In practice we may give a message about some expressions being assumed to be zero (this is already done in other situations in Reduce, GroebnerBasis, Limit etc.).

We do not have a square free factorization algorithm for exp-log functions hence we are not able to distinguish between multiple real roots, clusters of close real roots and clusters of complex roots close to the real axis. In our algorithm this problem appears as the problem of deciding whether a function is zero at a root of its derivative, and is solved by use of heuristic zero testing (with unproven assumptions returned in a warning message). In particular, multiple real roots of exp-log functions are represented as simple roots of elements of their Fourier sequences.

Algorithm: Finding roots of analytic functions in bounded regions

For arbitrary holomorphic functions we can only find roots in bounded regions of the complex plane. There are three main components to this algorithm, and each of the three components may fail. First, we need to be able to count the number of roots in the region. For elementary functions we can do it in a provably correct way by using complex interval arithmetic to compute the winding number. For other functions we need to rely on correctness of numeric integration (in this case Reduce gives a warning message). Next, we need to find approximations of all roots in the region. This is done using FindRoot with starting points provided either by computational geometry (at the moment we use post-processing of ContourPlot, Ulisses is going to write a more direct method) or by the Kravanja-Van Barel algorithm. Finally, we need to validate the approximations. We use a criterion based on Rouche's theorem. The criterion allows to prove that a given circle (represented by an approximate number) contains exactly a certain number of roots, multiplicities counted. Again, the criterion requires complex interval arithmetic, so the results are only provably correct for elementary functions. For other functions correctness depends on correctness of error bounds given by significance arithmetic (by default we use approximations with at least 30 digits of precision, so significance arithmetic is expected to give good error bounds). The algorithm fully succeeds if the number of validated roots, counted with multiplicities, equals the number of roots computed in the first step. If we find and validate some roots, but their number is less than the computed total number of roots or we failed to compute the total number of roots, Reduce gives a message and returns the roots it has found.