search

sagemath / sage

Main repository of SageMath
https://www.sagemath.org
Other
1.45k stars 482 forks source link

Integration returns nonsensical answer #36662

Open beew opened 1 year ago

beew commented 1 year ago

Steps To Reproduce

a, b = var('a b')

integrate(cos(x)/(1 - (b^2)*(sin(x - a))^2)^(3/2), x)

Expected Behavior

It should return 

((2 - b^2)*sin(x) + b^2*sin(2*a - x) )/(2*(1 -b^2)*sqrt(1 - (b^2)*(sin(a - x))^2))

or something equivalent. This can be verified either by differentiation or integrating manually.

Mathematica gives

(-b^2 Sin[2 a - x] + (-2 + b^2) Sin[x])/((-1 + b^2) Sqrt[4 - 2 b^2 + 2 b^2 Cos[2 (a - x)]])

Which is easily seen to be equivalent

Actual Behavior

sage gave

2*sqrt(-4*b^2*tan(1/2*a)^4*tan(1/2*x)^2 + 8*b^2*tan(1/2*a)^3*tan(1/2*x)^3 - 4*b^2*tan(1/2*a)^2*tan(1/2*x)^4 + tan(1/2*a)^4*tan(1/2*x)^4 - 8*b^2*tan(1/2*a)^3*tan(1/2*x) + 16*b^2*tan(1/2*a)^2*tan(1/2*x)^2 + 2*tan(1/2*a)^4*tan(1/2*x)^2 - 8*b^2*tan(1/2*a)*tan(1/2*x)^3 + 2*tan(1/2*a)^2*tan(1/2*x)^4 - 4*b^2*tan(1/2*a)^2 + tan(1/2*a)^4 + 8*b^2*tan(1/2*a)*tan(1/2*x) - 4*b^2*tan(1/2*x)^2 + 4*tan(1/2*a)^2*tan(1/2*x)^2 + tan(1/2*x)^4 + 2*tan(1/2*a)^2 + 2*tan(1/2*x)^2 + 1)*(((b^8*tan(1/2*a)^25 + 10*b^8*tan(1/2*a)^23 - b^6*tan(1/2*a)^25 + 44*b^8*tan(1/2*a)^21 - 10*b^6*tan(1/2*a)^23 + 110*b^8*tan(1/2*a)^19 - 44*b^6*tan(1/2*a)^21 + 165*b^8*tan(1/2*a)^17 - 110*b^6*tan(1/2*a)^19 + 132*b^8*tan(1/2*a)^15 - 165*b^6*tan(1/2*a)^17 - 132*b^6*tan(1/2*a)^15 - 132*b^8*tan(1/2*a)^11 - 165*b^8*tan(1/2*a)^9 + 132*b^6*tan(1/2*a)^11 - 110*b^8*tan(1/2*a)^7 + 165*b^6*tan(1/2*a)^9 - 44*b^8*tan(1/2*a)^5 + 110*b^6*tan(1/2*a)^7 - 10*b^8*tan(1/2*a)^3 + 44*b^6*tan(1/2*a)^5 - b^8*tan(1/2*a) + 10*b^6*tan(1/2*a)^3 + b^6*tan(1/2*a))*tan(1/2*x)/(b^8*tan(1/2*a)^24 + 12*b^8*tan(1/2*a)^22 - 2*b^6*tan(1/2*a)^24 + 66*b^8*tan(1/2*a)^20 - 24*b^6*tan(1/2*a)^22 + b^4*tan(1/2*a)^24 + 220*b^8*tan(1/2*a)^18 - 132*b^6*tan(1/2*a)^20 + 12*b^4*tan(1/2*a)^22 + 495*b^8*tan(1/2*a)^16 - 440*b^6*tan(1/2*a)^18 + 66*b^4*tan(1/2*a)^20 + 792*b^8*tan(1/2*a)^14 - 990*b^6*tan(1/2*a)^16 + 220*b^4*tan(1/2*a)^18 + 924*b^8*tan(1/2*a)^12 - 1584*b^6*tan(1/2*a)^14 + 495*b^4*tan(1/2*a)^16 + 792*b^8*tan(1/2*a)^10 - 1848*b^6*tan(1/2*a)^12 + 792*b^4*tan(1/2*a)^14 + 495*b^8*tan(1/2*a)^8 - 1584*b^6*tan(1/2*a)^10 + 924*b^4*tan(1/2*a)^12 + 220*b^8*tan(1/2*a)^6 - 990*b^6*tan(1/2*a)^8 + 792*b^4*tan(1/2*a)^10 + 66*b^8*tan(1/2*a)^4 - 440*b^6*tan(1/2*a)^6 + 495*b^4*tan(1/2*a)^8 + 12*b^8*tan(1/2*a)^2 - 132*b^6*tan(1/2*a)^4 + 220*b^4*tan(1/2*a)^6 + b^8 - 24*b^6*tan(1/2*a)^2 + 66*b^4*tan(1/2*a)^4 - 2*b^6 + 12*b^4*tan(1/2*a)^2 + b^4) - (b^8*tan(1/2*a)^26 + 9*b^8*tan(1/2*a)^24 - 2*b^6*tan(1/2*a)^26 + 34*b^8*tan(1/2*a)^22 - 22*b^6*tan(1/2*a)^24 + b^4*tan(1/2*a)^26 + 66*b^8*tan(1/2*a)^20 - 112*b^6*tan(1/2*a)^22 + 13*b^4*tan(1/2*a)^24 + 55*b^8*tan(1/2*a)^18 - 352*b^6*tan(1/2*a)^20 + 78*b^4*tan(1/2*a)^22 - 33*b^8*tan(1/2*a)^16 - 770*b^6*tan(1/2*a)^18 + 286*b^4*tan(1/2*a)^20 - 132*b^8*tan(1/2*a)^14 - 1254*b^6*tan(1/2*a)^16 + 715*b^4*tan(1/2*a)^18 - 132*b^8*tan(1/2*a)^12 - 1584*b^6*tan(1/2*a)^14 + 1287*b^4*tan(1/2*a)^16 - 33*b^8*tan(1/2*a)^10 - 1584*b^6*tan(1/2*a)^12 + 1716*b^4*tan(1/2*a)^14 + 55*b^8*tan(1/2*a)^8 - 1254*b^6*tan(1/2*a)^10 + 1716*b^4*tan(1/2*a)^12 + 66*b^8*tan(1/2*a)^6 - 770*b^6*tan(1/2*a)^8 + 1287*b^4*tan(1/2*a)^10 + 34*b^8*tan(1/2*a)^4 - 352*b^6*tan(1/2*a)^6 + 715*b^4*tan(1/2*a)^8 + 9*b^8*tan(1/2*a)^2 - 112*b^6*tan(1/2*a)^4 + 286*b^4*tan(1/2*a)^6 + b^8 - 22*b^6*tan(1/2*a)^2 + 78*b^4*tan(1/2*a)^4 - 2*b^6 + 13*b^4*tan(1/2*a)^2 + b^4)/(b^8*tan(1/2*a)^24 + 12*b^8*tan(1/2*a)^22 - 2*b^6*tan(1/2*a)^24 + 66*b^8*tan(1/2*a)^20 - 24*b^6*tan(1/2*a)^22 + b^4*tan(1/2*a)^24 + 220*b^8*tan(1/2*a)^18 - 132*b^6*tan(1/2*a)^20 + 12*b^4*tan(1/2*a)^22 + 495*b^8*tan(1/2*a)^16 - 440*b^6*tan(1/2*a)^18 + 66*b^4*tan(1/2*a)^20 + 792*b^8*tan(1/2*a)^14 - 990*b^6*tan(1/2*a)^16 + 220*b^4*tan(1/2*a)^18 + 924*b^8*tan(1/2*a)^12 - 1584*b^6*tan(1/2*a)^14 + 495*b^4*tan(1/2*a)^16 + 792*b^8*tan(1/2*a)^10 - 1848*b^6*tan(1/2*a)^12 + 792*b^4*tan(1/2*a)^14 + 495*b^8*tan(1/2*a)^8 - 1584*b^6*tan(1/2*a)^10 + 924*b^4*tan(1/2*a)^12 + 220*b^8*tan(1/2*a)^6 - 990*b^6*tan(1/2*a)^8 + 792*b^4*tan(1/2*a)^10 + 66*b^8*tan(1/2*a)^4 - 440*b^6*tan(1/2*a)^6 + 495*b^4*tan(1/2*a)^8 + 12*b^8*tan(1/2*a)^2 - 132*b^6*tan(1/2*a)^4 + 220*b^4*tan(1/2*a)^6 + b^8 - 24*b^6*tan(1/2*a)^2 + 66*b^4*tan(1/2*a)^4 - 2*b^6 + 12*b^4*tan(1/2*a)^2 + b^4))*tan(1/2*x) - (b^8*tan(1/2*a)^25 + 10*b^8*tan(1/2*a)^23 - b^6*tan(1/2*a)^25 + 44*b^8*tan(1/2*a)^21 - 10*b^6*tan(1/2*a)^23 + 110*b^8*tan(1/2*a)^19 - 44*b^6*tan(1/2*a)^21 + 165*b^8*tan(1/2*a)^17 - 110*b^6*tan(1/2*a)^19 + 132*b^8*tan(1/2*a)^15 - 165*b^6*tan(1/2*a)^17 - 132*b^6*tan(1/2*a)^15 - 132*b^8*tan(1/2*a)^11 - 165*b^8*tan(1/2*a)^9 + 132*b^6*tan(1/2*a)^11 - 110*b^8*tan(1/2*a)^7 + 165*b^6*tan(1/2*a)^9 - 44*b^8*tan(1/2*a)^5 + 110*b^6*tan(1/2*a)^7 - 10*b^8*tan(1/2*a)^3 + 44*b^6*tan(1/2*a)^5 - b^8*tan(1/2*a) + 10*b^6*tan(1/2*a)^3 + b^6*tan(1/2*a))/(b^8*tan(1/2*a)^24 + 12*b^8*tan(1/2*a)^22 - 2*b^6*tan(1/2*a)^24 + 66*b^8*tan(1/2*a)^20 - 24*b^6*tan(1/2*a)^22 + b^4*tan(1/2*a)^24 + 220*b^8*tan(1/2*a)^18 - 132*b^6*tan(1/2*a)^20 + 12*b^4*tan(1/2*a)^22 + 495*b^8*tan(1/2*a)^16 - 440*b^6*tan(1/2*a)^18 + 66*b^4*tan(1/2*a)^20 + 792*b^8*tan(1/2*a)^14 - 990*b^6*tan(1/2*a)^16 + 220*b^4*tan(1/2*a)^18 + 924*b^8*tan(1/2*a)^12 - 1584*b^6*tan(1/2*a)^14 + 495*b^4*tan(1/2*a)^16 + 792*b^8*tan(1/2*a)^10 - 1848*b^6*tan(1/2*a)^12 + 792*b^4*tan(1/2*a)^14 + 495*b^8*tan(1/2*a)^8 - 1584*b^6*tan(1/2*a)^10 + 924*b^4*tan(1/2*a)^12 + 220*b^8*tan(1/2*a)^6 - 990*b^6*tan(1/2*a)^8 + 792*b^4*tan(1/2*a)^10 + 66*b^8*tan(1/2*a)^4 - 440*b^6*tan(1/2*a)^6 + 495*b^4*tan(1/2*a)^8 + 12*b^8*tan(1/2*a)^2 - 132*b^6*tan(1/2*a)^4 + 220*b^4*tan(1/2*a)^6 + b^8 - 24*b^6*tan(1/2*a)^2 + 66*b^4*tan(1/2*a)^4 - 2*b^6 + 12*b^4*tan(1/2*a)^2 + b^4))/(4*b^2*tan(1/2*a)^4*tan(1/2*x)^2 - 8*b^2*tan(1/2*a)^3*tan(1/2*x)^3 + 4*b^2*tan(1/2*a)^2*tan(1/2*x)^4 - tan(1/2*a)^4*tan(1/2*x)^4 + 8*b^2*tan(1/2*a)^3*tan(1/2*x) - 16*b^2*tan(1/2*a)^2*tan(1/2*x)^2 - 2*tan(1/2*a)^4*tan(1/2*x)^2 + 8*b^2*tan(1/2*a)*tan(1/2*x)^3 - 2*tan(1/2*a)^2*tan(1/2*x)^4 + 4*b^2*tan(1/2*a)^2 - tan(1/2*a)^4 - 8*b^2*tan(1/2*a)*tan(1/2*x) + 4*b^2*tan(1/2*x)^2 - 4*tan(1/2*a)^2*tan(1/2*x)^2 - tan(1/2*x)^4 - 2*tan(1/2*a)^2 - 2*tan(1/2*x)^2 - 1)

