Function Singularities and Their Applications

Adam Strzeboński, Wolfram Research

Join the Conversation #WolframTechConf

Abstract

Mathematica 11 has built-in knowledge about singularities and branch cuts of mathematical functions. In my talk I will discuss how this information is used in various parts of the system, including visualization, equation solving, root counting, global optimization, computation of limits, and definite integration.

Function Singularities

The built-in knowledge about “singularities” of mathematical functions includes the following information.

• Is the function locally analytic?

• Location of branch cuts.

• Location of “definition boundaries” -- these include: boundaries of the domain, boundaries of the set of real arguments for which the function is real-valued, and boundaries of piecewise cases.

• Location of poles.

• Location of essential singularities.

• Location of branch points.

The singularity sets are described as pairs {equation, condition}. The solution set of equation && condition includes the singularity set (it may be a strict inclusion).

The description can be propagated through composition of functions.

Sing(f((Subscript[g, 1]),…,(Subscript[g, n])))=Sing(f((Subscript[x, 1]),…,(Subscript[x, n])))(Subscript[|, (Subscript[x, 1])=(Subscript[g, 1]),…,(Subscript[x, n])=(Subscript[g, n])]) ⋃ Sing((Subscript[g, 1]))⋃…⋃Sing((Subscript[g, n]))

Poles of an inner function may become essential singularities or removable singularities of the composed function. FunctionSingularities assumes the “worst case scenario”, hence poles and removable singularities may be reported as essential singularities, and removable singularities may be reported as poles.

In Mathematica 11 there is no “official” direct interface, but the knowledge is used internally in various parts of the system.

