Hello,
I used sage to find the root (named P2plus) of an equation. The variable is 
P1plus and there are 3 parameters, M1, M2 and gamma. The expression of the 
root (P2plus) being complicated, and as I am interested in its behavior for 
small values of P1plus, I asked Sage for the Taylor expansion (2nd order) 
of P2plus. The problem is, when I plot the initial expression of P2plus and 
its Taylor expansion, one can clearly see that the expressions differ near 
P1plus=0. In other words, the Taylor expansion given by Sage seems to be 
erroneous.
The root P2plus being a fraction, I looked at the Taylor expansions of the 
numerator and denominator. It appears that the Taylor expansion of the 
numerator is erroneous. But I can't figure out why.

Here is my code. Is there any error in it or am I doing something wrong 
while asking for the Taylor expansion ? Any help would be greatly 
appreciated.
Thanks
Maxime



P1plus  = var('P1plus')
M1      = var('M1')
M2      = var('M2')
gamma   = var('gamma')


P2plus_init = (M1^3*M2^2*gamma + M1^2*M2^3*gamma + 2*M1*M2^3*P1plus*gamma - 
M1^3*M2^2 + M1^3*M2*gamma - M1^2*M2^3 - 2*M1*M2^3*P1plus - M1*M2^3*gamma + 
2*M1*M2^2*P1plus*gamma - M1^3*M2 - M1^2*M2*gamma + M1*M2^3 - 
2*M1*M2^2*P1plus - M1*M2^2*gamma + 3*M1^2*M2 + 3*M1*M2^2 + 4*M1*M2*P1plus + 
2*M1^2 + 4*M1*P1plus - 2*M1 - 2*M2 + sqrt((M1^4*P1plus^2 + 2*(2*P1plus^3 - 
P1plus^2)*M1^3 + (4*P1plus^4 - 4*P1plus^3 + P1plus^2)*M1^2)*M2^6*gamma^4 + 
((P1plus^2 + 4*P1plus + 1)*M1^4 + 2*(2*P1plus^3 + 6*P1plus^2 - 1)*M1^3 + 
(4*P1plus^4 + 12*P1plus^3 + P1plus^2 - 8*P1plus + 1)*M1^2 + 2*(4*P1plus^4 - 
5*P1plus^2 + 2*P1plus)*M1)*M2^6 + 2*((P1plus^2 + 6*P1plus - 1)*M1^4 + M1^5 
+ (4*P1plus^3 + 12*P1plus^2 - 2*P1plus - 3)*M1^3 + (4*P1plus^4 + 
12*P1plus^3 + P1plus^2 - 12*P1plus + 5)*M1^2 + 2*(4*P1plus^4 - 5*P1plus^2 + 
4*P1plus - 1)*M1)*M2^5 + (M1^6 - 2*(2*P1plus^2 + 4*P1plus + 7)*M1^4 + 
2*M1^5 - 2*(8*P1plus^3 + 24*P1plus^2 + 2*P1plus - 5)*M1^3 - (16*P1plus^4 + 
48*P1plus^3 - 12*P1plus - 13)*M1^2 - 8*(4*P1plus^4 - 5*P1plus^2 + 2)*M1 + 
4)*M2^4 + 8*(2*P1plus - 1)*M1^3 - 2*((2*(P1plus^2 + P1plus)*M1^4 + 
(8*P1plus^3 + 3*P1plus^2 - 2*P1plus)*M1^3 + 2*(4*P1plus^4 - P1plus)*M1^2 + 
(4*P1plus^4 - 5*P1plus^2 + 2*P1plus)*M1)*M2^6 + (M1^4*P1plus^2 + 
2*(2*P1plus^3 - P1plus^2)*M1^3 + (4*P1plus^4 - 4*P1plus^3 + 
P1plus^2)*M1^2)*M2^5 - 2*(M1^4*P1plus^2 + 2*(2*P1plus^3 - P1plus^2)*M1^3 + 
(4*P1plus^4 - 4*P1plus^3 + P1plus^2)*M1^2)*M2^4)*gamma^3 + 2*(M1^6 - 
(4*P1plus^2 + 18*P1plus + 5)*M1^4 - 3*M1^5 - (16*P1plus^3 + 48*P1plus^2 - 
2*P1plus - 17)*M1^3 - 4*(4*P1plus^4 + 12*P1plus^3 + P1plus^2 - 9*P1plus + 
2)*M1^2 - 2*(16*P1plus^4 - 20*P1plus^2 + 10*P1plus + 3)*M1 + 4)*M2^3 + 
4*M1^4 + 4*(4*P1plus^2 - 4*P1plus + 1)*M1^2 + (((6*P1plus^2 + 12*P1plus + 
1)*M1^4 + 2*(12*P1plus^3 + 15*P1plus^2 - 4*P1plus - 1)*M1^3 + (24*P1plus^4 
+ 24*P1plus^3 - 2*P1plus^2 - 16*P1plus + 1)*M1^2 + 6*(4*P1plus^4 - 
5*P1plus^2 + 2*P1plus)*M1)*M2^6 + 2*((3*P1plus^2 + 6*P1plus - 1)*M1^4 + 
M1^5 + (12*P1plus^3 + 8*P1plus^2 - 2*P1plus - 1)*M1^3 + (12*P1plus^4 + 
4*P1plus^3 + 3*P1plus^2 - 8*P1plus + 1)*M1^2 + 2*(4*P1plus^4 - 5*P1plus^2 + 
2*P1plus)*M1)*M2^5 + (M1^6 - 2*(6*P1plus^2 + 4*P1plus + 3)*M1^4 + 2*M1^5 - 
2*(24*P1plus^3 + 16*P1plus^2 - 6*P1plus - 1)*M1^3 - (48*P1plus^4 + 
16*P1plus^3 - 8*P1plus^2 - 12*P1plus - 1)*M1^2 - 8*(4*P1plus^4 - 5*P1plus^2 
+ 2*P1plus)*M1)*M2^4 + 2*(M1^6 - (4*P1plus^2 - 2*P1plus + 1)*M1^4 - M1^5 - 
(16*P1plus^3 - 8*P1plus^2 + 2*P1plus - 1)*M1^3 - 4*(4*P1plus^4 - 4*P1plus^3 
+ P1plus^2)*M1^2)*M2^3 + (M1^6 + (4*P1plus^2 + 1)*M1^4 - 2*M1^5 + 
8*(2*P1plus^3 - P1plus^2)*M1^3 + 4*(4*P1plus^4 - 4*P1plus^3 + 
P1plus^2)*M1^2)*M2^2)*gamma^2 + (M1^6 + (4*P1plus^2 + 13)*M1^4 - 10*M1^5 + 
16*(P1plus^3 + 3*P1plus^2 + 1)*M1^3 + 4*(4*P1plus^4 + 12*P1plus^3 - 
3*P1plus^2 + 4*P1plus - 8)*M1^2 + 8*(4*P1plus^4 - 5*P1plus^2 - 2*P1plus + 
1)*M1 + 4)*M2^2 + 4*(2*(P1plus^2 + 3*P1plus + 2)*M1^4 - M1^5 + (8*P1plus^3 
+ 24*P1plus^2 + 2*P1plus - 3)*M1^3 + 2*(4*P1plus^4 + 12*P1plus^3 + P1plus^2 
- 6*P1plus - 1)*M1^2 + 2*(8*P1plus^4 - 10*P1plus^2 + 2*P1plus + 1)*M1)*M2 - 
2*(((2*P1plus^2 + 6*P1plus + 1)*M1^4 + (8*P1plus^3 + 17*P1plus^2 - 2*P1plus 
- 2)*M1^3 + (8*P1plus^4 + 16*P1plus^3 - 10*P1plus + 1)*M1^2 + 3*(4*P1plus^4 
- 5*P1plus^2 + 2*P1plus)*M1)*M2^6 + ((3*P1plus^2 + 12*P1plus - 2)*M1^4 + 
2*M1^5 + 2*(6*P1plus^3 + 11*P1plus^2 - 2*P1plus - 2)*M1^3 + (12*P1plus^4 + 
20*P1plus^3 + 3*P1plus^2 - 20*P1plus + 6)*M1^2 + 2*(8*P1plus^4 - 
10*P1plus^2 + 6*P1plus - 1)*M1)*M2^5 + (M1^6 - 2*(3*P1plus^2 + 4*P1plus + 
5)*M1^4 + 2*M1^5 - 2*(12*P1plus^3 + 22*P1plus^2 - 2*P1plus - 3)*M1^3 - 
(24*P1plus^4 + 40*P1plus^3 - 6*P1plus^2 - 12*P1plus - 5)*M1^2 - 
4*(8*P1plus^4 - 10*P1plus^2 + 2*P1plus + 1)*M1)*M2^4 + 2*(M1^6 - 
(4*P1plus^2 + 8*P1plus + 3)*M1^4 - 2*M1^5 - (16*P1plus^3 + 20*P1plus^2 - 
7)*M1^3 - 2*(8*P1plus^4 + 8*P1plus^3 + 2*P1plus^2 - 7*P1plus + 1)*M1^2 - 
(16*P1plus^4 - 20*P1plus^2 + 6*P1plus + 1)*M1)*M2^3 + (M1^6 + (4*P1plus^2 + 
5)*M1^4 - 6*M1^5 + 4*(4*P1plus^3 + 5*P1plus^2 - 2*P1plus + 1)*M1^3 + 
4*(4*P1plus^4 + 4*P1plus^3 - 3*P1plus^2 - 1)*M1^2 + 4*(4*P1plus^4 - 
5*P1plus^2 + 2*P1plus)*M1)*M2^2 + 2*(2*(P1plus^2 - P1plus + 1)*M1^4 - M1^5 
+ (8*P1plus^3 - 4*P1plus^2 + 2*P1plus - 1)*M1^3 + 2*(4*P1plus^4 - 
4*P1plus^3 + P1plus^2)*M1^2)*M2)*gamma) - 2*M2^2)/((M1 + 2*P1plus - 
1)*(M2*gamma - M2 - 2)*(M2^2*gamma - M2^2 + 2)*M1)



