∫   ∑  ∇ γ  

Symbolic Computation with Mathematica

Adam Strzebonski, Devendra Kapadia, Itai Seggev, Oleksandr Pavlyk

Wolfram Research, Inc.

Plan for the Workshop

The following is a broad plan for topics to be discussed in this workshop:

Algebra and Optimization

Calculus Operations

Vector Calculus and Special Functions

Problems from Calculus

Formulation of problems for algebraic solvers

Main classes of problems that are solvable algorithmically

• Systems of linear equations (linear algebra)
• Systems of polynomial equations and inequations in complex variables (Groebner bases)
• Systems of polynomial equations and inequalities in real variables (cylindrical algebraic decomposition)
• Real univariate equations and inequalities involving tame elementary functions (real root isolation)

Mathematica solvers also implement several algorithms for solving certain classes of Diophantine problems (linear, binary quadratic equations, Thue equations, etc.). For other systems Mathematica tries various heuristics in order to reduce the problem to a combination of solving polynomial systems and using known solutions to simple transcendental problems.

Example from StackExchange

f=y E^(x-z);
g=9 x^2+4 y^2+36 z^2-36;
h=x y+y z-1;

Find the stationary points of f over the constraints g==0 ∧h==0.

Direct solution (unsuccessful)

gradf=Grad[f,{x,y,z}];
gradg=Grad[g,{x,y,z}];
gradh=Grad[h,{x,y,z}];
TimeConstrained[sol=Solve[gradf==L gradg+M gradh && g==0 && h==0,{x,y,z,L,M},Reals],300]

Two-step solution (successful)

eq1=Eliminate[gradf==L gradg+M gradh && g==0 && h==0,{L,M}]

sol1=Solve[eq1,{x,y,z},Reals];

y E^(x-z)/.sol1//N

Polynomial formulation

gradf

Since stationary point conditions depend only on directions of gradients and not on their magnitudes, one can divide gradients by nonzero functions.

sol2=Solve[gradf/E^(x-z)==L gradg+M gradh && g==0 && h==0,{x,y,z,L,M},Reals];

y E^(x-z)/.sol2//N

Formulation of problems for exact optimization

Main classes of optimization problems that are solvable algorithmically

• Linear objective functions and constraints (linear programming)
• Polynomial objective functions and constraints (cylindrical algebraic decomposition)
• Real univariate objective functions and constraints involving tame or periodic elementary functions (differentiation/limits/real root isolation)
• Some analytic objective functions and analytic equation constraints in a bounded box (Lagrange multipliers/equation solving)

For other systems Mathematica tries various heuristics in order to reduce the problem to a solvable form.

Example

Find the distance from the point p={3,3,3} to the torus given by the parametrization

p={3,3,3};
par={Cos[a](2+Cos[b]),Sin[a](2+Cos[b]),Sin[b]};

ParametricPlot3D[Evaluate[par],{a,0,2Pi},{b,0,2Pi}]

Direct solution (unsuccessful)

MinValue[{Norm[par-p],0≤a<2Pi&&0≤b<2Pi},{a,b}]

Reduce to analytic objective function in a bounded box

The “analytic objective function in a bounded box” method is not used for the last input, because the use of square root in the objective function adds an inequality constraint. The constraint is always true, but the method is used only if there are no inequalities other than box bounds (checked syntactically).

MinValue[{(par-p).(par-p),0≤a<2Pi&&0≤b<2Pi},{a,b}]//Timing

Specifying the bounds on a and b is necessary.

MinValue[(par-p).(par-p),{a,b}]

Reduce to a polynomial problem

eqs=Join[Thread[{x,y,z}-par],{Sin[a]^2+Cos[a]^2-1,Sin[b]^2+Cos[b]^2-1}];
GroebnerBasis[eqs,{x,y,z},{Sin[a],Cos[a],Sin[b],Cos[b]}]//Timing

poly=%[[2,1]]

MinValue[{Norm[{x,y,z}-p],poly==0},{x,y,z}]//Timing

Region formulations

Trigonometric parametric description with no parameter bounds

The heursitcs successfully reduce the problem, but it takes a long time.

R1=ParametricRegion[par,{a,b}]

RegionDistance[R1,{3,3,3}]//Timing

N[%]

