Simplifying expressions

Question

Let's define a function $\Omega(x,y)$

r1 = Sqrt[(x - 0.5)^2 + y^2];
r2 = Sqrt[(x + 0.5)^2 + y^2];
p1 = 1/r1^3 + λ/r2^3;
q1 = 1/2*(1/r1^3 - λ/r2^3);
Ax = -y*p1;
Ay = x*p1 - q1;
Ω = (x*Ay - y*Ax) + 1/2*(x^2 + y^2);

Now I want to obtain the first and the second derivatives of $\Omega$ at the most simplified expressions. After conducting the calculations by hand I derived

$\Omega_x = x + 2xp_1 - q1 + 3y^2q2 - 3x(x^2 + y^2)p_2 + 6x^2q_2 - 3xs_2$

$\Omega_y = y + 2yp_1 + 3xyq_2 - 3y(x^2 + y^2)p_2$

where

p2 = (1/r1^5 + λ/r2^5);
q2 = 1/2*(1/r1^5 - λ/r2^5);
s2 = 1/4*(1/r1^5 + λ/r2^5);

I tried to obtain the same results using Mathematica with

Ωx = D[Ω, x];
Ωy = D[Ω, y];

rule = {Sqrt[(x - 0.5)^2 + y^2] -> r1, 
Sqrt[(x + 0.5)^2 + y^2] -> r2, 
(1/r1^3 + λ/r2^3) -> p1, 1/2*(1/r1^3 - λ/r2^3) -> 
q1, (1/r1^5 + λ/r2^5) -> p2, 1/2*(1/r1^5 - λ/r2^5) -> 
q2, 1/4*(1/r1^5 + λ/r2^5) -> s2};

FullSimplify[Ωx] /. rule
FullSimplify[Ωy] /. rule

Unfortunately the program fails to provide the same elegant expressions obtained by hand. Is there a way to achieve this with Mathematica? For the derivatives of the second order the calculations are much more complicated and therefore I need the program to compute them in simplified form.

Ωxx = D[Ωx, x];
Ωxy = D[Ωx, y];
Ωyy = D[Ωy, y];

So, my question is how to obtain the expressions of the three second order derivatives in a simplified form as the ones, derived by hand regarding the first order derivatives?

Answer 1

Pattern matching of Mathematica isn't semantic but syntactic, so your replacement rules won't apply to the expression unless the FullForm of them matches. For more information you may want to read this. (Sadly the package given in that page isn't strong enough yet and hasn't been updated for years.)

Since what you want is to simplify the second order derivatives, I modified the rule to make it more general:

Clear[r1, r2, p1, p2, q1, q2]
rule = {Sqrt[(x - 0.5)^2 + y^2] -> r1, 
   Sqrt[(x + 0.5)^2 + y^2] -> r2, (1/r1^n + λ/r2^n) -> p[n], 
   1/2*(1/r1^n - λ/r2^n) -> q[n]};

patt = r1 | r2;
rulemid = (a : patt)^n_ :> a[n];

replace = Simplify[# /. MapAt[#^2 &, rule, {1 ;; 2, All}] /. 
       s : (patt^2)^Rational[_, 2] :> Simplify[s, r1 > 0 && r2 > 0] /. 
      rulemid /. (First@
           Solve[Equal @@@ #, Union@Cases[#, patt[_], Infinity]] &@(rule[[3 ;; 4]] /. 
           rulemid) /. n -> -n /. (a : patt)[n] :> a[n_])] &;


replace@Ωx
(* -3. (0. + x^3 p[5] + x (-0.333333 - 0.666667 p[3] + 
   (0.25 + y^2) p[5]) + 0.333333 q[3] - 2. x^2 q[5] - 1. y^2 q[5]) *)
replace@Ωy
(* y (1 + 2 p[3] - 3 x^2 p[5] - 3 y^2 p[5] + 3 x q[5]) *)