# Looking for a possible Lagrangian for a PDE using Mathematica

I have the following second order PDE for which I would like to find the Lagrangian (if one exists): $$u_{rr}+u_{zz}+u_{r}/r-u/r^2+u(u_r ^2+u_z ^2+uu_r/r)=0$$The only way that I know to do this is by trial and error, guessing a Lagrangian and then applying EulerEquations[ ] to it to see if the result is the PDE. Or alternatively, can one show that no Lagrangian exists for this? Any ideas?

Since this is a Mathematica site, I will go about this with an ansatz and see what I get. Turns out I get something close to OPs equation, but currently not exactly OPs equation. I will not use the Variational Methods Package that OP refers to, except in a comment at the end.

(* notation *)
independentVars={r,z};
dependentVars={u[r,z]};

(* these determine the ansatz, one can adjust them *)
coefficientFunctions={1,r,1/r};
maxDegree=3;

(* auxiliary notation *)
mons=Join[dependentVars,Flatten[Outer[D,dependentVars,independentVars],1]];
poly=DeleteDuplicates[Times@@@Tuples[Join[{1},mons],maxDegree]];
allTerms=Flatten[Outer[Times,coefficientFunctions,poly]];


To see what the last three mean, just type poly and hit shift-enter to see. The strategy will be to make an ansatz for the Lagrangian density, given by a big linear combination of terms with coefficients a[n], and then see if we can match the resulting Euler-Lagrange equations to what OP wants.

Let me call $$B$$ the left hand side of the equation that OP has written down. One question is whether we need to produce $$B$$ exactly, or if it would be good enough to produce $$MB$$ for some multiplier $$M$$. As long as $$M \neq 0$$, the equation $$B=0$$ and $$MB=0$$ are of course equivalent. I am going to assume that we have the freedom to introduce such a factor $$M$$ here. I will use the following multiplier, which I picked after some experiments:

multiplier=Exp[u[r,z]^2];


We continue with code:

(* Lagrangian density ansatz with parameters a[n] *)
ansatzLagrangian=multiplier*Sum[a[n]*allTerms[[n]],{n,1,Length[allTerms]}];
avars=Table[a[n],{n,1,Length[allTerms]}];

