# Numerical approximation of symbolic equation

What I would like to do is the following:

sage: var('theta')
theta
sage: theta=25*units.angles.degree
sage: print theta
25*degree
sage: a=sin(theta)
sage: print a
sin(25*degree)


At this point, I type 'a.' and tab and in the popup list of functions I find numerical_approx. So I type in:

sage: a.numerical_approx()
Traceback (click to the left of this block for traceback)
...
TypeError: cannot evaluate symbolic expression numerically


Okay, it would appear that adding the units to theta turns it into a symbolic variable which turns my expression a into a symbolic expression. But what I would like to do is enter my angle theta in either degrees or radians and get a numerical approximation of the answer based on the units of the input. Is there a clean way to do this?

edit retag close merge delete

Sort by » oldest newest most voted

In the show-your-work spirit, here was my line of thought. First, I thought that maybe all that was missing was adding a .n method to degree, to get degree to pretend as if it were pi/180, but that didn't work like I was hoping. Then I realized that this was silly: pi already behaves the way I want it to, I just want to define a new pi-like constant with name "degree" and value "pi/180". So looking at how pi is implemented, it subclasses Constant and then takes the expression from that. So in a similar spirit:


from sage.symbolic.constants import Constant

class ConstantFromExpression(Constant):
def __init__(self, name, expr):
conversions = dict(maxima=repr(maxima.coerce(expr)))
Constant.__init__(self, name,conversions=conversions)
self._expr = expr
def __float__(self): return float(self._expr)
def __complex__(self): return complex(self._expr)
def _mpfr_(self, R): return R(self._expr)
def _real_double_(self, R): return R(self._expr)

def NamedExpression(name, expr):
return ConstantFromExpression(name,expr).expression()


which should allow you to do what you want:


sage: degree = NamedExpression('degree', pi/180)
sage: degree
degree
sage: float(degree)
0.017453292519943292
sage:
sage: sin(25*degree)
sin(25*degree)
sage: a = sin(25*degree)
sage: a.n()
0.422618261740699
sage: RealField(100)(a)
0.42261826174069943618697848965
sage:
sage: sin(90*degree)
sin(90*degree)
sage: simplify(sin(90*degree))
1
sage: sin(45*degree)
sin(45*degree)
sage: simplify(sin(45*degree))
1/2*sqrt(2)


There are lots of conversions and features I didn't implement, of course, but I had to do at least the maxima one to get simplify to work correctly. I don't know enough about the maxima coercion repr to trust it entirely -- should probably put some tests in, to make sure the result tests equal to the input. But for this trivial case, it seems to work okay.

Maybe a NamedExpression wrapper would be useful?

more

Thank you for all the hard work. This definitely looks promising. I will experiment with it to see how it goes. It seems to me that if this sort of thing was implemented generally it would make sage a *much* more attractive package.

( 2010-12-31 12:03:50 -0500 )edit

Let me know if there's something you need it to do that it doesn't. And if it works, come back here and accept my answer. :^)

( 2011-01-01 23:33:27 -0500 )edit

If this works out, presumably all the dimensionless angle units in sage.symbolic.units should be treated this way - hopefully it wouldn't be too hard to do so. DSM, if so, this would be worth raising on sage-devel, since there was a lot of discussion of the best way to implement units there when it happened.

( 2011-01-03 14:21:56 -0500 )edit

@kcrisman: ISTM modifying the above to make it a real unit so that sin(5*degree) would work but sin(5*degree**2) didn't would be a bit tough though.

( 2011-01-04 15:21:34 -0500 )edit

Good point. What do unit packages in general programming languages do in these cases?

( 2011-01-04 17:09:55 -0500 )edit

Is there any way to split up a symbolic expression like '25*degree' into units and value? some function like:

sage: t = 25*units.angles.degree;
sage: get_numeric(t)
25.0
sage: t/get_numeric(t)
degree


This will make possible the usage of all functions, designed for numeric arguments, with units of measurement,after some simple redefinition. As the units work like symbolic variables, we can set them to unity to get a numerical value:

sage: float(t(units.angles.degree=1))
25.0


For that, one has to know the units involved. it can be achieved by arguments() function:

sage: t.arguments()
(degree,)


Putting all together:

def get_numeric(t)
dic={};
for arg in t.arguments():
dic[arg]=1;
return (float(t(dic)));


Don't know about the performance of such method.

more