Problem with extracting coefficients in a formal expression
I am computing recursively some big formal power series $f$ that depends on many different variables. It is a polynomial in $p_1,p_2,p_3,q_1,q_2,q_3,r_1,r_2,r_3$, whose coefficients are rational functions in $Q_1,Q_2,Q_3,h,a,L$. Eventually I set $Q_1=-1,Q_2=0,Q_3=1$.
I use variables $y_1,y_2,\dots$ as unknown coefficients that I want to compute by acting on $f$ by some operator $W_{2,3}$ and comparing the coefficients of $W_{2,3}(f)$ as a polynomial in $p_1,p_2,p_3,q_1,q_2,q_3,r_1,r_2,r_3$. Let's say that I want to find a coefficient of the quadratic term $p_1q_1$. In order to do so I use a function
.coefficient(q1,1).coefficient(p1,1)
but I have noticed that when I do this, I do not get a correct answer. Moreover, when I change the order of variables and I use a function
.coefficient(p1,1).coefficient(q1,1)
I have got a different answer, and still not correct. I am not very in programming, so I cannot figure out what is the reason, but it seems that there must by a problem with how the computer treats my functions. This is because When I act on this function by an operator that annihilates the function, the program gives me a huge expansion which is not equal to zero even if I use commands like expand or simplify.
Here is my code:
var('L,Q1,Q2,Q3,h,a,p1,p2,p3,q1,q2,q3,r1,r2,r3,y1,y2,y3,y4,y5,y6,y7,y8,y9,y10,y11,y12,y13,y14,y15,y16,y17,y18,y19,y20,y21,y22')
def pow(n):
return "p"+str(n)
def poww(n):
return eval(pow(n))
def qow(n):
return "q"+str(n)
def qoww(n):
return eval(qow(n))
def row(n):
return "r"+str(n)
def roww(n):
return eval(row(n))
def J1(f,i,n):
return h^2*f.diff(poww(i))
def J2(f,i,n):
return h^2*f.diff(qoww(i))
def J3(f,i,n):
return h^2*f.diff(roww(i))
def J12(f,i,n):
return h^2*((Q2-h*a)*f.diff(poww(i))+(Q1-h*a)*f.diff(qoww(i))+h^2*sum([(f.diff(poww(k))).diff(qoww(i-k)) for k in range (1,i)])+sum([k*(poww(k)*f.diff(qoww(i+k))+qoww(k)*f.diff(poww(i+k))) for k in range (1,n-i)]))
def J13(f,i,n):
return h^2*((Q3-h*a)*f.diff(poww(i))+(Q1-h*a)*f.diff(roww(i))+h^2*sum([(f.diff(poww(k))).diff(roww(i-k)) for k in range (1,i)])+sum([k*(poww(k)*f.diff(roww(i+k))+roww(k)*f.diff(poww(i+k))) for k in range (1,n-i)]))
def J23(f,i,n):
return h^2*((Q3-h*a)*f.diff(qoww(i))+(Q2-h*a)*f.diff(roww(i))+h^2*sum([(f.diff(qoww(k))).diff(roww(i-k)) for k in range (1,i)])+sum([k*(qoww(k)*f.diff(roww(i+k))+roww(k)*f.diff(qoww(i+k))) for k in range (1,n-i)]))
def JJ12(f,i,n):
return J12(f,i,n).subs(Q1=-1,Q2=0,Q3=1)
def JJ13(f,i,n):
return J13(f,i,n).subs(Q1=-1,Q2=0,Q3=1)
def JJ23(f,i,n):
return J23(f,i,n).subs(Q1=-1,Q2=0,Q3=1)
def W23(f,i,n):
return J12(f,i,n)+J13(f,i,n)+J23(f,i,n)-h*a*(i+1)*(J2(f,i,n)+2*J3(f,i,n))
def WW23(f,i,n):
return JJ12(f,i,n)+JJ13(f,i,n)+JJ23(f,i,n)-h*a*(i+1)*(J2(f,i,n)+2*J3(f,i,n))
f = (3*a*h - Q1 + Q3)*(2*a*h - Q1 + Q2)*(2*a*h - Q2 + Q3)+L*(2*a*h - Q2 + Q3)/h^2*p1-L*(4*a*h - Q1 + Q3)/h^2*q1+L*(2*a*h - Q1 + Q2)/h^2*r1-(y7+y8)/2*p1^2-(y7+y9)/2*q1^2-(y8+y9)/2*r1^2+y7*p1*q1+y8*p1*r1+y9*q1*r1+y10*p2-(y10+y12)*q2+y12*r2
g=f.subs(y7 = (L^2*Q1^3 - L^2*Q1^2*Q2 - 2*(24*L^2*a^3 - L^2*a)*h^3 + (L^2*Q1 - L^2*Q2)*Q3^2 - (11*L^2*Q1^2 - 8*L^2*Q1*Q2 + 3*L^2*Q3^2 - 2*(7*L^2*Q1 - 4*L^2*Q2)*Q3)*a*h - (L^2*Q1 - L^2*Q2 - 8*(5*L^2*Q1 - 2*L^2*Q2 - 3*L^2*Q3)*a^2)*h^2 - 2*(L^2*Q1^2 - L^2*Q1*Q2)*Q3)/((72*a^4 - 18*a^2 + 1)*h^8 - (6*(17*Q1 - 10*Q2 - 7*Q3)*a^3 - (12*Q1 - 5*Q2 - 7*Q3)*a)*h^7 - (12*Q1^3 - 19*Q1^2*Q2 + 7*Q1*Q2^2 + 5*(Q1 - Q2)*Q3^2 - (17*Q1^2 - 24*Q1*Q2 + 7*Q2^2)*Q3)*a*h^5 + ((53*Q1^2 - 59*Q1*Q2 + 12*Q2^2 - (47*Q1 - 35*Q2)*Q3 + 6*Q3^2)*a^2 - 2*Q1^2 + 2*Q1*Q2 - Q2^2 + 2*Q1*Q3 - Q3^2)*h^6 + (Q1^4 - 2*Q1^3*Q2 + Q1^2*Q2^2 + (Q1^2 - 2*Q1*Q2 + Q2^2)*Q3^2 - 2*(Q1^3 - 2*Q1^2*Q2 + Q1*Q2^2)*Q3)*h^4), y8 = -(4*L^2*a*h - L^2*Q1 + L^2*Q3)/(7*(Q1 - Q3)*a*h^5 - (12*a^2 - 1)*h^6 - (Q1^2 - 2*Q1*Q3 + Q3^2)*h^4), y9 = (L^2*Q1^2*Q2 - L^2*Q3^3 - 2*(24*L^2*a^3 - L^2*a)*h^3 + (2*L^2*Q1 + L^2*Q2)*Q3^2 - (3*L^2*Q1^2 + 8*L^2*Q1*Q2 + 11*L^2*Q3^2 - 2*(7*L^2*Q1 + 4*L^2*Q2)*Q3)*a*h - (L^2*Q2 - L^2*Q3 - 8*(3*L^2*Q1 + 2*L^2*Q2 - 5*L^2*Q3)*a^2)*h^2 - (L^2*Q1^2 + 2*L^2*Q1*Q2)*Q3)/((72*a^4 - 18*a^2 + 1)*h^8 - (6*(7*Q1 + 10*Q2 - 17*Q3)*a^3 - (7*Q1 + 5*Q2 - 12*Q3)*a)*h^7 - (5*Q1^2*Q2 + 7*Q1*Q2^2 + (17*Q1 + 19*Q2)*Q3^2 - 12*Q3^3 - (5*Q1^2 + 24*Q1*Q2 + 7*Q2^2)*Q3)*a*h^5 + ((6*Q1^2 + 35*Q1*Q2 + 12*Q2^2 - (47*Q1 + 59*Q2)*Q3 + 53*Q3^2)*a^2 - Q1^2 - Q2^2 + 2*(Q1 + Q2)*Q3 - 2*Q3^2)*h^6 + (Q1^2*Q2^2 - 2*(Q1 + Q2)*Q3^3 + Q3^4 + (Q1^2 + 4*Q1*Q2 + Q2^2)*Q3^2 - 2*(Q1^2*Q2 + Q1*Q2^2)*Q3)*h^4), y10 = (12*L^2*a^2*h^2 + 2*L^2*Q1*Q2 - L^2*Q2^2 - 2*L^2*Q1*Q3 + L^2*Q3^2 - 4*(L^2*Q1 + L^2*Q2 - 2*L^2*Q3)*a*h)/((72*a^4 - 18*a^2 + 1)*h^6 - (6*(17*Q1 - 10*Q2 - 7*Q3)*a^3 - (12*Q1 - 5*Q2 - 7*Q3)*a)*h^5 - (12*Q1^3 - 19*Q1^2*Q2 + 7*Q1*Q2^2 + 5*(Q1 - Q2)*Q3^2 - (17*Q1^2 - 24*Q1*Q2 + 7*Q2^2)*Q3)*a*h^3 + ((53*Q1^2 - 59*Q1*Q2 + 12*Q2^2 - (47*Q1 - 35*Q2)*Q3 + 6*Q3^2)*a^2 - 2*Q1^2 + 2*Q1*Q2 - Q2^2 + 2*Q1*Q3 - Q3^2)*h^4 + (Q1^4 - 2*Q1^3*Q2 + Q1^2*Q2^2 + (Q1^2 - 2*Q1*Q2 + Q2^2)*Q3^2 - 2*(Q1^3 - 2*Q1^2*Q2 + Q1*Q2^2)*Q3)*h^2), y12 = -(12*L^2*a^2*h^2 + L^2*Q1^2 - L^2*Q2^2 - 4*(2*L^2*Q1 - L^2*Q2 - L^2*Q3)*a*h - 2*(L^2*Q1 - L^2*Q2)*Q3)/((72*a^4 - 18*a^2 + 1)*h^6 - (6*(7*Q1 + 10*Q2 - 17*Q3)*a^3 - (7*Q1 + 5*Q2 - 12*Q3)*a)*h^5 - (5*Q1^2*Q2 + 7*Q1*Q2^2 + (17*Q1 + 19*Q2)*Q3^2 - 12*Q3^3 - (5*Q1^2 + 24*Q1*Q2 + 7*Q2^2)*Q3)*a*h^3 + ((6*Q1^2 + 35*Q1*Q2 + 12*Q2^2 - (47*Q1 + 59*Q2)*Q3 + 53*Q3^2)*a^2 - Q1^2 - Q2^2 + 2*(Q1 + Q2)*Q3 - 2*Q3^2)*h^4 + (Q1^2*Q2^2 - 2*(Q1 + Q2)*Q3^3 + Q3^4 + (Q1^2 + 4*Q1*Q2 + Q2^2)*Q3^2 - 2*(Q1^2*Q2 + Q1*Q2^2)*Q3)*h^2))/((3*a*h - Q1 + Q3)*(2*a*h - Q1 + Q2)*(2*a*h - Q2 + Q3))
gg = g.subs(Q1=-1,Q2=0,Q3=1)
ggg = gg+y1*p3+y2*q3+y3*r3+y4*p2*p1+y5*p2*q1+y6*p2*r1+y7*q2*p1+y8*q2*q1+y9*q2*r1+y10*r2*p1+y11*r2*q1+y12*r2*r1+y13*p1^3+y14*q1^3+y15*r1^3+y16*p1^2*q1+y17*p1^2*r1+y18*q1^2*p1+y19*q1^2*r1+y20*r1^2*p1+y21*r1^2*q1+y22*p1*q1*r1
WW23(ggg,1,7).coefficient(q1,1).coefficient(p1,1)
The last command produces
-4*(h^2*y18 + h^2*y22)*a*h - (2*a*h*y16 + 2*(a*h + 1)*y18 - y8)*h^2 - (a*h*y22 + 2*(a*h - 1)*y18)*h^2 - (2*(a*h - 1)*y16 + (a*h + 1)*y22 - y11)*h^2
which is incorrect. If I change the order and use
WW23(ggg,1,7).coefficient(p1,1).coefficient(q1,1)
I have
-4*(h^2*y18 + h^2*y22)*a*h - (2*a*h*y16 + 2*(a*h + 1)*y18 - y4)*h^2 - (a*h*y22 + 2*(a*h - 1)*y18 - y10)*h^2 - (2*(a*h - 1)*y16 + (a*h + 1)*y22)*h^2
Both expressions are different and they are not giving me a coefficient of the quadratic term p1*q1 in WW23(ggg). What is the reason?