Simplify`FunctionSingularities[Sin[1/x]+Log[x-1]+Tan[x]/(x-2)+InverseEllipticNomeQ[x],x]

Visualization

The knowledge of function singularities is used to automatically determine plot exclusions.

Mathematica 10:

Plot[Tan[x],{x,-3Pi,3Pi}]

Mathematica 11:

Plot[Tan[x],{x,-3Pi,3Pi}]

Mathematica 10:

Plot[Gamma[x],{x,-5,5}]

Mathematica 11:

Plot[Gamma[x],{x,-5,5}]

Mathematica 10:

Plot[Re[Gamma[1+Sqrt[Tan[x^2+1]]]],{x,-5,5}]

Mathematica 11:

Plot[Re[Gamma[1+Sqrt[Tan[x^2+1]]]],{x,-5,5}]

Mathematica 10:

Plot3D[Tan[Sin[x]+Log[y^2+1]],{x,-5,5},{y,-5,5},PlotStyle->Opacity[0.7]]

Mathematica 11:

Plot3D[Tan[Sin[x]+Log[y^2+1]],{x,-5,5},{y,-5,5},PlotStyle->Opacity[0.7]]

Solving univariate transcendental equations

Holomorphic functions

Reduce and related functions find roots of a univariate holomorphic function in a complex rectangle using the following heuristic.
1) Find bignum root approximations using a numeric method.
2) Use interval arithmetic (for elementary functions) or significance arithmetic (for non-elementary functions) and a criterion based on Rouche’s theorem to prove existence and multiplicities of roots in disks represented by the approximations.
3) Compute the winding number of the function along the boundary of the rectangle using interval arithmetic (for elementary functions) or numeric integration (for non-elementary functions) to prove that all roots have been found.

We can use the above heuristic, if we can show that the function is holomorphic in the given rectangle. We can do this using our knowledge of function singularities -- we can check that the rectangle does not intersect the branch cuts and does not contain poles or essential singularities.

sols=Solve[FresnelC[z]-z==1&& -2≤ Re[z]≤2&& -2≤Im[z]≤2, z]

f=FresnelC[x+I y]-x-I y-1;
Show[{ContourPlot[{Re[f],Im[f]},{x,-2,2},{y,-2,2},Contours->0],Graphics[{Red,PointSize[Large],Point[{Re[z],Im[z]}/.sols]}]}]

Meromorphic functions

If the function is meromorphic in the given rectangle we can use the same heuristic, provided that we show that the function is meromorphic and we can count the poles of the function in the rectangle.

The function is meromorphic if the rectangle does not intersect the branch cuts and does not contain essential singularities. FunctionSingularities provides equations and inequalities that poles need to satisfy. If we can solve these equations and find multiplicities of the poles, then we can determine the number of roots.

sols=Solve[Tan[z]-Log[z+3]==z^2&& -2≤ Re[z]≤2&& -2≤Im[z]≤2, z]

f=Tan[z]-Log[z+3]-z^2/.z->x+I y;
Show[{ContourPlot[{Re[f],Im[f]},{x,-2,2},{y,-2,2},Contours->0],Graphics[{Red,PointSize[Large],Point[{Re[z],Im[z]}/.sols]}]}]

Show[{Plot3D[Abs[f],{x,-2,2},{y,-2,2}],Graphics3D[{Red,PointSize[Large],Point[{Re[z],Im[z],0}/.sols]}]}]

Root counting

For counting roots of a transcendental function in a complex rectangle CountRoots uses FunctionSingularities in a similar way Solve does. First we need to prove that the function is meromorphic (the rectangle does not intersect the branch cuts and does not contain essential singularities). If the rectangle contains poles, we need to find the poles and compute their multiplicities. Then the number of roots is equal to the winding number plus the number of poles (counted with multiplicities). As in Solve, we compute the winding number using interval arithmetic (for elementary functions) or numeric integration (for non-elementary functions).

Example: elementary function

CountRoots[Tan[z^3+1]-E^z,{z,-3-3I,3+3I}]

f=Tan[z^3+1]-E^z/.z->x+I y;
Plot3D[Abs[f],{x,-3,3},{y,-3,3},PlotPoints->100]

Example: special function

CountRoots[BesselJ[2,z^2-1]-2z,{z,-3-3I,3+3I}]

f=BesselJ[2,z^2-1]-2z/.z->x+I y;
ContourPlot[{Re[f],Im[f]},{x,-3,3},{y,-3,3},Contours->0]

Example: periodic functions

CountRoots uses FunctionPeriod to recognize periodicity. For periodic functions it only needs to compute roots within one period and in the “remainder” rectangles.

CountRoots[Tan[x],{x,-10^10-I,10^10+I}]

CountRoots[JacobiCD[z/2,2],{z,0,10^10+10^10 I}]//AbsoluteTiming

Univariate optimization

Suppose we want to optimize a univariate function f over the solution set S of the given constraints, and suppose that we can symbolically compute the derivative g of f. It is not sufficient to find the roots of g in S, since g is the derivative of f only “generically”. More precisely, g is the derivative of f outside the set of singular points of f and g, where the singular points may be poles, essential singularities, branch points, endpoints of intersections of branch cuts with the real line, or boundary points of piecewise cases. We optimize f over the intersection of the real domain of f with S, so also need to know the boundary points of the real domain of f.

Example: maximum at a regular point

max = Maximize[{E^-Abs[x] - Erf[Tan[x]], x^2 <= 9}, x]

Show[{Plot[E^-Abs[x]-Erf[Tan[x]],{x,-3,3}],Graphics[{Red,PointSize[Large],Point[{x/.max[[2]],max[[1]]}]}]}]

Example: maximum at piecewise case boundary

max=Maximize[{2E^-Abs[x]-Erf[Tan[x]],x^2≤9},x]

Show[{Plot[2E^-Abs[x]-Erf[Tan[x]],{x,-3,3}],Graphics[{Red,PointSize[Large],Point[{x/.max[[2]],max[[1]]}]}]}]

Example: maximum at a singular point

max=Maximize[{ArcCos[Sin[5x^2]]-x^2,x^2≤4},x]

Show[{Plot[ArcCos[Sin[5x^2]]-x^2,{x,-2,2}],Graphics[{Red,PointSize[Large],Point[{x/.max[[2]],max[[1]]}]}]}]

Example: maximum at the boundary of the real domain

max=Maximize[{ArcCos[2Sin[5x^2]]-x^2,x^2≤4},x]

Show[{Plot[ArcCos[2Sin[5x^2]]-x^2,{x,-2,2}],Graphics[{Red,PointSize[Large],Point[{x/.max[[2]],max[[1]]}]}]}]

Limit computation

Here knowledge of function singularities is used to implement a very simple criterion -- if the function f is continuous at c then

(lim_(x->c))f(x)=f(c)

To show that f is continuous at c it suffices to show that c does not satisfy any of the conditions returned by FunctionSingularities. This can be done very fast compared to the general method used by Limit, which is based on computation of series.

Example: nested entire function example

This takes 119385. seconds in V11, but in the development version implementing the criterion it takes 0.00003 seconds on the same machine. This is a 3.9 billion times speedup.

Limit[Nest[Cos,x,17],x->0]//RepeatedTiming

Limit for 18 times nested Cos does not finish in a week in V11. With FunctionSingularities we can go up to about 20000.

f=Nest[Cos,x,20000];
(lim=Limit[f,x->0]);//Timing

N[lim,20]

The limitation here is not the speed of FunctionSingularities, but the size of function call stack needed by ReplaceAll.

f=Nest[Cos,x,10^6];
Simplify`FunctionSingularities[f,x]//Timing

Example: nested function with singularities

Here is an example with a function that does have singularities.

f=Nest[Tan,x,100];
LeafCount[Simplify`FunctionSingularities[f,x]]

(lim=Limit[f,x->1]);//Timing

N[lim,20]

The size for this example is limited by the speed of resolving numeric conditions with nested Tan. We can find the conditions for much larger examples.

f=Nest[Tan,x,10^4];
LeafCount[Simplify`FunctionSingularities[f,x]]//Timing

Example: composition of special functions

A “random” combination of special functions.

Limit[BesselJ[Gamma[Tan[Log[x]]],Exp[Beta[x,Erf[FresnelC[Sec[x^2-7]]]]]],x->2]//RepeatedTiming

N[%,20]

Limit[BesselJ[Gamma[Tan[Log[x]]],Exp[Beta[x,Erf[FresnelC[Sec[x^2-7]]]]]],x->2]//Timing

N[%,20]

Integration

Future projects making use of knowledge of function singularities

• Computation of definite integrals using the Newton-Leibniz theorem.

• Computation of Cauchy Principal Value integrals.

• Computation of Hadamard Finite Part integrals.