(* resulting Euler-Lagrange equations *)
EulerLagrangeEquations=Map[Function[{dv},
D[ansatzLagrangian,dv]
-Total[Map[D[D[ansatzLagrangian,D[dv,#]],#]&,
independentVars]]],dependentVars]//Expand;

(* left hand side of equation given by OP times multiplier *)
desiredEquation=r*multiplier*(
D[u[r,z],{r,2}]
+D[u[r,z],{z,2}]
+D[u[r,z],r]/r
-c[1]*u[r,z]/r^2
+u[r,z]*(c[2]*D[u[r,z],r]^2+c[3]*D[u[r,z],z]^2+c[4]*u[r,z]*D[u[r,z],r]/r)
);


I introduced an extra factor r, so the multiplier we are using is actually r*multiplier. Again, does not change the equation. We also introduced constants c[1] through c[4]. Now OP would want all four to be equal to 1 but that did not work in my first attempt, so I added these constants to see if I can at least get something similar. Now let us solve:

zeroGoal=(First[EulerLagrangeEquations]-desiredEquation)/multiplier//Expand;
zeroGoal=r^10*zeroGoal//Expand; (* make polynomial in r *)
eqs=Map[Last,CoefficientRules[zeroGoal,DeleteCases[Variables[zeroGoal],(a|c)[_]]]];

(* c[2]==1&&c[3]==1&&c[4]==0 *)


So the result is that with the ansatz that we made, there is no condition on c[1], therefore we can take c[1]=1 if we want. And we get c[2]=1 and c[3]=1 as requested, but we only get c[4]=0 with this ansatz. So this means OPs equation but without the $$u^2u_r/r$$ term.

Code to produce the Lagrangian density for that case:

ansatzLagrangian/.First[Solve[


The result is that a Lagrangian density for OPs equation, but without the term $$u^2u_r/r$$, is $$-\frac{r}{2} e^{u^2} \left( \frac{1}{r^2} + u_r^2 + u_z^2 \right)$$ This density is relative to $$dr\,dz$$, since $$r$$ and $$z$$ are our independent variables.

Comment 1. One can use the Variational Methods Package to calculate the Euler-Lagrange equations corresponding to the Lagrangian density above:

<<VariationalMethods;
With[{u=u[r,z]},
EulerEquations[
-r/2*Exp[u^2]*(1/r^2+D[u,r]^2+D[u,z]^2),
u,{r,z}]
]//Factor


Comment 2. We approached this with an ansatz. One can play with the various parameters above and try to improve. One may go about this in other ways, for example study "integrability conditions" for Euler-Lagrange equations.

Addendum. I thought I try a second approach, this time not based on an ansatz. I will try to state this addendum more or less independently from the answer above.

I want to see if we can find an auxiliary function $$M = M(r,z,u,u_r)$$, and a Lagrangian density of the form $$e^M L(r,z,u,u_r,u_z)$$, such that the associated Euler Lagrange equation is as follows: $$r e^{M(r,z,u,u_r)} \left( u_{rr} + u_{zz} + u_r/r - c_1 u/r^2 + u(c_2 u_r^2 + c_3 u_z^2 + c_4 uu_r/r) \right) = 0$$ Here $$c_1,c_2,c_3,c_4$$ are constants. OP wants $$c_1=c_2=c_3=c_4=1$$ ultimately. This addendum again comes to the conclusion (more convincingly than the answer above I would say) that this only works if $$c_2=c_3$$ and $$c_4=0$$.

Note that also this treatment is based on some (mild) assumptions. For example, I assume from the outset that $$M$$ does not depend on $$u_z$$ explicitly. Note in the displayed equation that the factor $$r$$ on the left is just for convenience, we could absorb it in $$M$$ but prefer to factor it out explicitly. The factor $$e^M$$ in the Lagrangian density could also be absorbed in $$L$$, but we prefer not to for convenience.

The following code implements this (there is some overlap with the code above, but this code is standalone):

independentVars={r,z};
dependentVars={u[r,z]};
On[Assert];
calculateConditions[replacements_,extrapolyvars_:{}]:=Module[{multiplier,lagrangian,desiredEquation,EulerLagrangeEquations,zeroGoal,polyvars},
multiplier=Exp[M[r,z,u[r,z],D[u[r,z],r]]]//.replacements;
lagrangian=multiplier*L[r,z,u[r,z],D[u[r,z],r],D[u[r,z],z]]//.replacements;
desiredEquation=r*multiplier*(D[u[r,z],{r,2}]+D[u[r,z],{z,2}]+D[u[r,z],r]/r-c[1]*u[r,z]/r^2+u[r,z]*(c[2]*D[u[r,z],r]^2+c[3]*D[u[r,z],z]^2+c[4]*u[r,z]*D[u[r,z],r]/r));
EulerLagrangeEquations=Map[Function[{dv},D[lagrangian,dv]-Total[Map[D[D[lagrangian,D[dv,#]],#]&,independentVars]]],dependentVars]//Expand;
zeroGoal=(First[EulerLagrangeEquations]-desiredEquation)/multiplier//Expand;
zeroGoal=zeroGoal/.{u[r,z]->u,D[u[r,z],r]->ur,D[u[r,z],z]->uz,
D[u[r,z],{r,2}]->urr,D[u[r,z],{z,2}]->uzz,D[u[r,z],r,z]->urz};
polyvars=Join[{urr,urz,uzz},extrapolyvars];
Assert[PolynomialQ[zeroGoal,polyvars]];
SortBy[Map[Last,CoefficientRules[zeroGoal,polyvars]],ByteCount]];


What this code does is tell us what differential equations $$M$$ and $$L$$ have to satisfy.

rep={};
First[calculateConditions[rep]]
(* -r-(L^(0,0,0,0,2))[r,z,u,ur,uz] *)


This tells us that the second derivative of $$L$$ with respect to $$u_z$$ must be $$-r$$, which implies that $$L$$ can be written as follows with some new functions $$L_1$$ and $$L_2$$:

AppendTo[rep,L->Function[{r,z,u,ur,uz},L1[r,z,u,ur]+L2[r,z,u,ur]*uz-r*uz^2/2]];


Given this, the next condition is

First[calculateConditions[rep,{uz}]]
(* 2 r (M^(0,0,0,1))[r,z,u,ur] *)


This says that $$M$$ cannot depend on $$u_r$$, so let us add that:

AppendTo[rep,M->Function[{r,z,u,ur},M1[r,z,u]]];


Given this, the next condition is

First[calculateConditions[rep,{uz}]]
(* -2 (L2^(0,0,0,1))[r,z,u,ur] *)


This says that $$L_2$$ cannot depend on $$u_r$$, so let us add that:

AppendTo[rep,L2->Function[{r,z,u,ur},L2a[r,z,u]]];


Given this, the next condition is

First[calculateConditions[rep,{uz}]]
(* r (M1^(0,1,0))[r,z,u] *)


This says that $$M_1$$ cannot depend on $$z$$, so let us add that:

AppendTo[rep,M1->Function[{r,z,u},M1a[r,u]]];


Given this, the next condition is

First[calculateConditions[rep,{uz}]]
(* -r-(L1^(0,0,0,2))[r,z,u,ur] *)


This tells us that the second derivative of $$L_1$$ with respect to $$u_r$$ must be $$-r$$, which implies that $$L$$ can be written as follows with some new functions $$L_{1a}$$ and $$L_{1b}$$:

AppendTo[rep,L1->Function[{r,z,u,ur},L1a[r,z,u]+L1b[r,z,u]*ur-r*ur^2/2]];


Given this, the next condition is

First[calculateConditions[rep,{uz,ur}]]
(* -r u c[2]+1/2 r (M1a^(0,1))[r,u] *)


This tells us how $$M_{1a}$$ depends on $$u$$, so let us add that:

AppendTo[rep,M1a->Function[{r,u},M1b[r]+c[2]*u^2]];


Given this, let us see all remaining equations (there are three):

calculateConditions[rep,{uz,ur}]//Simplify//TableForm
(*
r u (c[2]-c[3])
-u^2 c[4]+r (M1b^\[Prime])[r]
(u c[1])/r+2 u c[2] L1a[r,z,u]-L1b[r,z,u] (M1b^\[Prime])[r]+(L1a^(0,0,1))[r,z,u]-(L2a^(0,1,0))[r,z,u]-(L1b^(1,0,0))[r,z,u]
*)
`

The first forces $$c_2 = c_3$$. The second forces $$c_4 = 0$$, because $$M_{1b}$$ is only a function of $$r$$. It also forces $$M_{1b}$$ to be a constant. The last equation says that there is not restriction on $$c_1$$. So also with this more general calculation, we recovered the result from the answer above, and the answer to OPs original question is still negative.

• @Bobster I used the Lagrangian formalism with $r$ and $z$ as independent variables. The Lagrangian itself (as opposed to just the density) is obtained by multiplying with $drdz$ and integrating formally, which gives $-\frac{1}{2} \int e^{u^2}(1/r^2+u_r^2+u_z^2)\, rdrdz$. I moved the factor $r$ over to $rdrdz$ now, which we recognize as the volume element for polar (or cylindrical) coordinates. Commented Aug 1, 2022 at 4:39
• Can PossibleLagrangian be used for this PDE? I was not able to figure out how to do that. Commented Aug 1, 2022 at 18:05
• @Bobster I tried equations similar to yours, but with only $r$ as independent variable, on Wolfram Alpha. Example 1: Input "u''(r)+2u'(r)-u(r)=0" gives among other things the possible Lagrangian "ℒ(u', u, r) = 1/2 (e^(2 r) u^2 + e^(2 r) (u')^2)". Example 2: Input "u''(r)+u'(r)/r-u(r)/r^2=0" gives possible Lagrangian "ℒ(u', u, r) = 1/2 (u^2/r + r (u')^2)". Example 3: Input "u''(r)+u'(r)/r-u(r)/r^2+u(r)*(u'(r))^2=0" gives no Lagrangian. It is possible that this functionality is only available for linear equations. Commented Aug 2, 2022 at 5:24
• I found a very interesting fact about the existence of Lagrangians for PDEs or systems of PDEs in P. J. Olver's book "Applications of Lie Groups to DEs," Second Edition, page 333. It says that the Frechet derivative of the PDE (or system of PDEs) must be self adjoint for the Euler-Lagrange equations to exist. Also there are Mathematica programs for exactly this (and much more) in Gerd Baumann's book "Symmetry Analysis of DEs with Mathematica" (with CD, 1998). I need some time to look into this. Thank you for all of the work that you have done on this problem. Commented Aug 3, 2022 at 16:11
• @Bobster Those are nice references, thank you. I did not know about these Gerd Baumann books! I looked at Olver briefly. Note that the theorem that you refer to is for Euler-Lagrange EXPRESSIONS. There is an important detail here: Given a differential equation $\Delta = 0$, then even if $\Delta$ is not an Euler-Lagrange expression, it is still possible that some (nonzero) function times $\Delta$ is an Euler-Lagrange expression. Olver discusses this, roughly a page after the theorem, where he says that these Helmholtz conditions are "somewhat unsatisfactory" for these reasons. Commented Aug 3, 2022 at 18:01