Basic dimensional computation

Say, for example, you want to use DimPy to calculate the power that could be generated from a weir on the river Cam. Start by loading DimPy:

sage: from sage.dimpy import *

At the prospective site the river is 10 m wide, 1 m deep and flows at about 0.1 m/s. The change in height at the lock is about 1 m. Hence define the following variables:

sage: width = 10 * meter
sage: depth = 1 * meter
sage: speed = 0.1 * meter / second
sage: change_in_height = 1 * meter

This gives a rate of flow of:

sage: width * depth * speed
1.00000000000000 m^3 s^-1

To calculate the power it is necessary to know the strength of Earth’s gravity, which is pre-defined under the name g_n (there are lists of pre-defined units and constants in the modules quantity and additional_units):

sage: g_n
9.80665 m s^-2

Power is the rate of change of energy, in this case the potential energy of the water, which is the product of mass, change in height and g. The rate of flow calculated above is the volume rate of flow, but that can easily be changed by multiplying with the density.

sage: density = 1 * kilogram / litre
sage: power = width * depth * speed * density * change_in_height * g_n
sage: power
9806.65000000000 W

Power is indeed measured in watt; if the result had had a different unit that would have been an indication that something has gone wrong. But there are also different units for power, some of them more intuitive than the watt. The next section will explain how to display quantities in different units.

Unit conversions

There are different ways of telling DimPy which unit to display something in. Many widely used units are pre-defined and so you can type, for example:

sage: width.in_unit(foot)
32.8083989501 ft

But if you try a similar thing with a combination of units, say

sage: speed.in_unit(foot/second)
0.328083989501312 * (0.3048 m s^-1)

it doesn’t work as expected.

DimPy cannot figure out how to print a combination of units, but given a string containing a combination of pre-defined units, e.g. ‘ft/s’ it can find the corresponding value using the built in QuantityParser. While the main program only accepts the long forms of units (i.e. ‘meter’, ‘foot’, ‘second’, etc.), the parser also knows the standard unit symbols (e.g. ‘m’, ‘ft’, ‘s’, ...). So to set units that are a combination of other units, pass them as a string:

sage: speed.in_unit('ft/s')
0.328083989501312 ft/s

The number is the same as before, but the unit is printed as the string passed to the in_unit function.

Similarly, to display the power in units of kWh per day, type:

sage: power.in_unit('kWh/d')
235.359600000000 kWh/d

Now, you may be creative and decide that you want to know the power in horsepower instead. This unit is not pre-defined, so a string version is not going to work either. Instead you have to define the unit yourself. According to wikipedia, one mechanical horsepower is 745.6999 W and the symbol for horsepower is ‘hp’. As this is a physical unit, you can create it by arithmetic. What distinguishes units from other quantities is that they come with a display name telling the program how to print them:

sage: horsepower = 745.6999 * watt
sage: horsepower.set_dispname('hp')   # setting the display name
sage: power.in_unit(horsepower)
13.1509337737607 hp

Dimension checking

Whenever two quantities are added or subtracted, DimPy checks whether that operation is legal:

sage: 4 * meter + 3 * foot   # this is ok because they're both lengths
4.9144 m
sage: 4 * meter - 2 * second   # length - time is illegal
DimensionMismatchError: Subtraction, dimensions were (m, ) (s, )

Quantities can be used in trigonometric functions, exponentials and logarithms, but only if they are dimensionless:

sage: cos(2000 * meter / kilometer)  # this is dimensionless
sage: cos(2 * meter)  # this isn't
DimensionMismatchError: Argument of 'cos' must be dimensionless, dimensions were (m, )

Flydims - non-physical units

Defining new units via arithmetic works fine as long as you only work with physical units. But you may want to find the how much power each inhabitant of Cambridge would get from the weir if the power was distributed evenly. This quantitiy would have the unit ‘power per person’. However a ‘person’ is not a pre-defined unit and it can’t be constructed as a product of physical units.

This is where what DimPy calls ‘Flydims’ come in. Flydims are like physical units, except they are not restricted to the domain of physics. There can be Flydims for ‘person’, ‘pound sterling’, ‘house’, ...

A Flydim is defined by typing name_of_variable = Flydim('name to be printed'). The two names may or may not be the same, but in most cases it will be easiest and cause the least confusion to just use the same word in both cases. Note though that the variable name cannot contain spaces whereas the print name can.


sage: person = Flydim('person')
sage: population = 115200 * person
sage: power_per_person = power / population
sage: power_per_person
0.0851271701388889 W person^-1

DimPy doesn’t have a unit for ‘power per person’, so it prints a combination of SI base units and the Flydim. To make it print ‘kWh/d per person’ you can once again use the QuantityParser (it knows that ‘per’ means division).

sage: power_per_person.in_unit('kWh/d per person')
NameError: name 'person' is not defined

That didn’t work because ‘person’ is not a pre-defined unit, so the parser doesn’t automatically know about it. You have to tell it what you want the word ‘person’ in the string to refer to, using a keyword argument of the form name_in_string=name_of_variable:

sage: power_per_person.in_unit('kWh/d per person', person=person)
 0.00204305208333333 kWh/d per person

Vectors and matrices

DimPy supports vectors and matrices containing quantities, called QuantVector and QuantMatrix.

Say, you want to buy certain amounts of coal, gas and electricity and you are wondering how much energy you will get, what the CO_2 footprint is and how much this will cost.

