Ask Your Question
2

I cannot take sqrt of units

asked 2024-06-04 12:52:23 +0200

DylanDylanDylan gravatar image

updated 2024-06-07 09:42:13 +0200

FrédéricC gravatar image

I'm trying to use sagemath to evaluate equations that involve units. And when I hit a sqrt of a unit it won't evaluate it.

For example:

sage: sqrt(4*units.length.meter^2)
2*sqrt(meter^2)

The result should be 2*meter but I can't figure out how to force it to evaluate the sqrt of m^2. Even if I try to convert it like:

(sqrt(4*units.length.meter^2)).convert(units.length.meter)

I get:

ValueError: Incompatible units
edit retag flag offensive close merge delete

Comments

2

Like this

sage: sqrt(4*units.length.meter^2)
2*sqrt(meter^2)
sage: _.canonicalize_radical()
2*meter
FrédéricC gravatar imageFrédéricC ( 2024-06-04 20:11:09 +0200 )edit

Thanks. That does actually get me 1 step further. But this was a stripped down version of my overall equation. I now have another pesky part of the equation to deal with:

 sage: stuck=123.0*sqrt(units.length.meter*(units.length.meter/units.time.second^2)/(17.0*tan(0.1*pi) + 1))/cos(0.2*pi)
sage: stuck.canonicalize_radical()
123.0*meter/(second*sqrt(17.0*tan(0.1*pi) + 1)*cos(0.2*pi))
sage: stuck.canonicalize_radical().polynomial(RR)
TypeError: unable to coerce since the denominator is not 1

In some simpler equations polynomial(RR) helps. Is there something I'm doing wrong here constructing and evaluating these equations?

I can do a similar thing in Mathematica and it allows me to call N[] and it reliably gives me back 1 number and a unit.

DylanDylanDylan gravatar imageDylanDylanDylan ( 2024-06-05 05:42:54 +0200 )edit

Another alternative (after canonicalizing the radical) is to use an ExpressionTreeWalker to substitute numerical approximations (that should probably be built-in functionality, hidden inside a new method).

rburing gravatar imagerburing ( 2024-06-05 14:05:35 +0200 )edit

1 Answer

Sort by » oldest newest most voted
2

answered 2024-06-05 12:24:53 +0200

rburing gravatar image

updated 2024-06-05 13:54:10 +0200

Here is something that works on your second example:

def evaluate_with_units(expr):
    expr_canonical = expr.canonicalize_radical()
    var_names = [str(v) for v in expr_canonical.variables()]
    poly_ring = PolynomialRing(RR, names=var_names)
    num = expr_canonical.numerator().polynomial(ring=poly_ring)
    denom = expr_canonical.denominator().polynomial(ring=poly_ring)
    return SR((num.lc() / denom.lc()) * (num/num.lc()) / (denom / denom.lc()))

Example:

sage: stuck=123.0*sqrt(units.length.meter*(units.length.meter/units.time.second^2)/(17.0*tan(0.1*pi) + 1))/cos(0.2*pi)
sage: evaluate_with_units(stuck)
59.5254442433757*meter/second
edit flag offensive delete link more

Comments

Nice. This does work for my example! I don't mind defining a function for these if it works reliably. But I notice this requires me to define all the units I'm using with RR['meter,second']. Is there any way to avoid that? I tried without and I got an error KeyError: 'variable_names_recursive'.

I am hoping to use sagemath as a general purpose physics calculator with units like I have been using Mathematica but it seems like there are some fundamental pieces missing here.

DylanDylanDylan gravatar imageDylanDylanDylan ( 2024-06-05 12:34:55 +0200 )edit

I updated the function to determine the variable names automatically.

rburing gravatar imagerburing ( 2024-06-05 13:55:41 +0200 )edit
2

Thanks! That works for me even with my complex examples. One caveat I found is that when it contains units.acceleration.gravity it's not smart figuring out to convert that in to m/s^2 and multiply things out.

But I did find that if I just do convert() the whole thing beforehand it gives me back what I want. So finally I ended up with:

def evaluate_with_units(expr):
    expr = expr.convert()
    expr_canonical = expr.canonicalize_radical()
    var_names = [str(v) for v in expr_canonical.variables()]
    poly_ring = PolynomialRing(RR, names=var_names)
    num = expr_canonical.numerator().polynomial(ring=poly_ring)
    denom = expr_canonical.denominator().polynomial(ring=poly_ring)
    return SR((num.lc() / denom.lc()) * (num/num.lc()) / (denom / denom.lc
DylanDylanDylan gravatar imageDylanDylanDylan ( 2024-06-06 08:36:32 +0200 )edit

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

1 follower

Stats

Asked: 2024-06-04 12:52:23 +0200

Seen: 138 times

Last updated: Jun 05