`Series` with non-analytic terms

Question

I have a function which has a smooth non-analytic term, a simplified version of which is $$ f(\lambda) = \frac{1}{1 + k \lambda + e^{-1/\lambda}}. $$ I want a series expansion of this, in Mathematica, in the sense that the non-analytic term is simply discarded, i.e. in this example the first order result I want would be $$ f(\lambda) = 1 - k \lambda + \mathcal{O}(\lambda^2). $$ However, if I try to use Series as in

Series[1/(1 + k λ + Exp[-1/λ]), {λ, 0, 1}]

I get $$ \frac{e^{1/\lambda}}{1 + e^{1/\lambda}(1 + k\lambda + \mathcal{O}(\lambda^2))}. $$

Is there a way to get what I want? Notes:

• the above is just a simple example, in reality I need to be able to do this for a reasonably general set of functions.
• the non-analyticity is, at least for now, limited to the form $a(\lambda) e^{-b /\lambda}$, where $a$ and $b$ are possibly complicated expressions of various factors, and $a$ is an analytic function of $\lambda$.
• I will further want to specifically extract the coefficients of the series, i.e. $k$ in this example.
• a next step in the computation may also numerically set $\lambda = 0$, so there should be no $1/\lambda$ left in any of the subexpressions.

Update:

The above is just a simplified example, and it seems to be misleading. So let me give a brief overview of what I want to do with this, to clarify the constraints and why it's not feasible to just delete the problematic terms:

The larger goal here is to numerically solve a system of ODE's, where the function $f(\lambda)$ (actually, a set of four functions) appears as a potential. Further, the boundary conditions and further processing of the results depends on the first few coefficients of the Taylor series of $f(\lambda)$. In general, $f(\lambda)$ does not need to be (and is not) a polynomial in $\lambda$.

Note that when I write Taylor series above, I really mean the Taylor series, even if it does not converge to the original function. This is of course at the heart of the problem, since the Taylor series of the example above, when constructed by differentiation, is exactly what I need. However, Mathematica is too smart here, and gives me a series which keeps the non-analytic term. When actually solving the system of ODE's, the non-analytic part must be kept though, since the potentials will be evaluated at large $\lambda$.

Why I cannot simply manually remove these terms when expanding: solving this problem has been encapsulated in a package, which in principle takes a black-box set of potentials $f(\lambda)$, extracts the necessary coefficients, numerically solves the ODE's, and constructs the final observables we want to look at from those results. Currently we are actually tuning the potentials to get certain results, so it is important that trying out a potential does not need manual intervention.

Even further, this package is used also by other people than me, so it is preferable that I won't add any hidden features which might trap an unwary user: for example, if I silently remove exponentials by /. Exp[__/λ] -> 0 before expanding, and the user enters a potential which contains for example $e^{-(a\lambda^2 + b \lambda^3)/\lambda}$ (which will give zero with the above replacement), I will not only produce the wrong result, but I will do so without producing an obvious error! While this example is maybe a bit contrived, and could be solved by Simplifying before the replacement, I'm afraid there could be a more realistic example whose existence I'll realize only after the damage is done.

So, the various classes of solutions I can think of, in order of preference:

1. extract the coefficients of the Taylor series exactly according to the definition, i.e. smooth non-analytic parts get discarded. The more general the set of functions $f(\lambda)$ and forms of non-analyticity this works for, the better.
2. Come up with a replacement rule or some other manipulation applied to the potentials to remove specific non-analyticies, such as $e^{-a/\lambda}$ (where $a$ is a positive real number) in a way which is reasonably robust against accidentally removing needed terms. Bonus if there's a way to detect suspect cases so I can give a warning to the user.
3. Simply let the user enter a separate analytic form for the potential, i.e. have them manually do this work. This I can of course do immediately, but it makes the package more laborious to use.

Having written the above, I'm starting to lean toward just implementing option 1. directly:

taylorSeriesCoefficient[expr_, {x_, x0_, ord_}] := Limit[1/Factorial[ord] D[expr, {x, ord}], x -> x0]

but I'm still open for other suggestions?

Answer 1

It seems to me the main trouble is when we have the form terms + series, in which the terms are not incorporated into the series due to some sort of singularity. The following attempts to determine whether the terms may be dropped because they are O[λ]^n for the appropriate n; it uses Limit to do that.

