Re: Calculating with Units
- To: mathgroup at smc.vnet.net
- Subject: [mg94090] Re: [mg93997] Calculating with Units
- From: Daniel Lichtblau <danl at wolfram.com>
- Date: Thu, 4 Dec 2008 07:17:45 -0500 (EST)
- References: <200812011202.HAA05923@smc.vnet.net> <493467F2.3020709@wolfram.com>
Nikolaus Rath wrote:
> Hello,
>
> Until now, I have always defined my units manually in terms of
> the 4 base units Kilogram, Meter, Coulomb and Second:
>
> Newton = Kg Meter / Second^2;
> Joule = Newton Meter;
> Henry = Weber / Ampere;
> Weber = Volt Second;
> etc.
>
> This allows me to evaluate expressions like
>
> B = (\[Mu] Ii)/(2 \[Pi] r);
>
> B /. { r -> 20 Meter, \[Mu] -> 4 \[Pi] * 10^-7 Henry/Meter,
> Ii -> 2 Ampere}
>
> % / Tesla
>
> and directly see if the units in the final result are correct.
>
> But since all these units (and many more) are already defined in the
> Mathematica Units package, I was wondering if there isn't a way to use
> that package in the same way. Unfortunately, it doesn't seem to
> simplify the units, i.e. when I enter
>
> (Ampere Henry)/(Meter Meter) // Simplify // SI
>
> I get the same expression back. What I would like to get is either
> "Tesla" or "Kilo Gram / (Coulomb Second)". Is there a way to
> accomplish this?
>
>
> Best,
>
>
> -Nikolaus
> [...]
My first reply, in private email, was to the effect that GroebnerBasis
and PolynomialReduce could be used. This was correct, but after further
prodding I realized it needs considerable elaboration.
In looking over some slides from a 2007 talk I realized I had more or
less done this once before, and the fact that a similar question has
appeared on MathGroup suggest that it is not a completely obscure thing
to want. It might be almost completely obscure (ACO), but all the same
I'm using this as an excuse to post a response.
We start with a list of units and a generating set of defining
relations. I also use some one and two letter abbreviations for brevity.
units = {Farads, Ohms, Pascals, Watts, Webers, Teslas, Henries,
Newtons, Joules, Volts, Amperes, Coulombs, Seconds,
Kilograms, Meters};
{f, o, p, wa, we, t, h, n, j, v, a, c, s, k, m} = units;
relations = {f - (c^2*s^2)/(k*m^2), o - (k*m^2)/(c^2*s),
p - k/(m*s^2), -((k*m^2)/s^3) + wa, -((k*m^2)/(c*s)) + we,
-(k/(c*s)) + t, h - (k*m^2)/c^2, n - (k*m)/s^2,
j - (k*m^2)/s^2, -((k*m^2)/(c*s^2)) + v, a - c/s};
Put each over a common denominator and extract the numerator. This gives
polynomial (rather than rational function) relations. This will be
important because Groebnerbasis and PolynomialReduce behave best when
handling legitimate polynomial input.
polyrelations = Numerator[Together[relations]]
Here is a problem. We want to sensibly handle denominators. That is
solved by creating the reciprocal variables and corresponding relations.
recipunits = Map[recip, units];
reciprels = Map[#*recip[#] - 1 &, units];
allrelations = Join[Numerator[Together[relations]], reciprels];
allvars = Flatten[Transpose[{recipunits, units}]];
Since the question involves representing units monomials in terms of a
core set of charge, time, mass, and length (a typical scenario), we
ordered those last. Also we ordered their corresponding reciprocals
after other variables.
With all this we can compute a Groebner basis, using the ordered
variables to define the (lexicographic) term order.
There is a catch here. The input defines what is called a toric ideal.
GroebnerBasis has reasonably fast special case code for that. It also
has code that attempts to recognize such input. But that recognition
code is a bit simple and fails in this example. And the general case
code is quite slow in this case (several minutes thus far, and still no
joy). So we resort to an explicit invocation of the special case code.
gbt = GroebnerBasis`ToricGroebnerBasis[allrelations, allvars];
Now we look at the given example. I've pluralized the units to
correspond with my version.
B = (\[Mu]*Ii)/(2*Pi*r);
B2 = B /.
{r->20*Meters, \[Mu]->(4*Pi*(Henries/Meters))/10^7, Ii->2*Amperes}
Out[29]= (Amperes*Henries)/(50000000*Meters^2)
We want to show that B2/Teslas is dimensionless, that is, it is a
numerical constant. We start by converting to an explicit polynomial.
This is done by a replacement rule that changes units with negative
exponents into their reciprocals, with corresponding positive exponents.
We then use PolynomialReduce with respect to gbt (keeping the same
variable ordering as from computing gbt). Finally we convert explicit
reciprocal variables to units with negative exponents.
Last[PolynomialReduce[
B2/Teslas /. (z_)^(n_Integer) /; n < 0 :> recip[z]^(-n),
gbt, allvars]] /. recip[z_] -> 1/z
Out[33]//InputForm= 1/50000000
This did give a constant, as we wanted.
To see what are the "fundamental" units of B2, we could simply do
In[36]:= InputForm[ Last[PolynomialReduce[
B2 /. (z_)^(n_Integer) /; n < 0 :> recip[z]^(-n),
gbt, allvars]] /. recip[z_] -> 1/z]
Out[36]//InputForm= Kilograms/(50000000*Coulombs*Seconds)
This, and other methods, can can be used to rewrite units monomials in
terms of equivalent ones that minimize total degree, or minimize the
maximal degree of numerator vs denominator, or minimize the number of
units that appears (possibly making the degree large in this scenario).
The starting point for all methods of which I am aware is to augment
with reciprocal variables and relations.
Daniel Lichtblau
Wolfram Research
- References:
- Calculating with Units
- From: Nikolaus Rath <Nikolaus@rath.org>
- Calculating with Units