I may have a solution but it is not slick! See below

Doing numerical calculations with functions that use units can be very slow. I would like to have a function or operator that can turn my functions into unitless functions. I have a partial solution but it doesn't work in all situations.

The blackbody radiation equation:

BB[l_Quantity, T_Quantity] := (
  2 Quantity["PlanckConstant"] Quantity[
    "SpeedOfLight"]^2)/(l^5 Quantity[1, "Steradians"]) 1/(
      "PlanckConstant"] Quantity[
       "SpeedOfLight"]/(l Quantity["BoltzmannConstant"] T)] - 

Now if I want to do many calculations with this equation it is very slow as can be seen from just making a plot:

     BB[Quantity[l, "Micrometers"], Quantity[333, "Kelvins"]],
      "Microflicks"]], {l, .1, 15}
   , Frame -> True
   , FrameLabel -> {"\[Lambda] (\[Mu]m)", 
     "\[ScriptCapitalL] (\[Mu]flicks)"}
   , ImageSize -> Medium

Slow Plot of Blackbody spectral radiation

Please keep in mind that I want to do much more than just make plots quickly. I need to do Monte Carlo simulations where I will make compiled versions of my equations. Also in my real work, the equations are much more complicated and automatically generated.

The function that I have come up with that works some of the time is the following:

   outunits_] := ((UnitConvert[#, outunits] & /@ 
       fun) /. {UnitConvert[a_, outunits] :>  a}) /. {Quantity[q_, 
      u_] :> q};

I can apply this function to BB[l,T] to get:

BBnu[l_, T_] := 
     T] /. {l -> Quantity[l, "Micrometers"], 
     T -> Quantity[T, "Kelvins"]}, "Microflicks"]]

Where I have to tell my NoUnits[] function what units to use for each variable and the units I want for the final output. The new function BBnu[l,T] is just a function of numbers and gives a numerical output. Now the calculations go much faster!

AbsoluteTiming[Plot[BBnu[l, 333], {l, .1, 15}
  , Frame -> True
  , FrameLabel -> {"\[Lambda] (\[Mu]m)", 
    "\[ScriptCapitalL] (\[Mu]flicks)"}
  , ImageSize -> Medium

fast Blackbody plot

But now if I try my function NoUnits[] on something more complex like a linear combination of blackbodies it fails:

NoUnits[(a BB[l, T1] + b BB[l, T2]) /. {l -> 
    Quantity[l, "Micrometers"], T1 -> Quantity[T1, "Kelvins"], 
   T2 -> Quantity[T2, "Kelvins"]}, "Microflicks"]


I have tried what I have learned from: How to create fast functions based on units? Converting the units of an equation or expression (how to become 3 times heavier…) Yet these do not solve my full problem.

Edit The following are some stressing testing functions:

The blackbody equation again B[l,T]:

BB[l_Quantity, T_Quantity] := (
   2 Quantity["PlanckConstant"] Quantity["SpeedOfLight"]^2)/(
   l^5 Quantity[1, "Steradians"]) 1/(
       "PlanckConstant"] Quantity[
        "SpeedOfLight"]/(l Quantity["BoltzmannConstant"] T)] - 1);

A made up trasnsmsion equation trans[l]:

trans[l_Quantity] := 
  1/(1 + Exp[- 
       Quantity[1, ("Micrometers")^-1] (l - 
        Quantity[4, "Micrometers"])]) 1/(
   1 + Exp[Quantity[3, (
       "Micrometers")^-1] (l - Quantity[10, "Micrometers"])]);

