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
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
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.
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
Asked: 2024-06-04 12:52:23 +0100
Seen: 209 times
Last updated: Jun 05
Copyright Sage, 2010. Some rights reserved under creative commons license. Content on this site is licensed under a Creative Commons Attribution Share Alike 3.0 license.
Like this
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:
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.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).