Now coal is generally traded in kilograms, gas in cubic metres and electricity in kWh. The calorific value of coal is 29.3 GJ/t and that of natural gas is 38 MJ/m^3. The CO_2 emissions associated with the different types of energy production are about 2.4 kg/kg for coal, 2.01 kg/m^3 for gas and, in the UK, 0.58 kg/kWh for electricity. The prices are (very approximate and probably horribly out of date) $100 per short ton for coal, $0.5 per m^3 for gas and $0.12 per kWh for electricity. A short ton is 907.18474 kg.

First some preparation: most of the units are pre-defined, but dollar and short ton are not:

sage: dollar = Flydim('$')
sage: short_ton = 907.18474 * kilogram

Next create the conversion matrix. To create a QuantMatrix from a list of quantities, use the function qmatrix. (It is possible to create QuantMatrices directly, but that requires separating the numbers and units, which would be a lot of work here - see Another example using matrices instead).

sage: matrix = qmatrix([[29.3*10**9*joule/ton, 38*10**6*joule/metre**3, 1*kWh/kWh], [2.4*kilogram/kilogram, 2.01*kilogram/metre**3, 0.58*kilogram/kWh], [100*dollar/short_ton, 0.5*dollar/metre**3, 0.12*dollar/kWh]])
sage: matrix
                     kg^-1 $       m^-3 $     m^-2 kg^-1 s^2 $
m^2 kg s^-2 $^-1   29300000.0    38000000.0         1.0
     kg $^-1           2.4          2.01     1.61111111111e-07
        1        0.110231131092      0.5     3.33333333333e-08

The reason why the QuantMatrix prints like this is to show its structure: the numbers and the units are stored separately and instead of storing a unit for each element, two vectors of units are stored. The unit associated with an element can be found by multiplying the units for its row and column. (This means that not all lists of quantities can be converted into QuantMatrices.)

If instead you want to see the elements as they are, create a PrintMatrix that shows each number with its associated unit:

sage: printmatrix = PrintMatrix(matrix)
sage: printmatrix
[[      29300000.0 Gy      38000000.0 Pa                  1.0                ]
 [           2.4            2.01 m^-3 kg      1.61111111111e-07 m^-2 s^2     ]
 [ 0.110231131092 kg^-1 $    0.5 m^-3 $   3.33333333333e-08 m^-2 kg^-1 s^2 $ ]]

In a PrintMatrix you can change the units, e.g. if you’d prefer the upper left corner to be in MJ/kg, then type:

sage: printmatrix.set_preferred(0, 0, 'MJ/kg')

(the two integers are the indices of the element, starting at 0 as usual in Python).

If you want to set new units for most or all elements of the PrintMatrix, you can do that too by giving the above function a list of lists of units. This list must have the same shape as the matrix. If there are any elements whose units you don’t necessarily want to change, fill up the list with None in those places. Again, don’t forget the keyword argument when using self-defined units.

sage: printmatrix.set_preferred([[None, 'MJ/m^3', 'MJ/kWh'], ['kg/kg', 'kg/m^3', 'kg/kWh'], ['USD/kg', 'USD/m^3', 'USD/kWh']], USD=dollar)
sage: printmatrix
[[       29.3 MJ/kg       38.0 MJ/m^3   3.6 MJ/kWh  ]
 [       2.4 kg/kg        2.01 kg/m^3   0.58 kg/kWh ]
 [ 0.110231131092 USD/kg  0.5 USD/m^3  0.12 USD/kWh ]]

In this case, you cannot make dollar print as ‘$’ because all keywords must be valid Python variable names.

Now if you buy 1 kg of coal, 1 m^3 of natural gas and 1 kWh of electricity:

sage: quantities = QuantVector([1, 1, 1], [kilogram , metre**3, kWh])
sage: matrix * quantities
(70900000.0 J  4.99 kg  0.730231131092 $)

The second and third elements of the vector look good, but the first one would be nicer if it was in kWh. There is no direct access to the quantities in a QuantVector, so you cannot directly assign it a unit like for single quantities. It is possible to convert the QuantVector into a PrintMatrix and then convert to other units, but DimPy also has other ways of dealing with units: if it hasn’t been told a specific unit to print a quantity in, it looks for one in a register and it attempts to find the unit with the most compact numerical representation. Only the SI units and some additional units are registered by default, but you can add or remove units from the registers via register_new_unit(unit) and unregister_unit(unit).

sage: register_new_unit(kWh)
sage: matrix * quantities
(19.6944444444 kWh  4.99 kg  0.730231131092 $)

Here, kWh gives the more compact representation, so it is automatically chosen over joule.

The choice of ‘kWh’ as a unit is now available for all results:

sage: quantities2 = QuantVector([1, 1, 10], [kilogram , metre**3, kWh])
sage: matrix * quantities2
(28.6944444444 kWh  10.21 kg  1.81023113109 $)

Another example using matrices

Consider a pendulum connected to a spring which is fixed to a wall.

|  spring      O | fixed wheel
|                |
                |_| mass

The potential energy of the pendulum depends on the extension x of the spring (a length) and the displacement \theta of the pendulum (an angle).

Assume there is an extremum of the energy function and consider the Taylor expansion about this point. Then the energy relative to the extremum is E = y M y /2 where y is the vector (x, \theta) and M is the matrix of second derivatives.

sage: y = QuantVector([0.01, 0.02], [metre, radian])
sage: y
(0.01 m  0.02)
sage: M = QuantMatrix([[100, 200], [200, 300]], [[joule/metre, joule/radian], [1/metre, 1/radian]])
sage: M
             m^-1    1
 m kg s^-2  100.0  200.0
m^2 kg s^-2 200.0  300.0
sage: 1/2 * y * M * y
0.105 J