A madeup source equation that includes an interpolation function with compound units Ladd[l,T]:

   fun = Interpolation[{{4/300, 0}, {5/300, 
       1 10^4}, {6/300, .5 10^6}, {8/300, 8 10^6}, {10/300, 
       4 10^6}, {12/300, 2 10^6}, {14/300, 1 10^6}, {15/300, 0}}, 
     InterpolationOrder -> 3]
  Ladd[l_Quantity, T_Quantity] := Piecewise[{
       "Kilograms" ("Meters")^-1 ("Seconds")^-3 ("Steradians")^-1], 
      l/T <= Quantity[4, "Micrometers"]/Quantity[300, "Kelvins"]}
     , {Quantity[
         UnitConvert[l, "Micrometers"]/UnitConvert[T, "Kelvins"]]], 
       "Kilograms" ("Meters")^-1 ("Seconds")^-3 ("Steradians")^-1], 
      Quantity[4, "Micrometers"]/Quantity[300, "Kelvins"] < l/T <= 
       Quantity[15, "Micrometers"]/Quantity[300, "Kelvins"]}
     , {Quantity[0, 
       "Kilograms" ("Meters")^-1 ("Seconds")^-3 ("Steradians")^-1], 
      l/T > Quantity[15, "Micrometers"]/Quantity[300, "Kelvins"]}
     }, Quantity[0, 
     "Kilograms" ("Meters")^-1 ("Seconds")^-3 ("Steradians")^-1]]

Another made up function: LaddTwo[l]

   fun = Interpolation[{{4, 0}, {5, 1 10^4}, {6, .5 10^6}, {8, 
       8 10^6}, {10, 4 10^6}, {12, 2 10^6}, {14, 1 10^6}, {15, 0}}, 
     InterpolationOrder -> 3]
  LaddTwo[l_Quantity] := Piecewise[{
       "Kilograms" ("Meters")^-1 ("Seconds")^-3 ("Steradians")^-1], 
      l <= Quantity[4, "Micrometers"]}
     , {Quantity[
       fun[QuantityMagnitude[UnitConvert[l, "Micrometers"]]], 
       "Kilograms" ("Meters")^-1 ("Seconds")^-3 ("Steradians")^-1], 
      Quantity[4, "Micrometers"] < l <= Quantity[15, "Micrometers"]}
     , {Quantity[0, 
       "Kilograms" ("Meters")^-1 ("Seconds")^-3 ("Steradians")^-1], 
      l > Quantity[15, "Micrometers"]}
     }, Quantity[0, 
     "Kilograms" ("Meters")^-1 ("Seconds")^-3 ("Steradians")^-1]]

and finally, a full system equation that brings it all together Lfull[l,T]:

Lfull[l_Quantity, T_Quantity, a_, b_] := 
  a trans[l] BB[l, T] + b trans[l] Ladd[l, T]+ a b LaddTwo[l];

I'm not trying to make this overly completed but this is the smallest full-featured example.

The goal would be to get this "units to unitless" operator to work on the Lfull[l,T,a,b] function.

My poor solution I'm still very interested in a better solution! I can think of some cases where this may fail. Update: It has not failed yet but I still do not like the implementation. I have updated my function complete with some documentation:

NoUnits::usage = 
  "NoUnits[fun, rp, outunits] Will take a function fun that is \
defined with units and convert it into just a numerical function. Of \
course the user is now responsible for providing the numbers in the \
correct \"Units\".   rp is a list of variables and the corresponding \
units the user wants to use for the resulting function.  The units on 
 the output is the outunits.
  \"MaxPower\" What maximum power to use when looking for match \
units. Default is 3. If you have something like l^5/a^2 in your \
equation you may wat to use 5 for your maximum power.
  fun[v_Quantity,m_Quantity,c_]:= c m \!\(\*SuperscriptBox[\(v\), \(2\
  Kinetic energy equation were we have to put in c=0.5;
  Test it for some inputs:
  Lets say we want just a numbers equation where we input mph, pounds \
and want Ergs out Then do:
  \!\(\*FractionBox[\(88523495162333\\\ c\\\ m\\\ \
\*SuperscriptBox[\(v\), \(2\)]\), \(97656250\)]\)
  Test it:
  eq/.{c\[Rule].5,m\[Rule] 10, v\[Rule] 55}
  Use the original equation:
Options[NoUnits] = {
   "MaxPower" -> Automatic
NoUnits[fun_, rp_List, outunits_, opts : OptionsPattern[NoUnits]] := 
  Module[{eq, list, u, pow, unit, a, b, mp},
   mp = OptionValue["MaxPower"];
   If[Not[NumericQ[mp] && IntegerQ[mp]], mp = 3];
   If[mp <= 0, mp = 3];
   eq = fun /. 
     Evaluate[Rule[#[[1]], Quantity[#[[1]], #[[2]]]] & /@ rp];
   eq = eq /. {Quantity[a_, b_] :> 
        UnitConvert[Quantity[a, b]]} /. {Quantity[a_, 
        "PureUnities"] :> a, Quantity[a_, "DimensionlessUnit"] :> a};
   eq = eq /. {Quantity[a_, b_] /; CompatibleUnitQ[b, outunits] :> 
       QuantityMagnitude[UnitConvert[Quantity[a, b], outunits]]};
   list = 
    DeleteDuplicates[Cases[eq, Quantity[a_, b_] :> b, \[Infinity]]];
   u = rp[[All, 2]];
   pow = Tuples[Range[-mp, mp], Length[u]];
   pow = SortBy[pow, Total[Abs[#]] &];
   pow = {u, #}\[Transpose] & /@ pow;
    unit = 
          Sequence @@ #] & /@ (pow /. {{a_String, b_} :> 
            Power[a, b]}), CompatibleUnitQ[#, ele] &, 1][[1]];
    eq = eq /. {Quantity[a_, b_] /; CompatibleUnitQ[b, unit] :> 
        QuantityMagnitude[UnitConvert[Quantity[a, b], unit]]};
    , {ele, list}];

Here is an example of the usage:


Instead of using UnitConvert, you can just divide the output by the relevant output unit. Here's one implementation of the idea, based on the OP's version:

ClearAll @ NoUnits;
SetAttributes[NoUnits, HoldAll];

NoUnits[fun_, quants : {__Quantity}, output : Quantity[1, out_]] := 
   Expand[fun/output /. Thread[quants[[All, 1]] -> quants]] /. x_ y_Piecewise :> PiecewiseExpand[x y] /. Quantity[x_, _] :> x


LFullNU[l_, T_, a_, b_] = NoUnits2[Lfull[l, T, a, b], {Quantity[l, "Micrometers"], Quantity[T, "Kelvins"]}, Quantity[1, "Microflicks"]];

Seems to work. I have as inputs the full Quantitys to make the replacements easier.

As a check, you can try

BBNU[l_, T_] = NoUnits2[BB[l, T], {Quantity[l, "Micrometers"], Quantity[T, "Kelvins"]}, Quantity[1, "Microflicks"]];

which yields the same result as in the OP.

  • $\begingroup$ I love the elegance of the solution but it didn't work. I did: Clear[LfullNU]; LfullNU[l_, T_, a_, b_] := Evaluate[NoUnits[ Lfull[l, T, a, b], {Quantity[l, "Micrometers"], Quantity[T, "Kelvins"]}, Quantity[1, "Microflicks"]]]; Then: UnitConvert[ Lfull[Quantity[6, "Micrometers"], Quantity[330, "Kelvins"], .2, .3], "Microflicks"] gave: Quantity[194.166, "Microflicks"] and LfullNU[6, 330, .2, .3] gave: 55601.9 $\endgroup$ – c186282 Nov 3 '20 at 23:45
  • $\begingroup$ @c186282 I've diagnosed what the problem is. I imagine it's the same sort of problem that you've been running into, but in this case it has to do with the Piecewise-defined functions. Essentially, I'm taking advantage of the fact that Mathematica simplifies products of Quantity's automatically if the one unit is simply a derived unit (as microflicks is derived from products of the base SI units). However, here, the derived units are inside the Piecewise-defined function, and Mathematica doesn't automatically push the 1/Quantity[1,"microflicks"] inside. I'll work on it. $\endgroup$ – march Nov 4 '20 at 0:49
  • $\begingroup$ @c186282. So, I can imagine that this could get klugier and klugier. The current version now works on your function, but I'm trying to think about places that it can break. I treated the Piecewise-defined function case separately, but I can imagine that there are other cases that will need to be expanded before the quantities can work out right. It's likely not to be very robust. $\endgroup$ – march Nov 4 '20 at 1:20
  • $\begingroup$ Talking about a kluged solution, have you seen mine? Thanks for the update! Your solution worked for me. I do want to test it some more because once I take a solution it will be placed deep in the heart of my code and I need to know I can trust it. The key thing I want to test is if I have two interpolation functions one uses nm and the other uses micrometers and both are functions of the same variable. This is common and I should have added it in my test set. Thanks again. $\endgroup$ – c186282 Nov 4 '20 at 13:00
  • $\begingroup$ Sorry, it took so long I just wanted to test it on a few more functions. Thanks, I love the simplicity of this solution it is spot on! $\endgroup$ – c186282 Nov 6 '20 at 12:00

Based on Mathematica documentation, and alternative solution could be implement a simple replacement rule as follows in the next function:

NoUnits[f_] := f /. q_Quantity :> QuantityMagnitude[q]

Then, implement the Plot as a pure function to increase drastically the execution time:

    {l, .1, 15}, Frame -> True, 
    FrameLabel -> {"\[Lambda] (\[Mu]m)", 
      "\[ScriptCapitalL] (\[Mu]flicks)"}, ImageSize -> Medium]] &@ NoUnits[BB[Quantity[l, "Micrometers"], Quantity[333, "Kelvins"]]]

The result will be:

Result of the code proposed

  • $\begingroup$ This doesn't answer the real question. The plotting was just an example of how slow units can be. Please see @march response which is much more elegant than my attempt at the solution. $\endgroup$ – c186282 Nov 5 '20 at 11:25
  • $\begingroup$ Am sorry, but indeed it solves the real question: try running the following code using my function NoUnits: NoUnits[Lfull[Quantity[l, "Micrometers"], Quantity[T, "Kelvins"], a, b]] and you will get a version of the equation with no units. $\endgroup$ – Juanito970 Nov 5 '20 at 19:35
  • $\begingroup$ You can use a more elaborate example using unit convert to test the function proposed: NoUnits[Lfull[UnitConvert[Quantity[l, "Micrometers"], "Meters"], Quantity[T, "Kelvins"], Quantity[a, "Kilograms"], Quantity[b, "Parsecs"]]] , the result will be an equation with no units, as you requested in the title of your question. $\endgroup$ – Juanito970 Nov 5 '20 at 19:39
  • $\begingroup$ Yes, it removes the units but I also want the correct numerical values out. I did: LfullNU[l_, T_, a_, b_] := Evaluate[NoUnits[ Lfull[Quantity[l, "Micrometers"], Quantity[T, "Kelvins"], a, b]/ Quantity[1, "Microflicks"]]] Then: UnitConvert[ Lfull[Quantity[6, "Micrometers"], Quantity[330, "Kelvins"], .2, .3], "Microflicks"] gave Quantity[194.166, "Microflicks"] and LfullNU[6, 330, .2, .3] gave 55601.9 $\endgroup$ – c186282 Nov 5 '20 at 22:13
  • $\begingroup$ I would suggest that you change the order of Evaluate because it alter the order of the replacement rule , for example, the following code does give the correct result: NoUnits[Evaluate[ UnitConvert[ Lfull[Quantity[l, "Micrometers"], Quantity[T, "Kelvins"], a, b], "Microflicks"]]], this can be tested as you stated LfullNU[6, 330, .2, .3] and gives 194.166; The correct numerical value $\endgroup$ – Juanito970 Nov 6 '20 at 5:54

The issue is that your definition of NoUnits uses Map which only maps at level one.

An alternative solution is:

  outunits_] := ((fun /. 
      HoldPattern[Quantity[x__]] :> 
       UnitConvert[Quantity[x], outunits]) /. {UnitConvert[a_, 
       outunits] :> a}) /. {Quantity[q_, u_] :> q}
  • $\begingroup$ This fixes the issue at the top level of the unit conversion but now it fails in the exponents. Please try your function on the linear combination of the blackbodies. $\endgroup$ – c186282 Oct 30 '20 at 20:05
  • $\begingroup$ @c186282 I am not sure what you mean. If I call my definition myNoUnits and run your code then I get a BBnu[l, T1] + b BBnu[l, T2] === myNoUnits[(a BB[l, T1] + b BB[l, T2]) /. {l -> Quantity[l, "Micrometers"], T1 -> Quantity[T1, "Kelvins"], T2 -> Quantity[T2, "Kelvins"]}, "Microflicks"] evaluating to True. Perhaps you should clarify in your question what the expected output should be for each case. $\endgroup$ – Natas Nov 2 '20 at 11:59
  • $\begingroup$ I copied your code above and the right-hand side of your expression in the comment and I get an error: Quantity::compat: PureUnities and Microflicks are incompatible units. Also in the exponent, I get Exp[$Failed] But you are right I should clarify my goal in the question. $\endgroup$ – c186282 Nov 2 '20 at 21:40
  • $\begingroup$ Interesting, what is your $Version? $\endgroup$ – Natas Nov 3 '20 at 6:19
  • $\begingroup$ I'm using version 11.3 on Linux. I have the latest version 12 available but I need to support 11.3 for a while longer. :( However, when I get a chance I'll test in 12. What version are you using? $\endgroup$ – c186282 Nov 3 '20 at 12:34

If already in the input of the actual notebook: KnownUnitQ["Microflicks"] gives True. So there is need for someaction to get this referentially connceted to other known unit entities.

Mathematica V 12.0.0 has FormularLookup and FormulaData for lot of the questions purposes.

FormulaLookup["planck's law"]

enter image description here

The unit system in Physics is in need to be closed both operatinal and pur informal. As shown on the documentation page for Quantity.

enter image description here (* True *)

Normal[Quantity[10, "Percent"]] (* 1/10 *)

UnitDimensions["PartsPerMillion"] (* {} *)

For example

N[Quantity[1, "GravitationalConstant"], 20] (* 1 G *)

UnitConvert[N[Quantity[1, "GravitationalConstant"], 20]] (* Quantity[6.674*10^-11, ("Meters")^3/("Kilograms" ("Seconds")^2)] *)

This input does the required job:

equation = FormulaData[{"PlanckRadiationLaw", "Wavelength"},
    {"T" -> Quantity[5000, "Kelvins"], 
     "\[Lambda]" -> Quantity[wl, "Micrometers"]}
    ][[2, 2]];

Plot[equation, {wl, 0.1, 5}, 
 AxesLabel -> {"Wavelength [\[Mu]m]", 
   "Spectral radiance [W \!\(\*SuperscriptBox[\(sr\), \
\(-1\)]\)\!\(\*SuperscriptBox[\(m\), \(-3\)]\)]"}]

enter image description here

Like with the Unit it is not necessary that the forumla is a Physics one. If the running kernel knows it suffices.

To works with Units or Quantity is already optimized as a built-in. It is generalzied to match the requirement of a physical CAS completely. You have QuantityUnit and QuantityMagnitude for operatinal optimization. QuantityMagnitude works always for separation of magnitude from unit at most for visualization purposes.

You have only this built-in at hand. The solutions of the others answerer make use of that by the use of Quantity.

On the documentation page of QuantityMagnitude is an section about how Mathematica deals with unknown units. It is interpreted semantically.

QuantityMagnitude[Quantity[100, "Yards"], "Meterss"]

(* 2286/25 *)

There is the built-in IndependentUnit represents a unit string with no relationship to other units within a Quantity.

KnownUnitQ returns True for valid IndependentUnit specifications:


(* True *)

Names of untis have to be unique. And can be checked with CompatibleUnitQ. UnitDimensions is the logical center of the system:

UnitDimensions (* {{"ElectricCurrentUnit", -1}, {"LengthUnit", 2}, {"MassUnit", 1}, {"TimeUnit", -3}} *)

If this properly resolves Your Flicks You can do everything even get the input QuantityUnits out of the equation. And You can check this with it too.




{{"AngleUnit", -2}, {"LengthUnit", -1}, {"MassUnit", 
  1}, {"TimeUnit", -3}}

This proves my first solution in depths.

But take are:


enter image description here

(* UnitDimensions["1"] *)

So Your target is


(* {} *)

**Numerical values are considered dimensionless!**

Simplifyication works this way:


Interrogation of units works this way:

Composition[QuantityUnit, UnitConvert, Quantity] /@ {"PlanckConstant",
   "BoltzmannConstant", "MolarGasConstant", "AvogadroNumber"}

{("Kilograms" ("Meters")^2)/("Seconds"), ("Kilograms" ("Meters")^2)/(
 "Kelvins" ("Seconds")^2), ("Kilograms" ("Meters")^2)/(
 "Kelvins" "Moles" ("Seconds")^2), "DimensionlessUnit"}

The selection of fundamental constants can be replace by Your interests.

Good practice to work with units in graphical representation is

s = Quantity[100, "Meters"] - Quantity[9.8, "Meters/Seconds^2"]*t^2;
factor = QuantityMagnitude[Quantity[1, "Meters"], "Feet"];
Plot[s*factor, {t, Quantity[0, "Seconds"], Quantity[3, "Seconds"]}, 
 FrameLabel -> Automatic, PlotTheme -> "Web"]

enter image description here

from this answers: how do i properly use quantities units in plots?.

To get more free space for units make use of getting useful units for combinations of physical constants like on wolframalpha.

This should do the job: specify set of base units to use in unitconvert as I understand the question.


(*a set of standard units that are used when not specified*)
siUnits = {"Seconds", "Meters", "Kilograms", "Kelvins", 
   "KelvinsDifference", "Amperes", "Candelas", "Moles", "Radians"};
siUnitDimensions = UnitDimensions[#][[1, 1]] & /@ siUnits;

makeUnitSystem::dependent = 
  "The unit system `1` is overcomplete. Please remove some unit.";
makeUnitSystem[] = Thread[siUnitDimensions -> siUnits];
makeUnitSystem[L_List] := 
 Module[{M, n, u},(*convert the desired unit system to base units*)
  M = Lookup[#, siUnitDimensions, 0] & /@ 
    Apply[Rule, UnitDimensions /@ L, {2}];
  If[MatrixRank[M] < Length[L], 
   Message[makeUnitSystem::dependent, L];
  (*check which base units cannot be expressed in this system*)
  n = Position[Diagonal[PseudoInverse[M].M], Except[1], {1}, 
    Heads -> False];
  (*extend the unit system if necessary*)
  If[Length[n] > 0, 
   Return[makeUnitSystem[Append[L, siUnits[[n[[1, 1]]]]]]]];
  (*find the compound units that represent the base units*)
  u = Times @@@ Transpose[L^Transpose[PseudoInverse[M]]];
  (*return replacement list*)Thread[siUnitDimensions -> u]]

  unitSystem_ /; VectorQ[unitSystem, Head[#] === Rule &]] := 
 UnitConvert[x, Times @@ Power @@@ (UnitDimensions[x] /. unitSystem)]

With this You may check whether the dimensionless parameters suit into the systemes international de unites for a complete system. This is inspired by the answer of Roman Maeder, @roman.

This makes a consistent dimensionless system of dimensions provable. I does not make the choice. It does not make any calculations faster. This is done internally and at fastest if consistency is True.

This is for converting units, an example:

SetAttributes[withUnits, HoldAll];
withUnits[code_] :=
       SetAttributes[Quantity, HoldRest];
       Quantity /: UnitConvert[arg_, Quantity[_, unit_]] :=
          UnitConvert[arg, unit];
       Quantity /: Times[0, Quantity[_, unit_]] :=
          Quantity[0, unit];
          m = Quantity[1, "Meters"], 
          s = Quantity[1, "Seconds"],       
          min =  Quantity[1, "Minutes"],
          km = Quantity[1, "Kilometers"]

withUnits[UnitConvert[1 m/s^2*(1 min)^2,km]]

from simpler input for the new unit support by @leonid-shifrin.

k[uqty_] := Quantity[uqty, "Kelvins"]
um[uqty_] := Quantity[uqty, "Micrometers"]

FormulaData[{"PlanckRadiationLaw", "Wavelength"}, 
{"T" -> k@5000, "\[Lambda]" -> um@\[Lambda]}]

gives a formular in dimensionless [Lambda]:

enter image description here

  • $\begingroup$ Thanks but not really what I'm looking for. As far as Microflicks, there was a time when Mathematica had it wrong: mathematica.stackexchange.com/questions/111818/… $\endgroup$ – c186282 Nov 4 '20 at 13:13
  • $\begingroup$ Why do You reference to this? The first comment already states this in not a bug! $\endgroup$ – Steffen Jaeschke Nov 5 '20 at 9:55

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.