ClearAll[normalSing];
normalSing::esing = "Unable to find limit of the singularity in ``";
normalSing[ser_, x_ -> x0_, opts : OptionsPattern[Limit]] /; !FreeQ[ser, SeriesData] :=
   Replace[ser,
    {expr_SeriesData :> (           (* already normal *)
       Sow["Normal", "Limit"];             (* optional debugging *)
       expr),
     expr_ :>                       (* try to normalize  series + terms *)
      Module[{res, lim = True},
       res = expr /. Plus[sd_SeriesData, terms__] :> Module[{},
           lim = Limit[Plus[terms]/x^(sd[[5]]/sd[[6]]), x -> x0, opts];
           Sow[lim, "Limit"];              (* optional debugging *)
           sd + lim /; FreeQ[lim, Limit | DirectedInfinity | ComplexInfinity | Interval]
           ];
       If[TrueQ@lim,                       (* optional debugging *)
        Sow["No normalization", "Limit"]];
       (* normalization successful if limit exists or no limit was taken (lim == True) *)
       res /; FreeQ[lim, Limit | DirectedInfinity | ComplexInfinity | Interval]
       ],
     expr_ :> (                     (* failed to normalize *)
       Message[normalSing::esing, expr]; 
       expr)}
    ];

Examples:

examples = {
   {1/(1 + k λ + Exp[-1/λ]), -1},   (* essential sing. term *)
   {1/(1 + k λ + Exp[-1/λ]), 1},    (* essential sing. term *)
   Exp[-2/λ]/(1 + k λ),             (* essential sing. factor *)
   1/(1 + k λ + Sin[1/λ]),          (* Interval[] limit *)
   {1/(1 + k λ + λ^Sqrt[2]), -1},   (* algebraic sing.: keep *)
   {1/(1 + k λ + λ^Sqrt[5]), -1},   (* algebraic sing.: drop *)
   {1/(1 + k λ + Log[λ]), -1},      (* log. sing.: handled by Series[] *)
   Exp[λ]/(1 + k λ)};               (* nonsingular *)

Flatten@{Row[Flatten@{#}, "; "],
      Check[
       {res, lim} = Reap[
         Replace[
          #,
          {{f_, dir_} :> normalSing[Series[f, {λ, 0, 1}], λ -> 0, Direction -> dir],
           f_ :> normalSing[Series[f, {λ, 0, 1}], λ -> 0]}
          ],
         "Limit"];
       {res, lim /. {
          {{s_String}} :> s,
          _?(FreeQ[#, Limit | DirectedInfinity | ComplexInfinity | Interval] &) ->
            "Normalized"}},
       {res, Row[{"Limit failed: ", lim[[1, 1]]}]},
       {normalSing::esing}
       ]} & /@ examples /. Rule -> List //
 Grid[
   Prepend[#, Style[#, Bold] & /@ {"f[λ]; direction  ", "Normalized series", "Status"}],
   Alignment -> Left, Dividers -> {None, All}] &

Answer 2

Perhaps the following replacement will do?

removeNonanalyticExp = Exp[a_] /; MatchQ[Exp[a+O[λ]], Exp[_SeriesData]] -> 0;

Here the non-analytic Exp gets removed:

1/(1+ k λ + Exp[-1/λ]) /. removeNonanalyticExp

1/(1+k λ)

And here the analytic Exp is retained:

1/(1+k λ + Exp[-(λ^2 + λ^3)/λ]) /. removeNonanalyticExp

1/(1+E^((-λ^2-λ^3)/λ)+k λ)

Answer 3

The problem is, that the function has a jump at lambda=0 and two different gradients coming from negative or positive, as you can see in the picture below.

In order not to loose this information, Series gives just the original function, as you see with //Normal//FullSimplify.

Substitute lambda by lambda + d and take the limit d->0 from both sides, to get the right tangente. (unfortunately Assumptions->lambda >0 does not work in Series)

    expr[lambda_] = 1/(1 + k lambda + Exp[-1/lambda])


    ser[d_] = Series[expr[lambda + d], {lambda, 0, 1}] // Normal

    (*   1/(1 + E^(-1/d) + d k) - (
         E^(1/d) (1 + d^2 E^(1/d) k) lambda)/(
         d^2 (1 + E^(1/d) + d E^(1/d) k)^2)     *)

    Limit[ser[d], d -> 0, Direction -> -1]

    (*   1 - k lambda    *)

    Limit[ser[d], d -> 0, Direction -> 1]

    (*   0    *)

    Plot[
         Evaluate[{1 - k lambda, expr[lambda]} /. k -> 3], {lambda, -2, 2}, 
         PlotStyle -> {Blue, Red}, PlotRange -> 5]

Answer 4

It is not quite clear, why do not you discard the exponential term, though you explicitly write that you want it to be discarded:

    Series[expr /. Exp[-1/\[Lambda]] -> 0, {\[Lambda], 0, 1}]

(*  SeriesData[\[Lambda], 0, {1, -k}, 0, 2, 1]  *)

Have fun!