I’m working on a problem involving Kelvin-Helmholtz instability and I’m trying to perform the linearization on Mathematica. For some reason when I get to my first order equation (ε × stuff + O(ε²) == 0) and I ask Mathematica to divide by epsilon, it refuses to carry the operation through. Here is my code:

Remove["Global`*"]; $Assumptions =.;
$Assumptions = {c > 0};

xx := {x, y, z} (*position vector*);
u := {ux[xx, t], uy[xx, t], uz[xx, t]} (*perturbation velocity*);
ρ1 := ϱ1[xx, t] (*perturbation density*);
P1 := p1[xx, t] (*perturbation pressure*);

(*Specify background flow:*)
v := {0, 0, vz[x]};
P0 := p0;
ρ0 := ϱ0[x];

(*Quantities of problem*)
SetAttributes[eps, Constant]
ρ := ρ0 + eps* ρ1;
P := P0 + eps* P1;
vv := v + eps*u;

(*Lagrangian (comoving) derivative*)
DDt[function_] := D[function,t] + Sum[ vv[[i]] D[#,xx[[i]]], {i,3}]& @(function)

(*Continuity, momentum and energy:*)
eqsOfMotion := {
   DDt[ρ] == -ρ Div[vv,xx],
   ρ DDt[vv] == -Grad[P,xx],
   DDt[P] == c^2 DDt[ρ]
   } ;

  eqsOfMotion /. eps -> 0], "-> background flow is valid solution"]

(*Since background flow is valid solution 
we have that every term in eqsOfMotion is 
proportional to eps, hence can divide through:*)

And that last line is where my problem arises. The 1/eps term refuses to pass through the parentheses and cancel with eps. Any ideas?

My version of Mathematica is 11.3.

  • 1
    $\begingroup$ See DivideSides in documentation, looks to have been introduced in 11.3. Equations do not automatically absorb arithmetic operations. $\endgroup$ – eyorble Oct 29 '20 at 18:37
  • $\begingroup$ Or if DivideSides is not available: eqsOfMotion2 = eqsOfMotion /. Equal[lhs_, rhs_] :> Equal[lhs/eps, rhs/eps] // Simplify $\endgroup$ – Bob Hanlon Oct 29 '20 at 19:04
  • $\begingroup$ @eyorble DivideSides does not quite work for me, as it adds in too much complexity, but by leafing through the documentation for it I found a method that works with ApplySides — thank you! $\endgroup$ – confusedandbemused Oct 30 '20 at 8:17

Following @eyorble's suggestion to look at DivideSides, I found that a method that does the trick is

ApplySides[#/eps&, eqsOfMotion]

(and apparently another valid solution replaces ApplySides with Map, but I couldn’t get it to work).

DivideSides does not work for me because it is too intelligent for its own good, and creates a proliferation of cases in which I have no interest.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.