Maxima -- GPL CAS based on DOE-MACSYMA / Bugs / #3922 Zeroth order Taylor expansion is wrong

Stavros Macrakis - 2022-01-19

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.

If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
- boud - 2022-01-19
  
  Nice work in finding a much more minimal working example. :) We'll see if the solution fixes my original expression...
  
  If you would like to refer to this comment somewhere else in this project, copy and paste the following link:

Stavros Macrakis - 2022-01-19

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.

If you would like to refer to this comment somewhere else in this project, copy and paste the following link:

Robert Dodier - 2022-01-20

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.

If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
- Stavros Macrakis - 2022-01-20
  
  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.
  
  If you would like to refer to this comment somewhere else in this project, copy and paste the following link:

boud - 2022-01-20

@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
If you would like to refer to this comment somewhere else in this project, copy and paste the following link:

boud - 2022-01-22

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.
If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
- Stavros Macrakis - 2022-01-22
  
  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.
  
  If you would like to refer to this comment somewhere else in this project, copy and paste the following link:

boud - 2022-01-22

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".

If you would like to refer to this comment somewhere else in this project, copy and paste the following link:

Stavros Macrakis - 2022-01-22

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.

Last edit: Stavros Macrakis 2022-01-22

If you would like to refer to this comment somewhere else in this project, copy and paste the following link:

boud - 2022-01-25

The context of this bug is our paper https://arxiv.org/abs/2201.09102 with the Maxima scripts at https://codeberg.org/boud/topoaccel . The key formula of interest is (21), summed over locations in S^3 that have been converted to an R^4 embedding space, giving (22), where in the code we write rr instead of xi. If we define the main formula

U: 1/rr - cot(rr)*(1-rr/%pi); /* remove the singularity */ taylor(U,rr,0,5); /* accepted by maxima */ dUdrr: diff(U,rr,1); /* maxima finds an expression for the derivative */ dUdrr, rr:0; /* undefined - 0 to neg exponent */ limit(dUdrr, rr,0); /* the limit at 0 exists */

then we do not have the zeroth order bug in this case:

taylor(U,rr,0,1); taylor(U,rr,0,0);

Both give 1/pi as the zeroth order. The bug occurs when we do a sum - the other parts of the sum should have singularities, but not nearrr=0. The particular instance of the bug can be seen 'in the wild' in (old) version be02b0f4:
https://codeberg.org/boud/topoaccel/commit/be02b0f47c87b1d927e4b6570bb164383cb92b3f

At line 120 here:
https://codeberg.org/boud/topoaccel/src/commit/be02b0f47c87b1d927e4b6570bb164383cb92b3f/max_Mass_S3_8.m#L120

add a second line below, requesting only the zeroth order instead of the first order:

print("Phi_0*V^1/3 for the integral convention : ",bfloat(taylor(potFull*FD_vol^(1/3),rr,0,0)));

and run the full script maxima < main_max_Mass.m > mylog.
The output will then be

Phi_0*V^1/3 for the integral convention : 4.50427948343915b-1 rr + 2.580762041484299b0 Phi_0*V^1/3 for the integral convention : 2.150635034570249b0

Here we do have the bug: 6/(4*%pi)^(1/3) approx 2.581 versus 5/(4*%pi)^(1/3) approx 2.151.
If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
- Stavros Macrakis - 2022-01-31
  
  As I was working through your example, I noticed that there's a bug in limit.
  See https://sourceforge.net/p/maxima/bugs/3932/, which is based on your dUdrr.
  Maxima is returning 0 for that limit, but 1/3 for limit(trigsimp(dUdrr)...).
  
  If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
- Stavros Macrakis - 2022-01-31
  
  Thanks for working with us to help track down the problem.
  Could you please supply a minimal stand-alone example that shows the problem?
  What are the steps to reproduce the problem starting from a fresh Maxima?
  
  If you would like to refer to this comment somewhere else in this project, copy and paste the following link:

Zeroth order Taylor expansion is wrong

Computer Algebra System written in Common Lisp

Group

Searches

Help

#3922 Zeroth order Taylor expansion is wrong

Discussion