Partial derivative operator from 2 variables multiplied by functions

Question

My main issue when I try to put all the equations in Mathematica ,I get different solution from what I expected should be,Epsilon 11,22,12 are functions and multiplied by Their partial derivative operator,Can anyone confirm if the written code is right?

Clear[op1, op2, op3]
Subscript[\[Epsilon], 11] = Subscript[\[Epsilon], rr]*Cos[\[Theta]]^2 + Subscript[\[Epsilon], \[Theta]\[Theta]]*Sin[\[Theta]]^2 - 
     Subscript[\[Epsilon], r\[Theta]]*Sin[2*\[Theta]] & ; 
op1 = Sin[\[Theta]]^2*D[#1, {r, 2}] + Cos[\[Theta]]^2*(1/r)*D[#1, r] + Cos[\[Theta]]^2*(1/r^2)*D[#1, {\[Theta], 2}] - 
     Sin[2*\[Theta]]*(1/r^2)*D[#1, \[Theta]] + Sin[2*\[Theta]]*(1/r)*D[#1, r]*D[#1, \[Theta]] & ; 
firstterm := TrigExpand[op1[Subscript[\[Epsilon], 11][f[r, \[Theta]]]] + 
     op1[f[r, \[Theta]]]*Subscript[\[Epsilon], 11][f[r, \[Theta]]]]; 
Subscript[\[Epsilon], 22] = Subscript[\[Epsilon], rr]*Sin[\[Theta]]^2 + Subscript[\[Epsilon], \[Theta]\[Theta]]*Cos[\[Theta]]^2 + 
     Subscript[\[Epsilon], r\[Theta]]*Sin[2*\[Theta]] & ; 
op2 = Cos[\[Theta]]^2*D[#1, {r, 2}] + Sin[\[Theta]]^2*(1/r)*D[#1, r] + Sin[\[Theta]]^2*(1/r^2)*D[#1, {\[Theta], 2}] + 
     Sin[2*\[Theta]]*(1/r^2)*D[#1, \[Theta]] - Sin[2*\[Theta]]*(1/r)*D[#1, r]*D[#1, \[Theta]] & ; 
secondterm := TrigExpand[op2[f[r, \[Theta]]]*Subscript[\[Epsilon], 22][f[r, \[Theta]]] + 
     op2[Subscript[\[Epsilon], 22][f[r, \[Theta]]]]]; 
SSSS = Expand[firstterm + secondterm]; 
Subscript[\[Epsilon], 12] = -0.5*Subscript[\[Epsilon], \[Theta]\[Theta]]*Sin[2*\[Theta]] + 0.5*Subscript[\[Epsilon], rr]*Sin[2*\[Theta]] + 
     Subscript[\[Epsilon], r\[Theta]]*Cos[2*\[Theta]] & ; 
op3 = 2*(Sin[\[Theta]]*Cos[\[Theta]]*D[#1, {r, 2}] - Sin[\[Theta]]*Cos[\[Theta]]*(1/r)*D[#1, r] - 
      Sin[\[Theta]]*Cos[\[Theta]]*(1/r^2)*D[#1, {\[Theta], 2}] - Cos[2*\[Theta]]*(1/r^2)*D[#1, \[Theta]] + 
      Cos[2*\[Theta]]*(1/r)*D[#1, r]*D[#1, \[Theta]]) & ; 
thirdterm = Expand[op3[f[r, \[Theta]]]*Subscript[\[Epsilon], 12][f[r, \[Theta]]] + op3[Subscript[\[Epsilon], 12][f[r, \[Theta]]]]]; 
TrigExpand[TrigExpand[thirdterm - SSSS]]

The solution that I get in Mathematica is like below

-((2*Subscript[\[Epsilon], r\[Theta]]*
      Derivative[0, 1][f][r, \[Theta]])/r^2) - 
   (1.*Subscript[\[Epsilon], rr]*Derivative[0, 2][f][r, \[Theta]])/
  r^2 - 
   (1.*Subscript[\[Epsilon], rr]*Derivative[1, 0][f][r, \[Theta]])/r + 
   (2*Subscript[\[Epsilon], r\[Theta]]*
    Derivative[0, 1][f][r, \[Theta]]*
    Derivative[1, 0][f][r, \[Theta]])/r - 
   1.*Subscript[\[Epsilon], \[Theta]\[Theta]]*
  Derivative[2, 0][f][r, \[Theta]]

Answer 1

Here is an approach that gives an expression similar to the one you want. The main differences between this approach and the one in the original post are

1. Variables $\epsilon_{r\theta}$, etc, must be explicit functions of $r,\theta$.
2. Variables $\epsilon_{11}$, etc, are not "functions", so no & at the end.
3. The mixed partial derivatives, $\partial^2/\partial r\partial\theta$ are second derivatives, not the product of first derivatives.
4. Subscripts are not used, so the code is easier to read.
5. Exact values, like 1/2, are used instead of approximate ones, like 0.5.

First, here are some preliminary definitions. The rules allow us to work with $f, g, h$ instead of the subscripted variables. The next definition tells Mathematica how to display the partial derivatives when we use TraditionalForm. These definitions will be used later.

ClearAll["Global`*"]
Remove[ϵ]
rules = {
   f -> Subscript[ϵ, rr],
   g -> Subscript[ϵ, θθ],
   h -> Subscript[ϵ, rθ]};
Derivative /: 
 MakeBoxes[Derivative[α__][f1_][vars__?AtomQ], 
  TraditionalForm] := 
 Module[{bb, dd, sp}, 
  MakeBoxes[dd, _] ^= 
   If[Length[{α}] == 1, "[DifferentialD]", "∂"];
  MakeBoxes[sp, _] ^= "[ThinSpace]";
  bb /: MakeBoxes[bb[x__], _] := RowBox[Map[ToBoxes[#] &, {x}]];
  TemplateBox[{ToBoxes[bb[dd^Plus[α], f1]], 
    ToBoxes[Apply[bb, 
      Riffle[Map[
        bb[dd, #] &, (Pick[{vars}, #]^Pick[{α}, #] &[
          Thread[{α} > 0]])], sp]]], 
    ToBoxes[Derivative[α][f1][vars]]}, "ShortFraction", 
   DisplayFunction :> (FractionBox[#1, #2] &), 
   InterpretationFunction :> (#3 &), Tooltip -> Automatic]]

Now, here is the main calculation. The variables $f, g,h$ are explicit functions of $r,\theta$.

ϵ11 = f[r, θ]*Cos[θ]^2 + g[r, θ]*Sin[θ]^2 - h[r, θ]*Sin[2*θ];
ϵ22 = f[r, θ]Sin[θ]^2 + g[r, θ]Cos[θ]^2 + h[r, θ]Sin[2θ];
ϵ12 = 1/2f[r, θ]Sin[2*θ] - 1/2g[r, θ]Sin[2*θ] + h[r, θ]Cos[2θ];
op1 = Sin[θ]^2D[#1, {r, 2}] +
    Cos[θ]^2 (1/r)D[#1, r] + 
    Cos[θ]^2(1/r^2)D[#1, {θ, 2}] -
    Sin[2*θ](1/r^2)D[#1, θ] + 
    Sin[2*θ](1/r)D[#1, r, θ] &;
op2 = Cos[θ]^2D[#1, {r, 2}] +
    Sin[θ]^2(1/r)D[#1, r] + 
    Sin[θ]^2(1/r^2)D[#1, {θ, 2}] + 
    Sin[2θ](1/r^2)D[#1, θ] - 
    Sin[2*θ](1/r)D[#1, r, θ] &;
op3 = 2(Sin[θ]Cos[θ]D[#1, {r, 2}] -
      Sin[θ]Cos[θ](1/r)D[#1, r] -
      Sin[θ]Cos[θ](1/r^2)D[#1, {θ, 2}] -
      Cos[2θ](1/r^2)D[#1, θ] + 
      Cos[2*θ](1/r)D[#1, r, θ]) &;
lhs = op1[ϵ11] + op2[ϵ22];
rhs = op3[ϵ12];
FullSimplify[rhs - lhs] /. rules // Expand // TraditionalForm

This expression is, of course, the left hand side of the desired equation minus the right hand side. The expression has the same 6 derivatives as the desired equation. The signs and the exponents of $r$ are all the same. The only difference that I see is the first two terms of this expression have the factor of 2 in the numerator instead of in the denominator. Perhaps there is a typo somewhere?