4
$\begingroup$

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?

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

2 Answers 2

4
$\begingroup$

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

enter image description here

To see why the 5-th term generates warnings,

expr1[[5]]//Denominator

enter image description here

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,

enter image description here

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.

$\endgroup$
3
  • 1
    $\begingroup$ 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. $\endgroup$
    – Lacia
    Dec 1, 2022 at 0:32
  • $\begingroup$ What if I take b4->1 instead of b4->0 in the above analysis? $\endgroup$
    – codebpr
    Dec 1, 2022 at 6:43
  • $\begingroup$ @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. $\endgroup$
    – Lacia
    Dec 3, 2022 at 7:06
1
$\begingroup$

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)} *)
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.