NDSolve and Quantity[]

Question

I'm working on a set of relatively basic ODEs of a physical pendulum and motor:

$$\left\{\text{jb} \theta ''(t\text{s})+\text{rb} \theta '(t\text{s})+\alpha \sin (\theta (t\text{s}))+\mu \tan ^{-1}\left(\frac{\theta '(t\text{s})}{\iota }\right)=\text{jw} \phi ''(t\text{s}),\tau u(t\text{s})=\text{jw} \left(\theta ''(t\text{s})+\phi ''(t\text{s})\right)+\text{rw} \phi '(t\text{s})\right\}$$

In an attempt to understand the system some more, I've updated the code to use real world values via Quantity[]:

displacement = Quantity[5 \[Pi]/180, "Radians"];
params = {\[Mu] -> 0.000795, rw -> 0.000039255383907286545`, 
   rb -> 0.00016799999999999996`, \[Iota] -> 0.1, 
   g -> Quantity[9.81, ("Meters")/("Seconds")^2], 
   mb -> Quantity[88, "Grams"], mw -> Quantity[129.31, "Grams"], 
   x -> Quantity[60.519, "Millimeters"], 
   y -> Quantity[61.816, "Millimeters"], 
   l -> Quantity[86.50, "Millimeters"], 
   k -> Quantity[90, "Millimeters"], 
   jw -> Quantity[0.0004032497194086029`, "Kilograms" ("Meters")^2], 
   jb -> Quantity[0.001412, "Kilograms" ("Meters")^2], \[Tau] -> 
    Quantity[33.5, ("Millinewtons" "Meters")/("Amperes")]};

eqn1 = g (l mb + k mw) Sin[\[Theta][Quantity[t, "Seconds"]]] + 
   jb (\[Theta]^\[Prime]\[Prime])[Quantity[t, "Seconds"]];
eqn2 = jw (\[Phi]^\[Prime]\[Prime])[Quantity[t, "Seconds"]];
deqns = {eqn1 == 
    jw \[Phi]''[Quantity[t, "Seconds"]] - 
     rb \[Theta]'[
       Quantity[t, "Seconds"]] - \[Mu] ArcTan[\[Theta]'[
        Quantity[t, "Seconds"]]/\[Iota]] + fx + fy , 
   eqn2 == \[Tau] u[Quantity[t, "Seconds"]] - 
     jw \[Theta]''[Quantity[t, "Seconds"]] - 
     rw \[Phi]'[Quantity[t, "Seconds"]]} ;
ics = {\[Theta][Quantity[0, "Seconds"]] == 
    displacement, \[Theta]'[Quantity[0, "Seconds"]] == 
    Quantity[0, "Radians"] , \[Phi][Quantity[0, "Seconds"]] == 
    0, \[Phi]'[Quantity[0, "Seconds"]] == 0};

And from here try to to use NDSolve[]

sol = NDSolve[{{deqns /. params, ics} /. {u[t] -> 0, 
      f[t] -> 0}}, {\[Theta][Quantity[t, "Seconds"]], \[Phi][
     Quantity[t, "Seconds"]]}, {t, 0, 60}];
(*pen ={k Sin[\[Theta][t]],-k Cos[\[Theta][t]]} /. sol;*)

However, I get the following error:

NDSolve::litarg: To avoid possible ambiguity, the arguments of the dependent variable in [Theta][ts] should literally match the independent variables.

Is this error because NDSolve can't work with Quantity[], or how can I give units with NDSolve?

The resulting unit should be in Radians.

Thank you for the help!

Answer 1

I was also having the same trouble when working with NDSolve and units. However, I think I've found an elegant workaround: non-dimensionalization.

Now, mechanics is not my art, so I don't dare to help you directly with your example. Instead, I propose to work with a simpler ODE: the decay of a chemical compound following a 1st order kinetics, which is basically the same as radioactive decay. By using this simple example, we can focus more on the method, rather than in its solution. Moreover, besides solving the problem numerically, this ODE has an analytic solution so you can compare both.

To begin with, you can model the decay in Mathematica as

mod = {c'[t] == -k c[t], c[0] == c0}

where $c'[t],k,c0$ represent the speed at which the concentration of the compound decreases wrt time, the reaction coefficient and the initial concentration of the compound. Recall that the negative sign is there, since we are modeling consumption (decay).

Then, you can assign dimensions to the quantities using the QuantitiVarible function (note that I use capital letters for the quantity in the function, otherwise Mathematica confuses the assigned symbol with the quantity and its answers become messy):

c = QuantityVariable["C", "Mass"/"Volume"];
c0 = QuantityVariable["C0", "Mass"/"Volume"];
k = QuantityVariable["K", 1/"Time"];
t = QuantityVariable["T", "Time"];

From here, you can tell Mathematica to scale the model for you:

smod = NondimensionalizationTransform[mod, {c, t}, {cs, ts}, 
   "PropertyAssociation"];

You can then ask Mathematica for the scaled model and the scaled parameters, so you can go back to the dimensional model at any time:

{smod["ReducedForm"], smod["DimensionalizationRules"]}

which yields

{{Derivative[1][cs][ts] == -cs[ts], 
  cs[0] == 1}, {Derivative[1][cs][ts] -> 
   Derivative[1][QuantityVariable["C",(("Mass")/("Volume"))]][
    QuantityVariable["T","Time"]]/(
   QuantityVariable["C0",("Mass")/("Volume")] QuantityVariable[
    "K",1/("Time")]), 
  cs[ts] -> 
   QuantityVariable["C",(("Mass")/("Volume"))][QuantityVariable[
    "T","Time"]]/QuantityVariable["C0",("Mass")/("Volume")], 
  cs[0] -> QuantityVariable["C",(("Mass")/("Volume"))][0]/
   QuantityVariable["C0",("Mass")/("Volume")], 
  cs -> QuantityVariable["C",("Mass")/("Volume")]/QuantityVariable[
   "C0",("Mass")/("Volume")], 
  ts -> QuantityVariable["K",1/("Time")] QuantityVariable[
    "T","Time"]}}

Now the fun part: you can actually NDSolve the scaled model directly without needing to go back to the dimensional one:

ssl = NDSolve[smod["ReducedForm"], cs, {ts, 0., 1.}]

which gives an InterpolatingFunction you can plot:

Plot[cs[ts] /. ssl, {ts, 0., 1.}]

As you see, you have a cleaner way to deal with the NDSolve issue.

Finally, some troubles unrelated to NDSolve, but rather to NondimensionalizationTransform: 1) if in the ODE there is a non-zero, dimensionless parameter, then the scaling fails. Mathematica cannot scale equations with dimensionless parameters. It will tell you that the equation is "dimensionally unbalanced" or something else. I cannot find an answer to that issue yet. But that's another topic. 2) To the best of my knowledge, you have no control whatsoever on the way Mathematica scales equations.

I believe this will work for your problem, but let me give you an advice: scale the model with symbolic variables and then, just before solving the ODE numerically, replace the symbols with numbers.

Hope the answer helps.