Trigonometric parametric description with parameter bounds

R2=ParametricRegion[{par,0≤a<2Pi&&0≤b<2Pi},{a,b}]

RegionDistance[R2,{3,3,3}]//Timing

N[%]

Polynomial implicit description

R3=ImplicitRegion[poly==0,{x,y,z}]

RegionDistance[R3,{3,3,3}]//Timing

N[%]

Solving for real variables vs. solving over the reals

Solve[(Sqrt[x]) (Sqrt[y])==(Sqrt[z])&&(Sqrt[x-y]) (Sqrt[x+y])==2 (Sqrt[z-1])&&(x^3)+(y^3)+(z^3)==7&&(x|y|z)∈Reals,{x,y,z}]//Timing

Solve[(Sqrt[x]) (Sqrt[y])==(Sqrt[z])&&(Sqrt[x-y]) (Sqrt[x+y])==2 (Sqrt[z-1])&&(x^3)+(y^3)+(z^3)==7,{x,y,z},Reals]//Timing

Symbolic Integration

Scope of Integrate

Integrate can compute:

A wide variety of  indefinite integrals

Definite integrals over finite and infinite intervals

Multiple integrals in an arbitrary number of dimensions

Integrals of piecewise functions and distributions

Integrals over regions, including areas and volumes

Integral of a Rational Function

The indefinite integral of this rational function has three parts - polynomial, rational, and logarithmic.

∫((3+15 x+37 (x^2)+36 (x^3)+14 (x^4))/((((1+x))^2) (4+7 x)))\[DifferentialD]x

We can use D to verify the result.

D[%,x]//Together

A Rational-Exponential integral

Now consider the following rational-exponential indefinite integral.

∫(((E^x) ((E^x) x-8 (E^(2 x)) x-15 (E^(3 x)) x-3 (E^(4 x)) x+3 (E^(5 x)) x+4 x-5 (E^x)-(E^(2 x))+(E^(3 x))-3))/((-3+(E^x)) (((1+(E^x)))^2) x))\[DifferentialD]x

To guess the answer, we replace exponentials by powers of a new variable t and then apply Apart to the integrand.

Apart[(((E^x) ((E^x) x-8 (E^(2 x)) x-15 (E^(3 x)) x-3 (E^(4 x)) x+3 (E^(5 x)) x+4 x-5 (E^x)-(E^(2 x))+(E^(3 x))-3))/((-3+(E^x)) (((1+(E^x)))^2) x))/. {(E^(a_. x)):>(t^a)}]

Exercise: Use the second argument of Apart to get the desired output for the above.

After evaluation, the integral has four parts - polynomial, fractional, logarithmic and transcendental.

∫(((E^x) ((E^x) x-8 (E^(2 x)) x-15 (E^(3 x)) x-3 (E^(4 x)) x+3 (E^(5 x)) x+4 x-5 (E^x)-(E^(2 x))+(E^(3 x))-3))/((-3+(E^x)) (((1+(E^x)))^2) x))\[DifferentialD]x

Once again, we can use D to verify the result.

D[%,x]//Together

A Gaussian Integral

The following Gaussian definite integral has an elegant closed form answer.

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

N[%]

We can use NIntegrate to verify the result.

NIntegrate[(x^12) (E^(-(x^2))),{x,0,∞}]

Exercise: Derive the symbolic result using differentiation under the integral sign.

Exact Solutions of Differential Equations

Scope of DSolve

DSolve can find exact solutions for the following:

Ordinary differential equations (linear, nonlinear, single, system)

Partial  differential equations  (general first-order PDEs, some higher order equations)

Differential-algebraic equations (linear constant coefficient systems)

Delay differential equations.

Hybrid differential equations (ODEs with events).

The Aging Spring

This ODE is a model for an aging spring. It has an exponential function in its coefficient and is solved by transforming it to Bessel’s equation.

AgingSpring=x''[t]+4 (E^(-(t/4))) x[t]==0;

The solution is expressed in terms of Bessel functions.

