Unit conversion

edit

It seems that the units end up in sqrt() expressions. Even if you have something like sqrt(second^2), Sage does not simplify this automatically. This is a bug caused by the fact that second is a symbolic variable and unless you explicitly state the domain for a variable it is assumed to be complex. For a complex variable x, we can't say that sqrt(x^2) = x. But if x > 0, we can:

sage: var('y', domain='positive')
y
sage: sqrt(y^2)
y
sage: sqrt(x^2)
sqrt(x^2)

There is also the fact that sqrt() only performs trivial simplifications automatically. This is not simplified for instance:

sage: sqrt(y^2/(x+1))
sqrt(y^2/(x + 1))

Because of another bug, even though we already declared y > 0, we need to state this to Maxima before we can simplify further:

sage: t = sqrt(y^2/(x+1))
sage: t.simplify_radical()
abs(y)/sqrt(x + 1)
sage: assume(y > 0)
sage: t.simplify_radical()
y/sqrt(x + 1)

You might be able to work around your problems by combining some of the tricks above. It would be great if someone files a ticket to declare units as positive.

For now, I find it more useful to define my own units instead of using the build-in Sage units. For example, I would rewrite this code as follows:

# Setup units
#    Define base units (kilogram, meter, second) and secondary units (lbf, inch, sec)
kilogram  = var('kilogram', domain='positive'); assume(kilogram,'real');
meter = var('meter', domain='positive'); assume(meter,'real');
second  = var('second', domain='positive'); assume(second,'real');
lbf  = 4.44822162*kilogram*meter/second^2
inch = 0.0254*meter
sec  = 1*second;
psi  = lbf / inch^2;
gal  = 231.0*inch^3;
min  = 60.0*sec;
gpm  = gal/min;
hole = gpm/sqrt(psi);

# Define the given constraints
FLext = 20000.0*lbf;
vext  = 15.0*inch/sec;
FLret = -10000.0*lbf;
vret  = -15.0*inch/sec;

# Assumed values
Ps = 3000.0*psi;
pv = 1.0;
pc = 1.4;

# System of 3 equations, 3 unknowns
ABE = var('A_BE');
ARE = var('A_RE');
kv = var('k_v');
eq1 = vext^2 == (Ps*ABE - FLext) * kv^2 / ABE^3 / (1 + pv^2 * pc^(-3));
eq2 = vret^2 == (FLret + Ps*ARE) * kv^2 / ARE^3 / (1 + pv^2 * pc^3);
eq3 = pc == ABE / ARE;
soln = solve([eq1, eq2, eq3], [ABE, ARE, kv], solution_dict=True);
sol = soln[1];
show(sol[ABE])
show(sol[ARE])
show(sol[kv])

# Showing the units of various variables
#     Use .full_simplify()._convert(RR) to provide simplified results
show(sol[kv].full_simplify()._convert(RR))
show(hole.full_simplify()._convert(RR))
show((sol[kv]/hole).full_simplify()._convert(RR) ) # this should be dimensionless

# This final result is now dimensionless, as desired.

I am still playing around, but the main idea is define one's own units (as shown at the beginning of the above example) using "domain" and "assume" to ensure they are positive real. Methods such as ".full_simplify()" and "._convert(RR)" can then be used to produce simple numerical results with all the units correctly cancelled out.