Mathematica flips the sign of a Maclaurin series

Question

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

It gives the answer as:

-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?

Answer 1

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]],
        Heads->False
    ];
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.

Answer 2

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)} *)