integrate(cos(x)/(1 - (b^2)*(sin(x - a))^2)^(3/2), x, algorithm="fricas")

returns

(2*(b^2*cos(a)^3 - b^2*cos(a))*cos(x)*sin(x) + (b^2*cos(a)^2 - (2*b^2*cos(a)^2 - b^2)*cos(x)^2 - 1)*sin(a) - (b^3*cos(a)*cos(x)*sin(a) - (b^3*cos(a)^2 - b)*sin(x))*sqrt(2*b^2*cos(a)*cos(x)*sin(a)*sin(x) - b^2*cos(a)^2 + (2*b^2*cos(a)^2 - b^2)*cos(x)^2 + 1))/(2*(b^5 - b^3)*cos(a)*cos(x)*sin(a)*sin(x) + b^3 - (b^5 - b^3)*cos(a)^2 - (b^5 - b^3 - 2*(b^5 - b^3)*cos(a)^2)*cos(x)^2 - b)

integrate(cos(x)/(1 - (b^2)*(sin(x - a))^2)^(3/2), x, algorithm="maxima")

returns

RuntimeError: ECL says: sign: argument cannot be imaginary; found %i

integrate(cos(x)/(1 - (b^2)*(sin(x - a))^2)^(3/2), x, algorithm="sympy")

