Zeroth order Taylor expansion is wrong

[rtoy]

Imported from SourceForge on 2024-07-04 18:08:41 Created by boud1 on 2022-01-19 19:46:06 Original: https://sourceforge.net/p/maxima/bugs/3922

---

DESCRIPTION:

The zeroth order Taylor expansion of a scientifically interesting expression is incorrectly evaluated if only the zeroth order is requested. The correct value requires evaluating to the first order and then ignoring the first order.

SYSTEM:

Debian GNU/Linux 11.2 (bullseye)
ii  maxima      5.44.0-3    amd64        Computer algebra system -- base system

Maxima version: "5.44.0"
Maxima build date: "2021-04-24 14:52:58"
Host type: "x86_64-pc-linux-gnu"
Lisp implementation type: "GNU Common Lisp (GCL)"
Lisp implementation version: "GCL 2.6.12"
User dir: "~/.maxima"
Temp dir: "/tmp"
Object dir: "~/.maxima/binary/5_44_0/gcl/GCL_2_6_12"
Frontend: false

The bug also occurs for

Maxima version: "5.38.1"
Lisp implementation version: "GCL 2.6.12"

REPRODUCE the bug:

/* The expression Q, depending on rr, x and y, yields a bug in
'taylor'. */

Q: (%pi^(2/3)*((-(sin(rr)*sqrt((-y^2)-x^2+1)
        *(1-acos(sin(rr)*sqrt((-y^2)-x^2+1))/%pi))
      /sqrt(1-sin(rr)^2*((-y^2)-x^2+1)))
    +(sin(rr)*sqrt((-y^2)-x^2+1)*(1-(%pi-acos(sin(rr)*sqrt((-y^2)-x^2+1)))/%pi))
    /sqrt(1-sin(rr)^2*((-y^2)-x^2+1))
    -(sin(rr)*y*(1-acos(sin(rr)*y)/%pi))/sqrt(1-sin(rr)^2*y^2)
    +(sin(rr)*y*(1-(%pi-acos(sin(rr)*y))/%pi))/sqrt(1-sin(rr)^2*y^2)
    -(sin(rr)*x*(1-acos(sin(rr)*x)/%pi))/sqrt(1-sin(rr)^2*x^2)
    +(sin(rr)*x*(1-(%pi-acos(sin(rr)*x))/%pi))/sqrt(1-sin(rr)^2*x^2)
    -(cos(rr)*(1-acos(cos(rr))/%pi))/sqrt(1-cos(rr)^2)
    +(cos(rr)*(1-(%pi-acos(cos(rr)))/%pi))/sqrt(1-cos(rr)^2)+1/rr+4/%pi))/4^(1/3);

/* The 0th order coefficient of the Taylor expansion around rr set to
   0 should be independent of whether the expansion is made to 0th or
   first order.
*/

/* This should give 0, but instead gives  (4 *%pi)^(-1/3). */
coeff(taylor(Q,rr,0,1),rr,0) - coeff(taylor(Q,rr,0,0),rr,0);

/* This should give 1, but instead gives  6/5. */
coeff(taylor(Q,rr,0,1),rr,0) / coeff(taylor(Q,rr,0,0),rr,0);

/* Is the greater expansion correct? Check with a few arbitrary
   values.  This shows that coeff(taylor(Q,rr,0,1),rr,0) is very
   likely to be correct and coeff(taylor(Qtest,rr,0,0),rr,0) wrong.
*/

coeff(taylor(Q,rr,0,1),rr,0), bfloat;
coeff(taylor(Q,rr,0,0),rr,0), bfloat;
Qtest: Q, rr: 1e-3, x:0, y:0;
taylor(Qtest,rr,0,5), bfloat;

DISCUSSION:

For the case of interest, rr is positive, but including assume(rr > 0) at the beginning doesn't solve the bug. Restricting to assume(rr < %pi/2) would also be acceptable, but doesn't solve the bug either. Together these would avoid any issues of cos and sin domains and invertibility.

[rtoy]

Imported from SourceForge on 2024-07-04 18:08:42 Created by macrakis on 2022-01-19 20:13:24 Original: https://sourceforge.net/p/maxima/bugs/3922/#7bf7

---

Thanks for the bug report!

A simpler case showing the same problem is R:(rr*cos(rr))/sqrt(1-cos(rr)^2). taylor(R,rr,0,0) => 0+... while taylor(R,rr,0,1) => 1+....

The problem is probably related to the square root. That's not an excuse for bad results, just a hint about where we might look to debug the problem.

[rtoy]

Imported from SourceForge on 2024-07-04 18:08:46 Created by macrakis on 2022-01-19 20:22:09 Original: https://sourceforge.net/p/maxima/bugs/3922/#d7f0

---

Interestingly, tlimit(R,rr,0) correctly yields 1, though it does ask the sign of sin(rr) -- result is the same regardless of the answer. Sadly, limit gets this wrong -- I will file a separate bug report.

[rtoy]

Imported from SourceForge on 2024-07-04 18:08:50 Created by boud1 on 2022-01-19 21:49:36 Original: https://sourceforge.net/p/maxima/bugs/3922/#7bf7/8336

---

Nice work in finding a much more minimal working example. :) We'll see if the solution fixes my original expression...

[rtoy]

Imported from SourceForge on 2024-07-04 18:08:53 Created by robert_dodier on 2022-01-20 06:52:41 Original: https://sourceforge.net/p/maxima/bugs/3922/#010c

---

Well, it looks like Q is discontinuous at rr = 0., and it's not removable. Is a Taylor series defined for such a function? Even if it is, it seems likely that Maxima's taylor function is going to have trouble with it.

I guess we could try to make headway by only looking at one-sided limits for the function and its derivatives -- I don't know how much trouble it would be, or whether it would fix the bug observed here.

At this point I would say the bug is that taylor has not detected that it cannot construct a valid Taylor series for this problem.

[rtoy]

Imported from SourceForge on 2024-07-04 18:08:57 Created by boud1 on 2022-01-20 11:01:45 Original: https://sourceforge.net/p/maxima/bugs/3922/#36b4

---

@robert_dodier I missed that: you're right: Q is discontinuous at rr = 0. Under the Wikipedia definition [1], what counts is whether Q is infinitely differentiable or not at rr=0. It looks to me like the limit as rr \rightarrow 0^+ is around 2.58 and as rr \rightarrow 0^- is -\infty, so I don't see how the lefthand and righthand first derivatives could be equal. So I agree that strictly speaking, the error in taylor is that it did not detect the non-existence of the Taylor series.

I forgot to say that the domain of interest is the positive reals, for rr. x and y are reals of max abs value 1. Since one-sided real derivatives exist (in general), it seems to me that one-sided real Taylor series should also exist.

The following

assume(sin(rr) > 0);
/* Q: ... as above */
tlimit(Q,rr,0,plus);

gives me 6/ ( 4^(1/3) * %pi^(1/3) ), which is the expected result for a real, positive, small (rr \ll \pi) result. This solves my personal, immediate problem, though not the general issue.

Should I file a separate wishlist item for extending taylor to optionally (opt-in) allow real, one-sided expansions?

As for the more general issue, I guess that the behaviour of taylor should depend on numer_pbranch for cases like this one.

[1] https://en.wikipedia.org/wiki/Taylor_series#Definition

[rtoy]

Imported from SourceForge on 2024-07-04 18:09:00 Created by macrakis on 2022-01-20 17:54:52 Original: https://sourceforge.net/p/maxima/bugs/3922/#010c/8ad6

---

R and hence Q are only discontinuous at rr=0 if you interpret sqrt as the "positive square root". But "positive square root" is not an analytic function, and taylor, limit, etc. work on analytic functions.

The kernel of the problem is R:(rr*cos(rr))/sqrt(1-cos(rr)^2). Numerical evaluation gives a negative result for rr=-.01 and a positive result for rr=.01 because it interprets sqrt as the positive square root. But as an analytic function of rr, R should stay on a consistent branch across rr=0. The simplification sqrt(x^2) => abs(x) is incorrect if you are working with analytical functions.

[rtoy]

Imported from SourceForge on 2024-07-04 18:09:04 Created by boud1 on 2022-01-22 16:51:05 Original: https://sourceforge.net/p/maxima/bugs/3922/#021a

---

The simplified expression can be made analytic across zero over the reals, it seems to me, using signum to compensate for not allowing analytic continuation of sqrt through zero.

