Mathematica users,

Let's see an example of three PDEs:

How can I eliminate variables q1[x,y,z] and q2[x,y,z] from a system of three PDEs in Mathematica. In final form, I need just one equation in function of p1[x,y,z]. Why this code is working well when I want to eliminate q1 and q2 from a1, a2 and a3, but doesn't work when I want to eliminate q1 and q2 from a1, a2, a3 and a4?

    D[p1[x,y,z],{y,2}]+D[q1[x,y,z], {x,1}]+D[q2[x,y,z],{y,1}])==0;




 rule=Flatten[{#[x,y,z]->#,Derivative[n_,m_,o_][#][x_,y_,z_]-># x^n y^m z^o}&/@{p1,q1,q2}];

 invRule=(x_->z_):>z Derivative[Sequence@@x][p1];




  • $\begingroup$ I'm not convinced that "working" example really works. You seem to replace mixed derivatives of the qj's with qj (not a derivative thereof) times monomials. Is that really correct for your example? (Are your qj functions known to be exponentials of sums of squares?) $\endgroup$ Apr 29, 2013 at 21:04
  • $\begingroup$ Crikey. I think I got something that works. Surprised am I. (Or maybe just seeing things.) Will post as response. $\endgroup$ Apr 29, 2013 at 21:05

1 Answer 1


We start with the given differential polynomials.

a1 = A1*D[p1[x, y, z], {z, 2}] - 
   A2*(D[p1[x, y, z], {x, 2}] + D[p1[x, y, z], {y, 2}] + 
      D[q1[x, y, z], {x, 1}] + D[q2[x, y, z], {y, 1}]);
a2 = A9*(A2*
       D[(D[q1[x, y, z], {x, 1}] + D[q2[x, y, z], {y, 1}]), {x, 1}] + 
      A2*(D[q1[x, y, z], {x, 2}] + D[q1[x, y, z], {y, 2}])) - 
   A4*(D[p1[x, y, z], {x, 1}] + q1[x, y, z]) - 
   A8*D[q1[x, y, z], {z, 2}];
a3 = A9*(A6*
       D[(D[q1[x, y, z], {x, 1}] + D[q1[x, y, z], {y, 1}]), {y, 1}] + 
      A7*(D[q2[x, y, z], {x, 2}] + D[q1[x, y, z], {y, 2}])) - 
   A1*(D[p1[x, y, z], {y, 1}] + q2[x, y, z]) - 
   A4*D[q2[x, y, z], {z, 2}];
a4 = A09*A006*
    D[(D[q2[x, y, z], {x, 1}] + D[q2[x, y, z], {y, 1}]), {y, 1}] + 
   A004*D[q1[x, y, z], {z, 2}];

Take three prolongations (that is, derivatives with respect to each variable in turn).

diffpolys = {a1, a2, a3, a4};
vars = {x, y, z};
d2 = Join[diffpolys, Flatten[Outer[D, diffpolys, vars]]];
d3 = Union[Join[d2, Flatten[Outer[D, d2, vars]]]];
d4 = Union[Join[d3, Flatten[Outer[D, d3, vars]]]];

pvars = Cases[Variables[d4], Derivative[__][p1][__]];
qvars = Cases[Variables[d4], Derivative[__][q1 | q2][__]];

(* Out[135]= 80 *)

(* Out[136]= 110 *)

110 unknowns in 80 polynomials is not promising. Undaunted, we forge ahead.

 gb = GroebnerBasis[d4, pvars, qvars, 
    MonomialOrder -> EliminationOrder, 
    CoefficientDomain -> RationalFunctions];]

(* Out[137]= {7.540000, Null} *)

So here it is.


{A006*A09*A1*A4*Derivative[0, 1, 2][p1][x, y, z] + A006*A09*A1*A8*Derivative[0, 1, 4][p1][x, y, z] - A006*A09*A2*A4*Derivative[0, 3, 0][p1][x, y, z] + ((-A006)*A09*A2*A8 - A006*A09*A1*A2*A9)*Derivative[0, 3, 2][p1][ x, y, z] + A006*A09*A2^2*A9*Derivative[0, 5, 0][p1][x, y, z] + (A006*A09*A1*A4 - A004*A2*A4)* Derivative[1, 0, 2][p1][x, y, z] + (A006*A09*A1*A8 + A004*A1*A2*A9)*Derivative[1, 0, 4][p1][x, y, z] - A006*A09*A2*A4*Derivative[1, 2, 0][p1][x, y, z] + ((-A006)*A09*A2*A8 - A006*A09*A1*A2*A9 - A004*A2^2*A9)* Derivative[1, 2, 2][p1][x, y, z] + A006*A09*A2^2*A9* Derivative[1, 4, 0][p1][x, y, z] + ((-A006)*A09*A2*A8 - 2*A006*A09*A1*A2*A9)* Derivative[2, 1, 2][p1][x, y, z] + 3*A006*A09*A2^2*A9* Derivative[2, 3, 0][p1][x, y, z] + ((-A006)*A09*A2*A8 - 2*A006*A09*A1*A2*A9 - A004*A2^2*A9)* Derivative[3, 0, 2][p1][x, y, z] + 3*A006*A09*A2^2*A9* Derivative[3, 2, 0][p1][x, y, z] + 2*A006*A09*A2^2*A9* Derivative[4, 1, 0][p1][x, y, z] + 2*A006*A09*A2^2*A9* Derivative[5, 0, 0][p1][x, y, z]}

This does seem to match your result. The fact that I do not understand why your method works is no indication of shortcoming (other than perhaps on my part).


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.