just returns the input.

Additional Information

No response

Environment

- **OS**: Ubuntu 22.04
- **Sage Version**: 10.1

Checklist

fchapoton commented 1 year ago

So what is the exact issue ? Which one of these is incorrect ? Note that sage tries all of these and returns the first that gives an answer. In this case, the answer comes from giac. The answer from sympy is just that it cannot do that, so it's fine.

fchapoton commented 1 year ago

Note that the derivative of fricas answer can be simplified to

cos(x)/(1/2*b^2*cos(-2*a + 2*x) - 1/2*b^2 + 1)^(3/2)
beew commented 1 year ago

So how do you know which one is correct if the answer depends on which algorithms you use? Also the fricas solution is presented in such horrendous form is there a way to simplify it to gives a compact answer like Mathematica gives?

miguelmarco commented 1 year ago

They all seem to be different representations of the same function. So they are all correct. One is just a nicer expression than the others.

One way to simplify it to a not so horrendous form is to use the .simplify_trig() (or .simplify_full()) methods.

beew commented 11 months ago

They all seem to be different representations of the same function. So they are all correct. One is just a nicer expression than the others.

One way to simplify it to a not so horrendous form is to use the .simplify_trig() (or .simplify_full()) methods.

 integrate(cos(x)/(1 - (b^2)*(sin(x - a))^2)^(3/2),x, algorithm="fricas").simplify_trig()

