Mathematica Asked by confusedandbemused on November 16, 2020
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[ρ]
} ;
Print[Simplify[
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:*)
eqsOfMotion/eps
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.
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.
Answered by confusedandbemused on November 16, 2020
Get help from others!
Recent Questions
Recent Answers
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP