Asymptotic Expansions

Adam Strzebonski and Devendra Kapadia, 
Wolfram Research, Inc.

Plan for the Talk

We will introduce three new functions for asymptotic expansions:

AsymptoticSolve - asymptotic solutions of algebraic equations

AsymptoticIntegrate - asymptotic expansions of integrals

AsymptoticDSolve - asymptotic solutions of differential equations

Asymptotic Solutions of Algebraic Equations

Code

AsymptoticSolve[f_,args___]:=With[{ans=Asymptotics`SeriesSolutions[If[ListQ[f],Subtract@@@f,Subtract@@f],args]},Normal[ans,SeriesData]/;FreeQ[ans,Asymptotics`SeriesSolutions|$Failed]]

AsymptoticSolve

AsymptoticSolve[eqn,{y,(Subscript[y, 0])},{x,(Subscript[x, 0]),n}]

computes series solutions of eqn in a neighbourhood of {(Subscript[x, 0]),(Subscript[y, 0])} up to the order (((x-(Subscript[x, 0])))^n).

Non-singular Points

Find the series solution in a neighbourhood of a non-singular point of a curve.

f=(E^(x+y))-Sin[x y]-1;
sols=AsymptoticSolve[f==0,{y,0}, {x, 0,7}]

Compare the asymptotic solution and a numeric solution.

nsol[xx_?NumberQ,i_]:=Block[{x=xx},y/.FindRoot[f,{y,sols[[i]]}]]

Plot[{nsol[x,1],sols[[1]]}, {x,-0.4,0.4}, PlotLegends ->{"numeric solution","asymptotic solution"}]

Singular Points

Find the Puiseux series solutions in a neighbourhood of a singular point of a curve.

f=(y^7) Cos[x]+((Sin[y])^5) Sin[x]-Tan[(y^3)] x-((Tan[x])^5);
AsymptoticSolve[f==0,{y,0},{x,0,4}]

Find just the solutions with real coefficients.

sols=AsymptoticSolve[f==0,{y,0},{x,0,4},Reals]

Compare the asymptotic solutions and numeric solutions.

nsol[xx_?NumberQ,i_]:=Block[{x=xx},y/.FindRoot[f,{y,sols[[i]]}]]

Plot[{sols[[1]],nsol[x,1],sols[[2]],nsol[x,2],sols[[3]],nsol[x,3]}, {x,0,0.5}, PlotLegends ->{"numeric solution 1","asymptotic solution 1","numeric solution 2","asymptotic solution 2","numeric solution 3","asymptotic solution 3"}]

Algebraic curves in n-dimensions

Find the Puiseux series solutions in a neighbourhood of a singular point of an algebraic curve in 3-dimensional space.

{f,g}={x^2-y^2+z x y+z^2 x+z,z^2 x^2+z y^2+x y+z y+z^2};
sols=AsymptoticSolve[{f==0,g==0},{{x,y},{0,0}},{z,0,4}]

Compare the asymptotic solutions and numeric solutions.

rsols={Append[sols[[1]],z],Append[sols[[2]]/.z->-z,-z],Append[sols[[3]]/.z->-z,-z],Append[sols[[4]],z]};
nsols3D[zz_?NumberQ,i_]:=Block[{z=zz},{x,y,rsols[[i,3]]}/.FindRoot[{f,g},{{x,rsols[[i,1]]},{y,rsols[[i,2]]}}]]

ParametricPlot3D[{rsols[[1]],nsols3D[z,1],rsols[[2]],nsols3D[-z,1],rsols[[3]],nsols3D[-z,3],rsols[[4]],nsols3D[z,4]},{z,0,0.3},PlotLegends ->{"numeric solution 1","asymptotic solution 1","numeric solution 2","asymptotic solution 2","numeric solution 3","asymptotic solution 3","numeric solution 4","asymptotic solution 4"}]

Possible extensions

• analytic curves in n-dimensions

• k-dimensional smooth surfaces in n-dimensions (power series; implicit function theorem)

• singular surfaces in n-dimensions (piecewise fractional power series)

• asymptotic expansions at infinity

Asymptotic Expansions of Integrals

AsymptoticIntegrate

AsymptoticIntegrate[expr, {t, a,b}, {ϵ, p, n}]

computes the first n terms in the asymptotic expansion of

Integrate[expr, {t, a, b}]

with respect to the parameter ϵ, around p.

Series Expansion

Obtain an asymptotic expansion for a complete elliptic integral. This relies on a series expansion of the integrand:

AsymptoticIntegrate[(1/(Sqrt[1-m ((Sin[θ])^2)])),{θ,0,(π/2)},{m,0,5}]

Compare with the series expansion of the exact result:

Integrate[(1/(Sqrt[1-m ((Sin[θ])^2)])),{θ,0,(π/2)},Assumptions->0<m<1]

Normal[Series[%,{m,0,5}]]

Integration by Parts

Obtain the first six terms in the asymptotic expansion for the integral of an exponential function from x to Infinity:

f[t_]:=((E^(-t))/t)

Plot[f[t],{t,30,32.5},PlotRange->{0,(1/(10^14.5))},Filling->Axis]

The result is obtained using repeated integration by parts:

int=AsymptoticIntegrate[f[t],{t,x,∞},{x,∞,6}]

Compare with a numerical approximation:

NIntegrate[f[t],{t,30,Infinity}]

int/. x->30.

Exponential and Oscillatory Integrals

Laplace Method

Consider an function involving an exponential kernel that has a maximum at x = 1:

f[x_]:=(E^(-ω ((x^2)-2 x))) (((1+x))^(5/2))

Plot the function for a specific value of ω :

Plot[f[x]/. {ω->3},{x,0,10},PlotRange->All,Filling->Axis,FillingStyle->Yellow]

Obtain the leading term in the expansion of the integral over {0,Infinity} using Laplace’s method:

AsymptoticIntegrate[f[x],{x,0,∞},{ω,∞,1}]

Complare with the numerical result:

%/. {ω->30.}

NIntegrate[Exp[-30 ((x^2)-2 x)] (((1+x))^(5/2)),{x,0,∞}]

Obtain a better approximation by computing an extra term:

AsymptoticIntegrate[f[x],{x,0,∞},{ω,∞,2}]

%/. {ω->30.}

Method of Stationary Phase

Compute the leading term approximation for the integral of an oscillatory function over {0, Pi}:

f[x_]:=(E^(I ω Sin[t]))

Plot the real part of the function for a particular value of the parameter:

Plot[Re[f[x]/. {ω->50}],{t,0,π},Filling->Axis,FillingStyle->Yellow]

The approximation is computed using the method of stationary phase:

int=AsymptoticIntegrate[f[t],{t,0,π},{ω,∞,1}]

Complare with the numerical approximation:

NIntegrate[Exp[I 500 Sin[t]],{t,0,π},MinRecursion->20,MaxRecursion->20]

int/. ω->500.

Series Solutions of Differential Equations

AsymptoticDSolve

AsymptoticDSolve[deqn, y[x], {x, p, n}]

computes a series solution of order n for deqn around p.

Ordinary Points

The following linear first-order ODE has an ordinary point x=0 :

