# Clear Guidelines for Symbolic Calculation With Units?

The following notebook performs a relatively simple physical calculation interacting position, velocity and acceleration of two charged point particles according to Weber electrodynamics.It is simple enough that the correctness of the physical dimensions is apparent by inspection.

Having read the Symbolic Calculations with Units section on Differentiation it wasn't clear to me what restrictions might apply during evaluation. This resulted in a number of failed attempts to get a dynamical system formulated.

So what guidelines am I violating here?

QuantityVector[name_String, components_List, unit_, vars_List : {}] :=
Module[{syms},
If[Length[vars] > 0,
syms = Quantity[Symbol[name <> #][Sequence @@ vars], unit] & /@
components,
syms = Quantity[Symbol[name <> #], unit] & /@ components];
syms]
Q = Quantity[QM, "Coulombs"];
q = Quantity[qM, "Coulombs"];
Txhertz = Quantity[TxhertzM, "Seconds^-1"];
Rxhertz = Quantity[RxhertzM, "Seconds^-1"];
rVec = QuantityVector["Txr", {"x", "y", "z"}, "Meters", {t}] -
QuantityVector["Rxr", {"x", "y", "z"}, "Meters", {t}]
vVec = D[rVec, t]
aVec = D[vVec, t]
rUnitVec = rVec/Norm[rVec]

c = Quantity["SpeedOfLight"];(*Speed of light in m/s*)

epsilon0 = Quantity["ElectricConstant"];
weberForce[
t_] = (Q q rUnitVec)/(4 Pi epsilon0 Norm[rVec]^2)*(1 + (1/
c^2) (vVec . vVec - 3/2 ((rUnitVec) . vVec)^2 +
Norm[rVec] rUnitVec . aVec));
(weberForce[t] /. {
Txrx -> Function[0 #],
Txry -> Function[0 #],
Txrz -> Function[Quantity[100, "Meters"] + Sin[# Txhertz]]
} /.
{
Rxrx -> Function[0 #],
Rxry -> Function[0 #],
Rxrz -> Function[Sin[# Rxhertz]]
} /. {t -> Quantity[1/10^9, "Second"], TxhertzM -> 1.3*10^9,
RxhertzM -> 1.3*10^8, qM -> 1, QM -> 1})[[3]]


Result: $$\frac{(100\text{m}+0.833924) \left(\frac{1}{4 \pi | 100\text{m}+0.833924| ^2}\text{C}^2\;\text{/(}\text{m}^2\varepsilon _0)\right) \left(\left(1/\text{c}^2\right) \left(4.78933\times 10^{16}\text{m}^2 \; /\text{s}^2+\frac{1}{2} (-3) \left(\frac{(100\text{m}+0.833924) \left(2.18845\times 10^8\text{m}/\text{s}\right)}{| 100\text{m}+0.833924| }+0\text{m}\right)^2+(| 100\text{m}+0.833924| \text{m}) \left(\frac{(100\text{m}+0.833924) \left(-1.62622\times 10^{18}\text{m}/\text{s}^2\right)}{| 100\text{m}+0.833924| }+0\text{m}\right)\right)+1\right)}{| 100\text{m}+0.833924| }$$

• (1) As per the instructions, please don't use [bugs] tag on new questions. (2) You probably forgot to include some part of the code (definition of Q and q). (3) Not sure if it is the same, but there is a built-in QuantityArray which you can use instead of QuantityVector. (4) You have to let D know that you are differentiating with respect to a Quantity, for example: D[Quantity[t, "Meters"], Quantity[t, "Seconds"]]. Commented Mar 25 at 16:05
• My apologies for failing to follow the instructions. My omissions in this post (now filled in after editing) aren't relevant to the actual result because what I forgot was to include the entirety of the original notebook, which did, indeed, have the definitions of all parameters, including t as a t -> Quantity[1/10^9, "Second"]. I further modified the notebook to head off some additional speculations as to the cause of the dropping of dimensions during differentiation. Commented Mar 25 at 18:39
• QuantityArray doesn't support the generation of dimensioned Functions of t for {x,y,z} which is necessary in any dynamical system model. Commented Mar 25 at 18:52
• I believe you have a mistake in your replacement function. You are replacing Txrx with 0, but Txrz with 100 m + Sin[...]. This is incompatible. You probably want to do Txrz -> Function[100 + Sin[# Txhertz]]. Then, rUnitVec will come out correctly as {0, 0, 1.}. The second problem is the way you represent time in your equations. Can you please look at Symbolic Calculations with Units? Commented Mar 26 at 10:24
• Yes, I recognized that this morning -- which is an error I introduced after having encountered incompatible units and incorrectly "correcting" it. The problem probably has something to do with your second comment regarding time representation. I did, originally, look at Symbolic Calculation with Units but apparently didn't comprehend it so I'll go back and do some homework on that. Commented Mar 26 at 15:45

A point of confusion contributing to my difficulties was the distinction between variables that are assigned Quantity values and variables that are assigned dimensionless QuantityMagnitude values. For instance, in the Integration section of the Symbolic Calculations with Units it says:

If the variable of integration is a Quantity object, then its QuantityMagnitude is used as the integration variable

In point of fact, when they give various examples of "symbols" they usually refer to dimensionless QuantityMagnitude values rather than dimensioned Quantity values.

Although the General section does show how one can use substitution to do dimensional checking of an expression involving symbols, the resulting output expression refers only to the symbols for the (dimensionless) QuantityMagnitude. Although I did look for some way of defining symbols to be constrained to particular dimensions, which one could then use in expressions that would automatically check for commensurability during additive operations and do automatic units conversions as needed , QuantityVariable, QuantityVariableDimensions, etc. did not provide this.

So below is a working script that does everything with QuantityMagnitude expressions using a naming convention where the variable names end in M to avoid running into trouble:

QuantityVector[name_String, components_List, unit_, vars_List : {}] :=
Module[{syms},
If[Length[vars] > 0,
syms = Quantity[Symbol[name <> #][Sequence @@ vars], unit] & /@
components,
syms = Quantity[Symbol[name <> #], unit] & /@ components];
syms]
Q = Quantity[QM, "Coulombs"];
q = Quantity[qM, "Coulombs"];
(*Txhertz=Quantity[TxhertzM,"Seconds^-1"];
Rxhertz=Quantity[RxhertzM,"Seconds^-1"];*)

tQuantity = Quantity[tM, "Seconds"];
subs1 = {
Txrx -> Function[0],
Txry -> Function[0],
Txrz -> Function[100 + TxAntennaLengthM*Sin[# 2 Pi TxhertzM] ],
Rxrx -> Function[0],
Rxry -> Function[0],
Rxrz -> Function[RxAntennaLengthM Sin[# 2 Pi RxhertzM]]};
subs1

subs2 = {TxhertzM -> 1.3*10^9, RxhertzM -> 1.3*10^8, qM -> 1, QM -> 1,
TxAntennaLengthM -> .00001, RxAntennaLengthM -> .001};
subs3 = {tM -> 1/10^9};
(*rVec=QuantityVector["r",{"x","y","z"},"Meters"]*)
(* Transmitter \
and Receiver (Tx Rx) relative position as a function of time *)

Txr = QuantityVector["Txr", {"x", "y", "z"}, "Meters", {tM}];
Rxr = QuantityVector["Rxr", {"x", "y", "z"}, "Meters", {tM}];
vRx = D[Rxr,
Quantity[tM,
"Seconds"]] (* velocity of the receiving antenna's charge *)

rVec = Txr - Rxr;
vVec = D[rVec, tQuantity];
aVec = D[vVec, tQuantity];
c = Quantity["SpeedOfLight"];(*Speed of light in m/s*)
epsilon0 =
Quantity["ElectricConstant"];
rUnitVec = rVec/Norm[rVec];
DVecs = (vVec . vVec - 3/2 ((rUnitVec) . vVec)^2 +
Norm[rVec] rUnitVec .
aVec); (* m^2/s^2 *)
DVecsDimensionless = (1 + DVecs/c^2);
weberForce[tM_] = (Q q rUnitVec)/(4 Pi epsilon0 Norm[rVec]^2)*
DVecsDimensionless;

(* Note carefully that *xhertz is previously defined in terms of \
*xhertM.  Therefore *xhertzM is what is being substituted. *)

plotWeber[tM_] = UnitConvert[weberForce[tM] /. subs1 /. subs2];
Plot[plotWeber[tM][[3]], {tM, 0, 1/10^8}]
Plot[vRx[[3]] /. subs1 /. subs2, {tM, 0, 1/10^8}]