1

Sadly I'm having trouble with ezunits and plotting again. This time in combination with if-clause. This works, but probably not as it should:

/*just works*/
M_K1q(n):= if n <= 500 then 0
    else (n-500)^2/20000;
wxdraw2d(explicit(qty(M_K1q(n)), n, 0, 1500));

/*kinda works*/
M_L(n):= (10+n/100`minute)`N*m; /*doesn't*/
M_L(n):= (10+(n`minute)/100)`N*m; /*doesn't*/
M_L(n):= (10`N*m)+n/100`N*m*minute; /*works*/
wxdraw2d(explicit(qty(M_L(n)), n, 0, 1500));

/*but both work that way*/
M_L(n):= (10+n/100`minute)`N*m;
M_L(n):= (10+(n`minute)/100)`N*m;
wxdraw2d(explicit(qty(qty(M_L(n))), n, 0, 1500));

But with units I cannot plot it, have tried many variations of it:

M_K1(n):= block(
    (if qty(n) <= 500 then 0
    else ((n`minute)-500)^2/20000)`N*m
);
M_K1(99`1/minute); /*works*/
wxdraw2d(explicit(qty(M_K1(n)), n, 450, 600)); /*doesn't*/

Input n should be in 1/minute and output should be N*m. Looks like qty() isn't working properly, since the units weren't removed completely:

draw2d (explicit): non defined variable in term: 5.0*10^-5*((realpart(501.7241379310345 ` minute)-500.0)^2-1.0*imagpart(501.7241379310345 ` minute)^2)
 -- an error. To debug this try: debugmode(true);

2022-05-29 EDIT:

Yes, using qty() would require entering n in 1/minute, which would be an acceptable workaround, if graphs and solving would work.

Damn it, I tried so many versions. That's the one with no unit at the 0 because it's "global":

M_K1(n):= (if n <= 500`1/minute then 0
    else (n`minute-500)^2/20000)`N*m;

Shouldn't the last one also work?

(%i102) M_K1(1000`1/minute);
(%o102) 25/2 ` N*m
(%i108) M_K1(2.5`1/s);
(%o108) 0 ` N*m
(%i111) M_K1((1000/60)`1/s)``N*m;
(%o111) (50/3 ` minute/s-500)^2/20000 ` N*m

Drawing graphs with that takes a minute or two and then crashes:

(%i24)  wxdraw2d(explicit(qty(M_K1(n)), n, 450, 600));
draw2d (explicit): non defined variable in term: 
realpart(if 455.1724137931034<=500.0 ` 1/minute then 0.0 else 5.0*10^-5*(455.1724137931034 ` minute-500.0)^2)
 -- an error. To debug this try: debugmode(true);

Shouldn't it work though?

It looks correct, but using your:

M_K2(n):=if n <= 500 ` 1/minute then 0 ` N*m
            else (n-500 ` 1/minute)^2/(45000 ` 1/(kg*m^2));

I get:

(%i40)  M_K2(1000`1/minute)``N*m;
(%o40)  1/648 ` N*m

That's not right!

(%i78)  M_K2(n):= if n <= 500`1/minute then 0`N*m
        else (n`minute-500)^2/45000`N*m;
(%o78)  M_K2(n):=if n<=500 ` 1/minute then 0 ` N*m else (n ` minute-500)^2/45000 ` N*m
(%i79)  M_K2(1000`1/minute);
(%o79)  50/9 ` N*m

That's more like it.

Also, it should be solvable:

(%i69)  res: dimensionally(solve(M_K2(n) = 10`N*m, n)); n_K2P2: rhs(res[1])``1/minute; float(%);
(%o67)  [(if n<=500 ` 1/minute then 0 ` N*m else (n^2 ` N*m*minute^2+(-1000*n) ` N*m*minute+250000 ` N*m)/45000)=10 ` N*m]
(%o68)  (600*kg*m^2)/s ` 1/minute
(%o69)  (600.0*kg*m^2)/s ` 1/minute

Should look like this:

(%i86)  res: dimensionally(solve(M_K2(n) = 10`N*m, n)); n_K2P2: rhs(res[2])``1/minute; float(%);
(%o84)  [n=(500-12*5^(5/2)) ` 1/minute,n=(12*5^(5/2)+500) ` 1/minute]
(%o85)  (12*5^(5/2)+500) ` 1/minute
(%o86)  1170.820393249937 ` 1/minute

If I get one thing working, another one breaks. Looks like one edge-case after the next.

Dux
  • 101
  • 4

2 Answers2

1

I think this formulation works okay.

(%i4) M_K1(n):=if qty(n) <= 500 then 0 ` N*m
            else (qty(n)-500)^2/20000 ` N*m;
(%i5) draw2d(explicit(qty(M_K1(n ` 1/minute)),n,450,600));

In M_K1, I'm being careful to return a result which has units on both branches of the if, and supplying an argument to M_K1 which has units. Note that M_K1 returns

(%i6) M_K1(400 ` 1/minute);
(%o6)                        0 ` N m

on the less than or equal to 500 branch, and

(%i7) M_K1(600 ` 1/minute);
                             1
(%o7)                        - ` N m
                             2

on the greater than 500 branch.

This formulation of M_K1 requires that you know the units of its argument n, since the function says qty(n) and then works with that number. So you couldn't change the units to 1/hour or something without also changing the constants in the function.

To make M_K1 behave consistently with different, equivalent units, one can make a couple of changes. The first is to write if n < 500 ` 1/minute then .... Ezunits does implement logical comparisons with units, but it's slow. E.g.

(%i8) if 375 ` 1/minute <= 500 ` 1/minute then foo else bar;
(%o8)                          foo

The other is to put units on the constants 500 and 20000 in the arithmetic expression. Then M_K1 returns a correct result even when the units for its argument are changed.

(%i9) M_K1(n):=if n <= 500 ` 1/minute then 0 ` N*m
            else (n-500 ` 1/minute)^2/(20000 ` 1/(kg*m^2));
                                 1
(%o9) M_K1(n) := if n <= 500 ` ------ then 0 ` N m
                               minute
                                                           1    2
                                              (n - 500 ` ------)
                                                         minute
                                         else -------------------
                                                           1
                                                 20000 ` -----
                                                             2
                                                         kg m

Quantity > 500 but units are smaller, so it's less than 500 Hz:

(%i10) M_K1(1000 ` 1/hour);
(%o10)                       0 ` N m

Quantity < 500 but units are bigger, so it's greater than 500 Hz:

(%i11) M_K1(1 ` kHz);
             1         2                         1    2
(%o11)    (----- ` kg m ) (1 ` kHz + (- 500) ` ------)
           20000                               minute

Oh, that's messy, does it have the dimensions we expect, namely energy?

(%i12) foo:expand(%)
                 2                                          2
       25    kg m       1        2     2      1     kHz kg m
(%o12) -- ` ------- + ----- ` kHz  kg m  + (- --) ` ---------
       2          2   20000                   20     minute
            minute
(%i13) dimensions(foo)
                                 2
                         3 length  mass
(%o13)                   --------------
                                 2
                             time
(%i14) dimensions(J)
                                2
                          length  mass
(%o14)                    ------------
                                 2
                             time

Good news, %o13 has same dimensions as J. (Minor wart: constant factor 3 appearing in %o13.) Okay, now convert to J:

(%i15) foo `` J
                            14161
(%o15)                      ----- ` J
                             288

Just for fun, nonstandard units of energy. lbf is "pound force", distinguished from lbm, "pound mass".

(%i16) foo `` foot*lbf
                  22578125000000000
(%o16)            ----------------- ` foot lbf
                   622569466070541
Robert Dodier
  • 16,905
  • 2
  • 31
  • 48
0

Thanks for your continued interest. I'll try to address the problems one by one. Overall the main thing I see is that units aren't used consistently.

(1)

M_K1(n):= (if n <= 500`1/minute then 0
    else (n`minute-500)^2/20000)`N*m;

n must have dimensions of 1/time -- okay. The "then" branch returns 0 ` N*m -- okay. On the "else" branch, I guess you want to convert n to 1/minute to match 500.

n`minute almost works. To be clear about what you want, you can say n`minute `` 1 or n/(1 ` 1/minute) `` 1 or n*(1 ` minute) `` 1 or qty(n `` 1/minute). (The foo `` 1 reduces a units expression which is actually nondimensional to a non-units expression.)

I think qty(n `` 1/minute) expresses the intent most clearly, but here is a formulation which says, "divide a unit by a compatible unit" -- I know this is conventional in applications.

(%i3) M_K1(n):= (if n <= 500`1/minute then 0
    else (n/(1 ` 1/minute)``1 - 500)^2/20000)`N*m;
                                  1
(%o3) M_K1(n) := (if n <= 500 ` ------ then 0
                                minute
                                       n                 2
                                  (---------- `` 1 - 500)
                                         1
                                   1 ` ------
                                       minute
                             else ------------------------) ` N m
                                           20000

Now I get:

(%i4) M_K1(1000`1/minute);
                            25
(%o4)                       -- ` N m
                            2
(%i5) M_K1(2.5`1/s);
(%o5)                        0 ` N m
(%i6) M_K1((1000/60)`1/s);
                            25
(%o6)                       -- ` N m
                            2

(2) About the plot, it's necessary to ensure the argument is dimensional.

draw2d(explicit(qty(M_K1(n ` 1/minute)), n, 450, 600));

seems to have the desired effect.

The dimensional comparison if n <= 500 ` 1/minute ... is correct, but it's slow. An equivalent formulation which is much faster is

M_K1a(n) := (if qty(n `` 1/minute) <= 500 then 0
    else (n/(1 ` 1/minute)``1 - 500)^2/20000)`N*m;
draw2d(explicit(qty(M_K1a(n ` 1/minute)), n, 450, 600));

(3) About M_K2, the result %o40 looks correct to me. I get:

(%i21) M_K2(1000 ` 1/minute);
                                    2
                          50    kg m
(%o21)                    -- ` -------
                          9          2
                               minute

That seems right because

(%i22) 500^2/45000;
                               50
(%o22)                         --
                               9

and then substituting seconds for minutes,

(%i23) subst(minute = 60*s, %o21);
                                     2
                            1    kg m
(%o23)                     --- ` -----
                           648     2
                                  s
(%i24) % `` J;
                              1
(%o24)                       --- ` J
                             648

The second formulation for M_K2 which you showed, %i78, is confusing seconds and minutes -- the numerator must be 1/minute^2 but the way it's written, the units are lost. Then N is introduced, which contains 1/time squared which must be measured in 1/second^2.

(4) About solve, the main problem is that it doesn't know what to do with if. We can look for solutions greater than 500 ` 1/minute via assume. (Ideally we wouldn't have to help out like that; it's a work-around.)

In a case like this, I'm inclined to say assume(foo > 500) and then work with foo ` 1/minute, because I'm doubtful about assume handling units correctly. However, I find that the more obvious assume(n > 500 ` 1/minute) has the expected effect.

(%i49) assume (n > 500 ` 1/minute);
                                    1
(%o49)             [n + (- 500) ` ------ > 0]
                                  minute

That statement isn't really correct, since 0 is not the same thing as 0 ` 1/minute. Let's forge ahead regardless. (We could be more careful by saying assume(%n > 500) and working with %n ` 1/minute here.)

(%i50) dimensionally (solve (M_K2(n) = 10`N*m, n));
                          1
(%o50) [n = ((- 1) ` -----------)
                     kg m minute
      5/2
 (12 5    ` sqrt(N) sqrt(kg) sqrt(m) minute + (- 500) ` kg m), 
              1            5/2
n = (1 ` -----------) (12 5    ` sqrt(N) sqrt(kg) sqrt(m) minute
         kg m minute
 + 500 ` kg m)]

Hmm, those units are very messy, let's make it more readable:

(%i51) map (lambda ([e], e `` 1/minute), %);
               1                  7/2      1
(%o51) [n `` ------ = (500 - 144 5   ) ` ------, 
             minute                      minute
                                1            7/2            1
                         n `` ------ = (144 5    + 500) ` ------]
                              minute                      minute
(%i52) float (%);
               1                                1
(%o52) [n `` ------ = (- 39749.22359499623) ` ------, 
             minute                           minute
                               1                            1
                        n `` ------ = 40749.22359499623 ` ------]
                             minute                       minute

Substituting a supposed solution back into M_K2:

(%i53) M_K2(40749.22359499623 ` 1/minute);
                                            2
                                        kg m
(%o53)             36000.00000000002 ` -------
                                             2
                                       minute
(%i54) % `` N*m;
(%o54)               10.00000000000001 ` N m
Robert Dodier
  • 16,905
  • 2
  • 31
  • 48