ode={y'[x]-x y[x]==0,y[0]==1};

Compute a Taylor series solution at the ordinary point :

AsymptoticDSolve[ode,y[x],{x,0,10}]

Plot the successive approximations given by the solution:

Plot[Evaluate[Accumulate[List@@%]],{x,0,3},PlotLegends->"Expressions"]

Regular Singular Points

The following linear second - order ODE has a regular singular point x = 0 :

ode=2 (x^2) y''[x]+3 x y'[x]-((x^2)+1) y[x]==0;

Compute a Frobenius series solution at the regular singular point:

AsymptoticDSolve[ode,y[x],{x,0,5}]

Irregular Singular Points

The Airy equation has an irregular singular point x = Infinity:

airyode=y''[x]-x y[x]==0;

Compute an asymptotic series solution at the irregular singular point:

AsymptoticDSolve[airyode,y[x],{x,∞,3}]

Compare with the expansions for the Airy functions at Infinity:

s1=Normal[Series[AiryAi[x],{x,∞,4}]]

s2=Normal[Series[AiryBi[x],{x,∞,4}]]

Visualize the approximations:

Plot[Evaluate[{AiryAi[x],AiryBi[x],s1,s2}],{x,-3,3},PlotLegends->{AiryAi[x],AiryBi[x],"s1","s2"},
 PlotStyle->Thickness[0.008]]

Nonlinear ODE

The following is a nonlinear first-order ODE:

eqns={3 ((y'[x])^2)+4 x y'[x]-y[x]+(x^2)==0,y[0]==1};

Compute a Taylor series approximation at x = 0:

Quiet[AsymptoticDSolve[eqns,y[x],{x,0,32}]]

Verify that the polynomial approximation is exact:

Simplify[eqns/. {y->Function[{x},Evaluate[%]]}]

System of ODEs

Set up a system of linear first-order ODEs with a regular singular point at t = 0:

r=(1/6) {{5,-6},{4,-6}};

Eigenvalues[r]

a={2,Sin[t]};

x[t_]={u[t],v[t]};

eqns=Thread[x'[t]==((r/t)+a).x[t]]

Compute a Frobenius series solution for the problem:

AsymptoticDSolve[eqns,x[t],{t,0,2}]

Perturbation Expansions

AsymptoticDSolve - Perturbation Solutions

AsymptoticDSolve[deqn, y, x, {ϵ, p, n}]

computes a perturbation solution of order n for deqn with respect to ϵ around p.

Duffing Equation

Consider the Duffing equation:

eqn={u''[x]+u[x]+ϵ ((u[x])^3)==0,u[0]==3,u'[0]==0};

Find a first-order Poincaré -Lindstedt perturbation expansion for the problem:

psol=AsymptoticDSolve[eqn,u[x],x,{ϵ,0,1}]

Plot the approximate solution:

Plot[psol/. {ϵ->(1/800)},{x,0,20},PlotPoints->100,WorkingPrecision->20]

Compare with the exact solution which is expressed in terms of Jacobian elliptic funtions:

dsol=DSolveValue[eqn/. {ϵ->(1/800)},u[x],x]

Plot[dsol // Evaluate, {x, 0, 20}, PlotPoints -> 100, WorkingPrecision -> 20]

Boundary Layer

The following is an example of a boundary layer (singular perturbation) problem:

eqn={ϵ y''[x]+2 y'[x]+y[x]==0,y[0]==0,y[1]==(1/2)};

Find a WKBJ approximation for the problem:

psol=AsymptoticDSolve[eqn,y[x],x,{ϵ,0,1}]

Compute numerical values for the solution:

Table[%/. {ϵ->(1/30)},{x,0,0.5,0.1}]

Compare with the exact solution:

dsol=DSolveValue[eqn,y[x],x]

N[Table[%/. {ϵ->(1/30)},{x,0,0.5,0.1}]]

Plot the exact and approximate solutions:

Plot[Evaluate[{psol,dsol+0.01}/. {ϵ->(1/30)}],{x,0,1},PlotStyle->{Red,Blue}]

Conclusion

Asymptotic expansions offer a powerful alternative to symbolic and numeric methods.

Our current focus is on expansions of algebraic/differential equations and integrals.

We also plan to develop Asymptotic solvers for sums and difference equations.

Other important objectives are to improve Series and develop general tools for Asymptotics.

Initialization

Code

<<"C:\\Users\\adams\\Desktop\\AsymptoticExpansions\\AsymptoticDSolve.wl"

<<"C:\\Users\\adams\\Desktop\\AsymptoticExpansions\\AsymptoticIntegrate.wl"