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