T: signum(rr)*(rr*cos(rr))/sqrt(1-cos(rr)^2);
assume(rr < 0, sin(rr) < 0);
Tp1: T, rr:-1e-3;
Tp2: T, rr:-1e-4;
printf(true,"rr < 0 gives ~a ~a~%", Tp1, Tp2)$
forget(rr < 0, sin(rr) < 0); assume(rr > 0, sin(rr) > 0);
Tp1: T, rr:1e-3;
Tp2: T, rr:1e-4;
printf(true,"rr > 0 gives ~a ~a~%", Tp1, Tp2)$

Then we have the taylor expansions to 0th, 1st and 5th order:

taylor(T,rr,0,0);
taylor(T,rr,0,1); 
taylor(T,rr,0,5);

So now we have T is analytic over the reals and taylor(T,rr,0,0) incorrectly gives 0 + ... instead of 1 + .... So this does look like a bug to me.

In terms of successive steps of Taylor series calculations, going from (r^2 + ...)^(1/2) to signum(r) r + ... might help to handle this - if the domain is restricted to the reals.

[rtoy]

Imported from SourceForge on 2024-07-04 18:09:08 Created by macrakis on 2022-01-22 18:27:58 Original: https://sourceforge.net/p/maxima/bugs/3922/#021a/74cb

---

I'm not sure what your point is here. As I said before, taylor ignores assumptions like rr>0, and the signum simplifies away in the presence of such an assumption.

As I said before, numerical evaluation is mathematically naive. It treats sqrt(1-cos(rr)^2) as though it were abs(sin(rr)) -- using the positive square root. Of course, that is not differentiable at rr=0, and therefore doesn't have a Taylor expansion there.

As an analytical function, it is either sin(rr) or -sin(rr) at rr=0, depending which branch you're on.

[rtoy]

Imported from SourceForge on 2024-07-04 18:09:12 Created by boud1 on 2022-01-22 20:50:13 Original: https://sourceforge.net/p/maxima/bugs/3922/#8a3d

---

I agree that numerical evaluation differs from normal maths (I like the URLs that you provided in some other maxima comments). Regarding what taylor does and doesn't do, I'm happy to create a separate wishlist item for an option to taylor which allows it to be evaluated in restricted domains (e.g. just the reals, or just the positive reals). Would taking into account the signs of variables require a fundamental rewrite of the whole of hayat.lisp?

I was thinking of defining the function g(x) := signum(x) * sqrt(x). Mathematically, this is infinitely differentiable at zero. However, the first and all greater derivatives are singular at zero. So even restricted to the reals, the Taylor series doesn't exist, it seems to me.

I guess this is the reason for taylor ignoring signum: signum is discontinuous at zero, the discontinuity is not removable, and the one-sided derivatives at zero are unequal.

On the other hand, mathematically, is the functionh(x) := signum(x) * abs(sin(x)) mathematically identical to sin(x)? I expect that there would be a mathematical problem if we want the chain rule of differentiation to hold, because the two individual functions are individually non-differentiable at zero, unless we take one-sided derivatives only. So the space of functions and some of the usual rules would have to allow special cases such as this one. Extending to https://en.wikipedia.org/wiki/Distribution_%28mathematics%29 would seem like overkill to me, and I expect would be difficult to handle in Maxima, though people who know Maxima better than me could judge that more realistically.

I wonder if the following discussion could be useful for a simpler approach than distributions: https://en.wikipedia.org/wiki/Exponentiation#Failure_of_power_and_logarithm_identities This includes warnings such as

"one has implicitly supposed that log ⁡exp ⁡z = z for complex values of z, which is wrong, as the complex logarithm is multivalued".

[rtoy]

Imported from SourceForge on 2024-07-04 18:09:15 Created by macrakis on 2022-01-22 21:12:12 Original: https://sourceforge.net/p/maxima/bugs/3922/#aeaf

---

You are welcome to work on one-sided derivatives if you think that's useful. In many cases, you can simply remove the one-sidedness and go from there. The right derivative of abs(x) at x=0 is the same as the derivative of x for x>=0, and the same as the derivative of -x for x<0.

Currently, Maxima gives the derivative of abs(x) as abs(x)/x, which is the same as the right and the left derivatives everywhere except for x=0, where it is undefined.