# Eliminating functions from system of PDE in Mathematica

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

• 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?) Commented Apr 29, 2013 at 21:04
• Crikey. I think I got something that works. Surprised am I. (Or maybe just seeing things.) Will post as response. Commented Apr 29, 2013 at 21:05

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