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 [email protected]
To unsubscribe from this group, send email to
[email protected]
For more options, visit this group at
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org
