Mathematica Stack Exchange is a question and answer site for users of Mathematica. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

Sorry if this question is too basic -- I'm not very familiar with Mathematica. I am interested in a way to systematically address the following sort of problem:

Suppose that $u=u(x,y)$ is a function of two variables on $\mathbb{R}^2$ that satisfies a Liouville type equation: $$\Delta u= e^{-2u}$$.

I can check by hand that the function $$v=\partial^2_{xy} u+\partial_x u \cdot \partial_y u$$ is harmonic, i.e. $$\Delta v=0.$$

My question is how can I make Mathematica verify this? I'm interested in more complicated examples I don't want to check by hand, but the principle is the same.

share|improve this question
Do you explicitly have $u$, or just the PDE? – J. M. Aug 23 '12 at 13:57
Just the PDE. It should be entirely a symbolic manipulation. The point is that for any $u$ satisfying the PDE the function $v$ is harmonic. You can verify this directly by differentiating. – oberon Aug 23 '12 at 13:59
so you are looking to automatically find a transformation that makes an arbitrary nonlinear PDE into the Laplace equation? – acl Aug 23 '12 at 14:05
@acl No that would be crazy. I'm looking for a way to start with a fixed PDE (the Liuoville equation say) and a solution $u$ to the PDE. I also start with a $v$ which is an explicit algebraic combination of derivatives of $u$ (say the $v$ I wrote down). I want Mathematica to check that $v$ is harmonic using the rules of differentiation and the PDE. – oberon Aug 23 '12 at 14:10
it's easy enough to differentiate the v and insert the explicit expression, but that gives a mass of symbols in which you'd need to replace $u_{xx}+u_{yy}$ by $\exp(-u)$ (not to mention the derivatives) and I don't see how to do that; I hope someone else will – acl Aug 23 '12 at 17:14
up vote 4 down vote accepted

The general idea, as best I understand, is to take your differential identities, generate more by differentiating (prolongation), treat all derivatives as algebraically independent variables, and try to show that the putative identity is in the (algebraic) ideal generated by the known polynomials in these variables. The trick is to know how far to prolong.

While there are algorithmic approaches (I believe this goes by the name "Cartan Kuranishi algorithm, or something to that effect), I'll just assume that it will suffice to get all second derivatives of our defining equations. Here is the code.

The independent ariables and initial differential "polynomials".

ivars = {x, y};
dpolys = {D[u[x, y], {x, 2}] + D[u[x, y], {y, 2}] - Exp[-2*u[x, y]],
   v[x, y] - (D[u[x, y], x, y] + D[u[x, y], x]*D[u[x, y], y])};

Get the Laplacian of v.

lapv = Total[Flatten[Outer[D[#1, {#2, 2}] &, {v[x, y]}, ivars]]];

Now get all first and second derivatives of these diff polys.

allpolys = Flatten[Join[dpolys, Outer[D, dpolys, ivars],
    Outer[D[#1, {#2, 2}] &, dpolys, ivars], D[dpolys, x, y]]];
allvars = Variables[allpolys];

I want to impose an ordering on the "variables". I'm not sure this is necessary but at least it will help to make sure we use the same term ordering when we attempt to reduce the Laplacian of v to zero.

uvars = Select[allvars, ! FreeQ[#, u] &];
duvars = Select[uvars, ! FreeQ[#, Derivative] &];

vvars = Select[allvars, ! FreeQ[#, v] &]; dvvars = 
 Select[vvars, ! FreeQ[#, Derivative] &];

derivorder[Derivative[a : __][_][__]] := Plus[a]

sortuvars = 
  Append[Sort[duvars, derivorder[#1] >= derivorder[#2] &], u[x, y]];
sortvvars = 
  Append[Sort[dvvars, derivorder[#1] >= derivorder[#2] &], v[x, y]];

Small problem: when we go to compute a Groebner basis, the exponentials will not correctly be seen as invertible paramaters. So I replace them up front with a new symbol. We will treat it as invertible. It turns out that this both makes the reduction step viable and also makes the basis computation hugely faster.

allpolys2 = allpolys /. E^(-2 u[x, y]) -> exp2u;

Compute a(n algebraic) Groebner basis for our polynomials. Again, while they involve derivatives, we treat those as being algebraically independent variables.

 gb = GroebnerBasis[allpolys2, Join[sortvvars, sortuvars], 
    MonomialOrder -> DegreeReverseLexicographic, 
    CoefficientDomain -> RationalFunctions];]

(* Out[615]= {0.11, Null} *)

Now check whether the Laplacian of v lives in this (algebraic) ideal. If so, that proves that its vanishing is a consequence of the original differential equations.

PolynomialReduce[lapv, gb, Join[sortvvars, sortuvars], 
  MonomialOrder -> DegreeReverseLexicographic, 
  CoefficientDomain -> RationalFunctions][[2]]

(* Out[619]= 0 *)

[Last step: Declare victory and go home.]

share|improve this answer
This sounds like the correct approach. I will investigate further and see if I can get it to work. Thanks! – oberon Aug 26 '12 at 9:56

I am fairly certain there is a better way to do this, but you can try the following which involves using the VectorAnalysis package:


\[CapitalDelta]u[x_, y_] := Exp[-2 u[x, y]]
v[x_, y_] := Derivative[1, 1][u][x,y] + (Derivative[1, 0][u][x,y])*(Derivative[0, 1][u][x, y])  

r1 = Laplacian[u[x, y], Cartesian[x, y, z]] -> \[CapitalDelta]u[x,y];
dx = D[Laplacian[u[x, y], Cartesian[x, y, z]], x] -> D[\[CapitalDelta]u[x, y], x];
dy = D[Laplacian[u[x, y], Cartesian[x, y, z]], y] -> D[\[CapitalDelta]u[x, y], y];
dxy = D[Laplacian[u[x, y], Cartesian[x, y, z]], x, y] -> D[\[CapitalDelta]u[x, y], x, y];

FullSimplify[Laplacian[v[x, y], Cartesian[x, y, z]]] //. {dx, dy, dxy,r1}
share|improve this answer

Your Answer


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.