Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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?

 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}])==0;


 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}]==0;


 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}]==0;

 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}]==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];

 ww1=Eliminate[{a1,a2,a4}/.rule,{q1,q2}]/.(x_==y_)->(x-y);

 ww2=CoefficientRules[ww1,{x,y,z}];

 t=Total[ww2/.p1->1/.invRule]==0;

 t
share|improve this question
    
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?) –  Daniel Lichtblau Apr 29 '13 at 21:04
    
Crikey. I think I got something that works. Surprised am I. (Or maybe just seeing things.) Will post as response. –  Daniel Lichtblau Apr 29 '13 at 21:05

1 Answer 1

up vote 3 down vote accepted

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][__]];

Length[d4]
(* Out[135]= 80 *)

Length[qvars]
(* Out[136]= 110 *)

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

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

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

So here it is.

gb

{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).

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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