# Taylor expansions
# -----------------

P2plus_taylor = taylor(P2plus_init,P1plus,0,2)

print "comparison between initial expression (black) and Taylor expansion 
(red)"
figP2plus =             plot(P2plus_init  
(M1=0.3,M2=0.7,gamma=1.4),P1plus,0,0.5,color='black')
figP2plus = figP2plus + 
plot(P2plus_taylor(M1=0.3,M2=0.7,gamma=1.4),P1plus,0,0.5,color='red',linestyle=':')
figP2plus.show()

num = P2plus_init.numerator  ()
den = P2plus_init.denominator()
num_taylor = taylor(num,P1plus,0,2)
den_taylor = taylor(den,P1plus,0,2)

print "comparison between initial expression (black) and Taylor expansion : 
numerator part only (red)"
figNum =          plot(num       
(M1=0.3,M2=0.7,gamma=1.4),P1plus,0,0.5,color='black')
figNum = figNum + 
plot(num_taylor(M1=0.3,M2=0.7,gamma=1.4),P1plus,0,0.5,color='red',linestyle=':')
figNum.show()

print "comparison between initial expression (black) and Taylor expansion : 
denominator part only (red)"
figDen =          plot(den       
(M1=0.3,M2=0.7,gamma=1.4),P1plus,0,0.5,color='black')
figDen = figDen + 
plot(den_taylor(M1=0.3,M2=0.7,gamma=1.4),P1plus,0,0.5,color='red',linestyle=':')
figDen.show()

-- 
-- 
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org



Reply via email to