sol=DSolve[{AgingSpring,x[0]==1,x'[0]==2},x,t]

The solution can be verified as follows.

{AgingSpring,x[0]==1,x'[0]==2}/.sol[[1]]//FullSimplify

Here is a plot of the solution.

Plot[Evaluate[x[t]/. sol[[1]]],{t,0,30},PlotStyle->{Thickness[0.0075],Red},Background->LightYellow]

Ball in a Square Box

The following models a ball in a square box, that changes direction upon impact with the side walls. Try varying the second argument of Rationalize.

sol=DSolve[{x'[t]==a[t],y'[t]==b[t],x[0]==0,y[0]==0,a[0]==1,b[0]==Rationalize[N[(Sqrt[2])],0.1],WhenEvent[((x[t])^2)==1,a[t]->-a[t]],WhenEvent[((y[t])^2)==1,b[t]->-b[t]]},{x,y},{t,0,100},DiscreteVariables->{a,b}];

ParametricPlot[{x[t],y[t]}/. sol,{t,0,100},Frame->True,FrameTicks->None,PlotRange->1,Axes->False,PlotPoints->300]

Symbolic Summation

Scope of Sum

Sum can compute:

A wide variety of  indefinite sums

Definite sums over finite and infinite intervals

Multiple sums in an arbitrary number of dimensions

Sums of piecewise functions, q-functions, etc.

Procedural vs Symbolic Summation

The sum of the squares of the first thousand natural numbers can be found as follows:

(∑_(k=1)^(10^3))(k^2)

The same answer can be obtained using symbolic summation:

(∑_(k=1)^n)(k^2)

%/. {n->(10^3)}

There is a built-in  threshold of (10^6) beyond which symbolic summation is typically invoked. This threshold can be controlled using Method -> “Procedural”/”Symbolic”.

Infinite Hypergeometric Series

Hypergeometric functions often have simple expressions for certain discrete ranges of their parameters.

Hypergeometric2F1[a,b,c,x]

Exercise: Find a value of the parameter ‘a’ for which the above expression is a polynomial in ‘x’.

Hypergeometric simplification is applied internally for evaluating series such as the following:

(∑_(n=0)^∞)(((((2/(5E)))^n) Binomial[2 n,n])/( (n+3)))

Before proceeding, let us verify the result numerically:

N[%]

NSum[((2/5)^(n )*Binomial[2*n, n])/(E^n*(n + 3)), {n, 0, Infinity}]

The internal conversion to HypergeometricPFQ may be seen as follows:

Sum[((2/5)^(n )*Binomial[2*n, n])/(E^n*(n + 3)), {n, 0, Infinity},Method->"InactivePFQ"]//TraditionalForm

Activate[%]

The Calculus Navigator

A poster for navigating Calculus functions in Mathematica is available here.

Vector Calculus

Calculus on Fields: Grad, Div, Curl, Laplacian  (9+)

Calculus on Parametric Regions: ArcLength, ArcCurvature, FrenetSerretSystem, Area, Volume (10+)

Coordinate Systems: CoordinateChartData, CoordinateTransformData, TransformedField, CoordinateTransform (9+)

Structure of a coordinate system:     ({{{, metric, ,, chart, ,, dim , }}, { , typically omitted if in Euclidean space,  , always needed,  , typically omitted as it can usually be inferred from other arguments,  }}), 

where metric and chart can each be either “name” or {“name”, params}.

Simple example (though note answer is an orthnormal basis):

Grad[f[r,θ], {r,θ},"Polar"]

Here we need to specify the dimension:

CoordinateChartData[{"Confocal",2}, "ParameterRangeAssumptions"]

A Euclidean chart with a parameter specified

ArcLength[{t,(t^2),(Sqrt[t])}, {t,.5,1.5}, {{"ProlateSpheroidal",{1}}}]

Assumptions

Integrate-like functions—ArcLength, Area, Volume—take and can really benefit from the Assumptions option.

D-like functions—Grad, Div, Curl, Laplacian, ArcCurvature, FrenetSerretSystem, TransformedField—do not take an Assumptions option.

Problem: Simplify the following formula for the length of a sine curve.

ArcLength[a Sin[b x],{x,c,d}]

Solution/Discussion

Mathematica defaults to assuming that symbolic parameters are complex.  When there are complicated conditions on real and imaginary parts, a good place to start is to assert that the parameters are real:

ArcLength[a Sin[b x],{x,c,d}, Assumptions->(a|b|c|d)∈Reals]

In this particular, we only need to put assumptions and a and b.  This can be seen from the condition d<0 || c>0, which implies that c and d are real.

ArcLength[a Sin[b x],{x,c,d}, Assumptions->(a|b)∈Reals]

Good assumptions can be so important to integration that ArcLength (and Area/Volume) will attempt to add useful assumptions automatically:

ArcLength[f[x],{x,c,d}]

ArcLength[{f[t],g[t]},{t,c,d}, {{"Elliptic",a}}]

Note that the condition that the endpoints are order isn’t just useful, it’s necessary for a correct answer:

ArcLength[f[x],{x,1,0}]

Homework (Needs 10.0.2): Simplify the following formula for the area of an ellipse with semi-major axes of length a and b.

Area[{r,v},{r,0,ArcSech[(Sqrt[1-((b^2)/(a^2))])]}, {v, 0, 2Pi},{{"Elliptic",(Sqrt[(a^2)-(b^2)])}}]

Homework solution

The presence of Abs[a] and Abs[b] strongly suggest the need to add the conditions a>0&& b>0 (they are lengths):

Area[{r,v},{r,0,ArcSech[(Sqrt[1-((b^2)/(a^2))])]}, {v, 0, 2Pi},{{"Elliptic",(Sqrt[(a^2)-(b^2)])}},Assumptions->a>0&&b>0]

Your exact starting point can affect the needed assumptions.  For example, here we only need the assumption b>0 (why?):

Area[{r,v},{r,0,ArcSech[(Sqrt[1-((b^2)/(a^2))])]}, {v, 0, 2Pi},{{"Elliptic",a(Sqrt[1-((b^2)/(a^2))])}},Assumptions->b>0]

Had we used a direct parametrization of an ellipse, our starting point would have been much simpler:

Area[{a  r Cos[θ], b  r Sin[θ]}, {θ,0,2Pi}, {r,0,1}]

But note that even this needs the assumptions on a and b to fully simplify:

Area[{a r Cos[θ], b r Sin[θ]}, {θ,0,2Pi}, {r,0,1},Assumptions->a>0&&b>0]

Finally, had we entered a disk region with two radii, the assumptions that a and b are positive would have been inferred from their meaning as radii:

Area[Disk[{0,0},{a,b}]]

Spherical Harmonics

Problem: Show that spherical harmonics are eigenfunctions of the Laplacian

It is a well known fact that the spherical harmonics are eigenfunctions of the Laplacian on the sphere.  That means

(Underscript[(∇^2), sphere])(Subscript[Y, l]^m)(θ,φ)==constant (Subscript[Y, l]^m)(θ,φ)

One way to show this is to divide by the Laplacian by the original function.  For example:

Table[((Laplacian[SphericalHarmonicY[l,0,θ,φ],{θ,φ},{"Standard",{"Sphere",1}}])/(SphericalHarmonicY[l,0,θ,φ])),{l,0,4}]//Simplify

However, the Laplacian looks like this:

Laplacian[SphericalHarmonicY[l,m,θ,φ],{θ,φ},{"Standard",{"Sphere",1}}]

And dividing by (Subscript[Y, l]^m)(θ,φ) doesn’t simplify it very much (as in, at all):

%/SphericalHarmonicY[l,m,θ,φ]

How can this be simplified?

Solution/discussion

With no assumptions, FullSimplify does a reasonable job:

ClearSystemCache[]
AbsoluteTiming[noAssums=TimeConstrained[FullSimplify[Laplacian[SphericalHarmonicY[l,m,θ,φ],{θ,φ},{"Standard",{"Sphere",1}}]/SphericalHarmonicY[l,m,θ,φ]],30]]

If we just put assumptions on the angles, it will take a really long time and I wasn’t willing to wait for more than 10 minutes:

ClearSystemCache[]
AbsoluteTiming[angleAssums=TimeConstrained[FullSimplify[Laplacian[SphericalHarmonicY[l,m,θ,φ],{θ,φ},{"Standard",{"Sphere",1}}]/SphericalHarmonicY[l,m,θ,φ],0<θ<Pi&&-Pi<φ<Pi],600]]

If we just put in the conditions on l and m, the run time increases but a factor of 3, but all we get is from the effort is a common demoniator:

ClearSystemCache[]
AbsoluteTiming[lmRangeAssums=TimeConstrained[FullSimplify[Laplacian[SphericalHarmonicY[l,m,θ,φ],{θ,φ},{"Standard",{"Sphere",1}}]/SphericalHarmonicY[l,m,θ,φ],-l<m<l],600]]

ClearSystemCache[]
AbsoluteTiming[lmRangeIntegerAssums=TimeConstrained[FullSimplify[Laplacian[SphericalHarmonicY[l,m,θ,φ],{θ,φ},{"Standard",{"Sphere",1}}]/SphericalHarmonicY[l,m,θ,φ],(l|m)∈Integers&&-l<m<l],600]]

Lest you think the assumption is {l,m}∈ℤ is superfluous, the (Subscript[Y, m]^l) are perfectly analytic in l and m:

SphericalHarmonicY[(Sqrt[2]),.35,Pi/2,0]

But if we combine the assumptions, we get something that returns quickly and is simpler:

ClearSystemCache[]
AbsoluteTiming[lmRangeAngleAssums=FullSimplify[Laplacian[SphericalHarmonicY[l,m,θ,φ],{θ,φ},{"Standard",{"Sphere",1}}]/SphericalHarmonicY[l,m,θ,φ],-l<m<l&&0<θ<Pi&&-Pi<φ<Pi]]

ClearSystemCache[]
AbsoluteTiming[fullAssums=FullSimplify[Laplacian[SphericalHarmonicY[l,m,θ,φ],{θ,φ},{"Standard",{"Sphere",1}}]/SphericalHarmonicY[l,m,θ,φ],(l|m)∈Integers&&-l<m<l&&0<θ<Pi&&-Pi<φ<Pi]]

In this particular case, the assumption the l and m are integers was not particularly helpful:

lmRangeAngleAssums===fullAssums

Bonus: how to get the desired answer

We need to apply the recurrence relation obeyed by LegendreP at the right moment.  That moment is the answer we got with the full assumptions.  The function Simplify`SimplifyRecurrence will do this for us:

Simplify`SimplifyRecurrence@fullAssums

We can now simplify down to an obvious constant:

Simplify[%,(l|m)∈Integers&&-l<m<l&&0<θ<Pi&&-Pi<φ<Pi]//TraditionalForm

Homework

Assuming the functions (Subscript[f, nlm])[r,θ,φ] :=(Subscript[R, nl])[r]SphericalHarmonicY[l,m,θ,φ] are eigenfunctions of the Laplacian in three-dimensional space, find the form of the functions (Subscript[R, nl[r]]).

Homework solution

We wish to show (((∇^2)(Subscript[f, nlm])[r,θ,φ])/( (Subscript[f, nlm])[r,θ,φ]))=λ for some constant λ which might depend on n, l, and m.  So we begin by substituing a generic function f[r,θ,φ]:

eqn = Expand[Laplacian[f[r,θ,φ],{r,θ,φ},"Spherical"]/f[r,θ,φ]==λ]

At this point we can either use brute force or be clever.

Brute force:

The obvious next step is to substitute in our form for the function f.  Because there are derivatives, it is important that we substitute it in a as a function f rather than as an expression f[r,θ,φ].  We therefore use pure function notation with {r,θ,φ} replaced by {#1,#2,#3}:

eqn/. f->(R[#1]SphericalHarmonicY[l,m,#2,#3]&)

We now have something that looks like similar to the spherical harmonic equations we saw above.  We therefore try to simplify it:

FullSimplify[%,(l|m)∈Integers&&-l<m<l&&0<θ<Pi&&-Pi<φ<Pi]

Progress, but we still haven’t gotten all the way to something that only depends on r.  So let us try SimplifyRecurrence followed by another Simplify as was done above:

Reqn =Simplify[Simplify`SimplifyRecurrence[%],(l|m)∈Integers&&-l<m<l&&0<θ<Pi&&-Pi<φ<Pi]

Or, rearranged to its traditional form, we see a spherical Bessel equation:

MapAt[Expand,- Simplify[First[Reqn] r R[r]],2]==0

As it is a well known 2nd-order ODE, DSolve can quickly return an answer in terms of known special functions:

DSolve[Reqn,R[r],r]

Spherical Bessel functions are really just particular combinations of “ordinary” Bessel functions.

Simplify[%//FunctionExpand]

Clever:

If we’re clever, we recognize that the first three terms are (1/(r^2)) times the eigenfunction equation for the Laplacian on the unit sphere:

First[eqn]-(1/(r^2))((Laplacian[f[r,θ,φ],{θ,φ},{"Standard",{"Sphere",1}}])/(f[r,θ,φ]))//Expand

Equivalently, the first three terms are the eigenfunction equation for the Laplacian on a sphere of radius r:

First[eqn]-((Laplacian[f[r,θ,φ],{θ,φ},{"Standard",{"Sphere",r}}])/(f[r,θ,φ]))//Expand

The function (Subscript[R, nl])will just come along for the ride and cancel on those three terms, and we’ll get the eigenvalue -l(l+1) found above.  Thus, we can rewrite eqn as

eqnNew = ((-l(l+1))/(r^2))+((2 Derivative[1,0,0][f][r,θ,φ])/(r f[r,θ,φ]))+((Derivative[2,0,0][f][r,θ,φ])/(f[r,θ,φ]))==λ

We know that all θ,φ dependence has been removed, so we can just replace f with R, again using pure function notation so that derivatives are computed correctly.

eqnNew /. f->(R[#]&)

This is (1/r) times the equation Reqn found above, which DSolve has no problem dealing with.

DSolve[%,R[r],r]

Some Problems from Calculus

P1: Evaluation of (∫_0^(π/2))log((sin^2)(x))dx

Source: math.stackexchange.com/questions/58654

Built-in

Integrate[Log[Sin[x]^2],{x,0,Pi/2}]

Derivation

Make a substitution Sin[x]==t

Log[Sin[x]^2]Dt[x]/.x->ArcSin[t]

Or use a nice function built especially for this task.

int=ChangeIntegrationVariables[(InactiveIntegrate[Log[((Sin[x])^2)], x, 0, (π/2)]),t==Sin[x],t]

Representing Log[t] as a parametric derivative (Feynman trick ?)

Limit[D[(t^(s-1)),s],s->1]

Assuming[s≥1,Activate[int/.Log[t]->(t^(s-1))]]

Limit[D[%,s],s->1,Direction->-1]

FunctionExpand[%]//Simplify

P2: Asymptotic expansion of (∫_n^∞)exp(-t) (t^n)dt

Source: math.stackexchange.com/questions/327525
Source: problem 11353 of AMM
Source: math.stackexchange.com/questions/71447

Built-in

Integrate[t^n Exp[-t],{t,n,Infinity},Assumptions->n>0]

Series[%,{n,Infinity,3},Assumptions->n>100]

Derivation

Change variables t==n (u+1).

int=Assuming[n>0,ChangeIntegrationVariables[Inactive[Integrate][t^n Exp[-t],{t,n,Infinity}],t==n+n u,u]]

factor=Exp[-n]n^(n+1);

int=Assuming[n>0,PullConstantOut[int,factor]]

Apply Watson’s lemma. The integrand attains maximum at the lower integration limit.

ϕ[u_]=FunctionExpand[Log[Exp[-n u](((1+u))^n)],u>0&&n>0]

{sol}=Solve[ϕ'[u]==0,u]

Develop the series exponent into series about the maximum.

ser=FullSimplify[Series[ϕ[u],{u,0,20}],n>0]

quadratic=Normal[ser+O[u]^3]

rest=ser-quadratic

Represent the integral as (∫_0^∞)Exp[-n u^2/2]((∑_k)(Subscript[c, k])(u^k))du

coeffs=Expand[CoefficientList[Exp[rest],u]]

ints=Table[Integrate[u^k Exp[quadratic],{u,0,Infinity},Assumptions->n>0],{k,0,Length[coeffs]-1}]

Series[Accumulate[coeffs ints],{n,Infinity,4}]//Differences//Take[#,-4]&

res=factor Series[coeffs.ints,{n,Infinity,3}]

P3: Evaluation of (∫_0^∞)(((tan^(-1))(a (sin^2)(x)))/(x^2))dx

Source: math.stackexchange.com/questions/174258

Built-in:

(∫_0^∞)((ArcTan[a ((Sin[x])^2)])/(x^2))\[DifferentialD]x

Derivation:

Key: develop a series with respect to parameter, integrate term-wise.

term=Assuming[n∈Integers&&n≥0,SeriesCoefficient[ArcTan[a Sin[x]^2]/x^2,{a,0,2n+1}]//FullSimplify]

Evaluate first few terms.

Table[{n,Integrate[term,{x,0,Infinity}]},{n,0,12}]

Guess a closed form.

integratedTerm[n_]=FindSequenceFunction[%,n]

Resum.

Sum[integratedTerm[n] (a^(2n+1)),{n,0,Infinity}]

P4: Find (lim_(x->∞)) ((∑_(n=1)^∞)(((((-1))^n) (x^(2 n-1)))/((2 n)! log(2 n))))

Source: math.stackexchange.com/questions/191008

Built-in:

Limit[(∑_(n=1)^∞)(((((-1))^n) (x^(2 n-1)))/((2 n)! Log[2 n])),x->∞]

Derivation:

Inactive[Integrate][(a^(-t)),{t,0,Infinity}]

Assuming[a>1,Activate[%]]

sumOfint=Inactive[Sum][(((-1))^n)((x^(2n-1))/((2n)!))(1/(Log[2n])),{n,1,Infinity}]/. (1/(Log[2n])):>Inactive[Integrate][(((2n))^(-t)),{t,0,Infinity}]

sumOfint2intOfsumRule={Inactive[Sum][c_.Inactive[Integrate][e_,iiter_],siter_]:>Inactive[Integrate][Inactive[Sum][c e,siter],iiter]};

intOfsum=sumOfint/.sumOfint2intOfsumRule

Using

eulerInt=Inactive[Integrate][(u^(t-1))Exp[-2n u],{u,0,Infinity}]

eulerInt==Assuming[n>0&&t>0,Activate[eulerInt]]

res1=intOfsum/.{(2^(-t)) (n^(-t)):>eulerInt/Gamma[t]}/.sumOfint2intOfsumRule

res2=Activate[res1,Sum]

exchangeIntegralsRule={Inactive[Integrate][c_.Inactive[Integrate][e_,iiter_],oiter_]:>Inactive[Integrate][Inactive[Integrate][c e,oiter],iiter]};

res3=res2/.exchangeIntegralsRule

res4=res3/.Inactive[Integrate][e_Times,{t,lb_,ub_}]:>Block[{f=Times@@Cases[e,ti_/;FreeQ[ti,t],{1}]},f Inactive[Integrate][e/f,{t,lb,ub}]]

For large x, the main contribution to the integral comes from large u. The small u part is averaged out to zero by fast oscillations.

Manipulate[Plot[(-1+Cos[(E^(-u)) x])/u,{u,0,10}],{x,{1,10,100,1000}}]

For large u the innermost integral approaches Exp[u], as justified by Euler-McLauren formula

nint[u_Real]:=NIntegrate[(u^(t-1))/Gamma[t],{t,0,Infinity}]

LogPlot[{nint[u],Exp[u]},{u,0,5}]

res5=res4/.Inactive[Integrate][(u^(t-1))/Gamma[t],{t,0,Infinity}]:>Exp[u]

Activate[res5]

Limit[%,x->Infinity]

Utilities

Clear[ChangeIntegrationVariables];
ChangeIntegrationVariables[int:Inactive[Integrate][e_,{t_,lb_,ub_}],eq_Equal,u_]:=Module[{sol,invsol,ne,nlb,nub},
sol=Solve[eq&&lb<t<ub,t,Reals];
If[Length[sol]≠1,Return[int]];
invsol=Solve[eq&&lb<t<ub,u,Reals];
If[Length[invsol]≠1,Return[int]];
sol=First[sol];invsol=First[invsol];
ne=e/.sol;
nlb=Limit[u/.invsol,t->lb,Direction->-1];
nub=Limit[u/.invsol,t->ub,Direction->1];
ne=Refine[ne D[t/.sol,u],nlb<u<nub];
Inactive[Integrate][ne,{u,nlb,nub}]
]

PullConstantOut[int:Inactive[Integrate][e_,{t_,lb_,ub_}],const_]:=const Inactive[Integrate][FullSimplify[e/const,lb<t<ub],{t,lb,ub}]