# Mathematica flips the sign of a Maclaurin series

I have the following Mathematica code:

Gdd = {{-((E^(2 AA[z]) L^2 g[z, B])/z^2), 0, 0, 0, 0}, {0, (
E^(2 AA[z]) L^2)/(z^2 g[z, B]), 0, 0, 0}, {0, 0, (
E^(2 AA[z]) L^2)/z^2, 0, 0}, {0, 0, 0, (
E^(B^2 z^2 + 2 AA[z]) L^2)/z^2, B/bb}, {0, 0, 0, -(B/bb), (
E^(B^2 z^2 + 2 AA[z]) L^2)/z^2}};
Gddinverse = {{-((E^(-2 AA[z]) z^2)/(L^2 g[z, B])), 0, 0, 0, 0}, {0, (
E^(-2 AA[z]) z^2 g[z, B])/L^2, 0, 0, 0}, {0, 0, (
E^(-2 AA[z]) z^2)/L^2, 0, 0}, {0, 0, 0, (
bb^2 E^(B^2 z^2 + 2 AA[z]) L^2 z^2)/(
bb^2 E^(2 B^2 z^2 + 4 AA[z]) L^4 + B^2 z^4), -((B bb z^4)/(
bb^2 E^(2 B^2 z^2 + 4 AA[z]) L^4 + B^2 z^4))}, {0, 0, 0, (
B bb z^4)/(bb^2 E^(2 B^2 z^2 + 4 AA[z]) L^4 + B^2 z^4), (
bb^2 E^(B^2 z^2 + 2 AA[z]) L^2 z^2)/(
bb^2 E^(2 B^2 z^2 + 4 AA[z]) L^4 + B^2 z^4)}};
AA[z_] = -a z^2;
g[z_, B_] =
1 - Integrate[z^3 Exp[-B^2 z^2 - 3 AA[z]], {z, 0, z}]/
Integrate[z^3 Exp[-B^2 z^2 - 3 AA[z]], {z, 0, zh}];
ff[z_, B_] = -(1/z^2) 2 E^(2 (B^2 z^2 + AA[z])) L^2 Sqrt[
1 + (B^2 E^(-2 B^2 z^2 - 4 AA[z]) z^4)/(
bb^2 L^4)] (g[z, B] (-2 + 2 B^2 z^2 + 3 z Derivative[1][AA][z]) +
z D[g[z, B], z]);
term[z_, B_] =
Sqrt[-Det[Gdd]] Gddinverse[[2, 2]] Gddinverse[[3, 3]] ff[z, B] //
PowerExpand // Simplify;
V2[z_] = c1 Log[z] (z^4 (1 + a2 z^2)) + (1 + b2 z^2 + b4 z^4);
eq2[z_] =
D[V2[z], {z, 2}] +
D[Log[term[z, B] ], z] D[V2[z],
z] - \[Omega]^2 ((Gddinverse[[1, 1]]/Gddinverse[[2, 2]])) V2[z];
seriesexpneq12 = Series[eq2[z], {z, 0.0, 1}] // Normal

-4 b2 + \[Omega]^2

But when I use 0.0000000000000001 instead of 0.0, I get the answer:

2 b2 + (1. + 1.*10^-32 b2 + 1.*10^-64 b4 - 3.68414*10^-63 c1 -
3.68414*10^-95 a2 c1) \[Omega]^2 + (2.*10^-16 b2 +
4.*10^-48 b4 + (-1.46365*10^-46 -
2.20048*10^-78 a2) c1) (-1.*10^-16 + z) \[Omega]^2

Why does Mathematica flip the sign of the leading order term in this case? Which answer should be considered the correct one?

• I suspect this is showing a bug in Series, and that the first is the incorrect one. Nov 30, 2022 at 16:45

The OP should report that when expanding near 0.0000000000000001 General::munfl will be generated and the numerical result is thus untrustable.

The reason is as follows.

The tools are

expr=eq2[z];
seriesCompare[expr_,order_:1]:=
Simplify@N[
Normal@Series[expr,{z,0,order}]-Chop@Normal@Series[expr,{z,10^-20.,order}]
];

symbolExtract[expr_,wrapper_:Identity,opts:OptionsPattern[]] :=
DeleteDuplicates@Cases[
expr,
$$symbol_Symbol?nonsystemQ:>wrapper[$$symbol],
{0,Infinity},
FilterRules[{opts},Options[Cases]],
];
nonsystemQ[symbol_Symbol] :=
Context@symbol=!="System`";
nonsystemQ[_] = False;

here expr is the expression from OP's code, seriesCompare compares the expansion near 0 and 10^-20, and symbolExtract extracts symbols in an expressions.

To separate the relevant part in expr we take some symbols at some analytic point by trial and error, until seriesCompare generates no-vanishing result,

expr1//symbolExtract
expr1=expr//ReplaceAll[{b4->0,zh->1,a->1(*B->1*),\[Omega]->1,L->1,bb->1,a2->1,c1->1}]//Simplify;

seriesCompare/@List@@expr1//TableForm

To see why the 5-th term generates warnings,

expr1[[5]]//Denominator

so actually there is a z-factor not being cancelled, and if we replace the 5-th term by

expr2=expr1[[5]]//Factor//Cancel//Simplify

expr1[[5]]//seriesCompare
expr2//seriesCompare//Simplify

the results will match,

To understand this difference, let's take

test=Sin[z]/z;
seriesCompare[test,2]
(*0. +1.20893 10^24 z^2*)

the coefficient of z^2 is (with x=10^-20)

e2=Series[test,{z,x,2}]//Normal//Coefficient[#,z,2]&
(*-(Cos[x]/x^2) + Sin[x]/x^3 - Sin[x]/(2 x)*)

this vanishes at $$O(x)$$

Series[e2,{x,0,2}]
(*-(1/6)+x^2/20+O[x]^3*)

but

e2/.x->10^-20//N
(*-1.20893 10^24*)

due to numerical precision.

• Another point should be stressed that Series is actually not Taylor nor Laurent series according to the documentation > Series can construct standard Taylor series, as well as certain expansions involving negative powers, fractional powers, and logarithms. seriesCompare[Log[z]] (*46.0517 - 1.\[CenterDot]10^20 z + Log[z]*) And there are log-terms in expr, although this is irrelevant to expr at order 1. Dec 1, 2022 at 0:32
• What if I take b4->1 instead of b4->0 in the above analysis? Dec 1, 2022 at 6:43
• @codebpr You can try it yourself. The purpose of specifying the parameters is just to simplify the analysis. Intuitively not all choices of parameters can reveal where the problem is located. Dec 3, 2022 at 7:06

Strangely enough both results seem to be correct or at least plausible.

ee = eq2[z] /. \[Omega] -> omega;
seriesexp1 = Normal[Series[ee, {z, 0, 1}]]
seriesexp2 = Normal[Series[ee, {z, 1/100, 1}]];
seriesexp2[[1 ;; 2]]

(* -4 b2 + omega^2

2 b2 + (3 b4)/2500 *)

Now a numeric check.

In[512]:= N[{seriesexp1, seriesexp2} /. {c1 -> 2, a2 -> 3, bb -> 5/11,
B -> 7/13, L -> 11/19, zh -> 13/16, b4 -> 17/3,
a -> 19/5, \[Omega] -> 1, b2 -> 2}, 20]

(* Out[512]= {-8.0000000000000000000 + omega^2, -7.9987469301487475467 +
1.0001999646629277488 omega^2 + (0.24853150921129158616 +
0.039987860268843152426 omega^2) (-0.010000000000000000000 + z)} *)