((2*cos(a)^2 - 1)*b^2*cos(x)^2*sin(a) - b^2*cos(a)^2*sin(a) - 2*(cos(a)^3 - cos(a))*b^2*cos(x)*sin(x) + (b^3*cos(a)*cos(x)*sin(a) - (b^3*cos(a)^2 - b)*sin(x))*sqrt(2*b^2*cos(a)*cos(x)*sin(a)*sin(x) + (2*cos(a)^2 - 1)*b^2*cos(x)^2 - b^2*cos(a)^2 + 1) + sin(a))/(b^5*cos(a)^2 - (cos(a)^2 + 1)*b^3 - ((2*cos(a)^2 - 1)*b^5 - (2*cos(a)^2 - 1)*b^3)*cos(x)^2 - 2*(b^5*cos(a)*sin(a) - b^3*cos(a)*sin(a))*cos(x)*sin(x) + b)

integrate(cos(x)/(1 - (b^2)*(sin(x - a))^2)^(3/2),x, algorithm="fricas").simplify_full()

((2*cos(a)^2 - 1)*b^2*cos(x)^2*sin(a) - b^2*cos(a)^2*sin(a) - 2*(cos(a)^3 - cos(a))*b^2*cos(x)*sin(x) + (b^3*cos(a)*cos(x)*sin(a) - (b^3*cos(a)^2 - b)*sin(x))*sqrt(2*b^2*cos(a)*cos(x)*sin(a)*sin(x) + (2*cos(a)^2 - 1)*b^2*cos(x)^2 - b^2*cos(a)^2 + 1) + sin(a))/(b^5*cos(a)^2 - (cos(a)^2 + 1)*b^3 - ((2*cos(a)^2 - 1)*b^5 - (2*cos(a)^2 - 1)*b^3)*cos(x)^2 - 2*(b^5*cos(a)*sin(a) - b^3*cos(a)*sin(a))*cos(x)*sin(x) + b)

None give a nice compact output.