Ideal gas law with dimensional quantities and a physical constant¶

Maxima is a system for symbolic and numerical computation. Let's use the ezunits package of Maxima to handle dimensional quantities in a problem involving the ideal gas law, pV = nRT.

Inspired by this problem on Stackoverflow: https://stackoverflow.com/questions/64279773/ezunits-not-so-ez-ideal-gas Thanks to user Dux for posing the problem.

In addition to ezunits, we will make use of the physical_constants package, which includes the molar gas constant as %R. Physical constants don't have an ordinary value, so they are carried around as symbols in calculations, until constvalue pulls out the value of the physical constant.

This notebook has been composd in Jupyter using the maxima-jupyter kernel. For more about maxima-jupyter, see: https://github.com/robert-dodier/maxima-jupyter

Problem statement¶

Fuel cell vehicle, molecular hydrogen fuel. 2 tanks, T = 20°C, p_rel = 700 bar, m_tank = 6 kg. What is the volume of one tank? Solve for the tank volume assuming ideal gas law applies.

Solution¶

In [1]:

load (ezunits);

Out[1]:

\[\tag{${\it \%o}_{1}$}\mbox{ /home/robert/maxima/maxima-code/share/ezunits/ezunits.mac }\]

In [2]:

load (physical_constants);

Out[2]:

\[\tag{${\it \%o}_{2}$}\mbox{ /home/robert/maxima/maxima-code/share/ezunits/physical\_constants.mac }\]

Inspect value of molar gas constant %R.

In [3]:

%R, constvalue;

Out[3]:

\[\tag{${\it \%o}_{3}$}\frac{1039309}{125000}\;\frac{\mathrm{J}}{\mathrm{mol}\,\mathrm{K}}\]

Declare unit conversion for bars. 1 bar = 100,000 pascals; this is a little less than one standard atmosphere, which is 101,325 pascals.

In [4]:

declare_unit_conversion (1`bar = 100000`Pa);

Out[4]:

\[\tag{${\it \%o}_{4}$}\mathbf{done}\]

Convert from degrees Celsius to kelvins.

In [5]:

T: 20`degC `` K;

Out[5]:

\[\tag{${\it \%o}_{5}$}\frac{5863}{20}\;\mathrm{K}\]

In [6]:

p_rel: 700`bar;

Out[6]:

\[\tag{${\it \%o}_{6}$}700\;\mathrm{bar}\]

In [7]:

m_tank: 6`kg;

Out[7]:

\[\tag{${\it \%o}_{7}$}6\;\mathrm{kg}\]

In [8]:

p: p_rel+1`bar;

Out[8]:

\[\tag{${\it \%o}_{8}$}701\;\mathrm{bar}\]

Convert from bars to newtons per square meter.

In [9]:

p: p``N/m^2;

Out[9]:

\[\tag{${\it \%o}_{9}$}70100000\;\frac{\mathrm{N}}{\mathrm{m}^2}\]

Molar mass of molecular hydrogen.

In [10]:

M_H_2:2*1.01`g/mol;

Out[10]:

\[\tag{${\it \%o}_{10}$}2.02\;\frac{\mathrm{g}}{\mathrm{mol}}\]

In [11]:

R: %R/M_H_2, constvalue;

Out[11]:

\[\tag{${\it \%o}_{11}$}4.116075247524752\;\frac{\mathrm{J}}{\mathrm{g}\,\mathrm{K}}\]

Solve ideal gas law for volume. dimensionally helps other functions handle dimensional quantities. Maxima prefers exact numbers to inexact, so a floating point value is replaced by a rational approximation.

In [12]:

dimensionally (solve (p*V = m_tank*R*T, V));

rat: replaced -7239.764752871287 by -1223527483/169001 = -7239.764752871285

Out[12]:

\[\tag{${\it \%o}_{12}$}\left[ V=\frac{1223527483}{11846970100000}\;\frac{\mathrm{m}^2\,\mathrm{kg}\,\mathrm{J}}{\mathrm{g}\,\mathrm{N}} \right] \]

In [13]:

V: rhs(%[1]);

Out[13]:

\[\tag{${\it \%o}_{13}$}\frac{1223527483}{11846970100000}\;\frac{\mathrm{m}^2\,\mathrm{kg}\,\mathrm{J}}{\mathrm{g}\,\mathrm{N}}\]

V is the total volume of two tanks. Find the volume of one tank.

In [14]:

V_tank: V/2;

Out[14]:

\[\tag{${\it \%o}_{14}$}\frac{1223527483}{23693940200000}\;\frac{\mathrm{m}^2\,\mathrm{kg}\,\mathrm{J}}{\mathrm{g}\,\mathrm{N}}\]

ezunits does not automatically convert units to basic units. Let's express V in its basic units.

In [15]:

fundamental_units (V_tank);

Out[15]:

\[\tag{${\it \%o}_{15}$}m^3\]

In [16]:

V_tank `` fundamental_units(V_tank);

Out[16]:

\[\tag{${\it \%o}_{16}$}\frac{1223527483}{23693940200}\;\mathrm{m}^3\]

In [17]:

float (%);

Out[17]:

\[\tag{${\it \%o}_{17}$}0.05163883561249133\;\mathrm{m}^3\]

Let's express the volume in more convenient units.

In [18]:

V_tank: % `` liter;

Out[18]:

\[\tag{${\it \%o}_{18}$}51.63883561249133\;\mathrm{liter}\]

We don't really need so many digits in the result.

In [19]:

fpprintprec: 4;

Out[19]:

\[\tag{${\it \%o}_{19}$}4\]

In [20]:

V_tank;

Out[20]:

\[\tag{${\it \%o}_{20}$}51.64\;\mathrm{liter}\]

Great, we'll let that be enough for now. Thanks for reading.

PS. This document was created as a Jupyter notebook using the maxima-jupyter kernel. It was converted to HTML via jupyter nbconvert --to html --template basic ezunits_ideal_gas_law_problem.ipynb and then the MathJax preamble was pasted in by hand.

Solving equations with dimensional quantities¶

Here's a basic example of solving equations in which the variables have units such as meters, kilograms, etc. This example is derived from this question at Stackoverflow. Thanks to user Dux for posting the question.

I'll make use of the add-on package ezunits which implements the representation of units, unit conversions, and other operations, one of which is a function which helps other functions, which don't know anything about units, to be able to handle units. This function is named dimensionally and you can use it by just wrapping it around another function. In this example, we'll apply dimensionally to solve.

I'll start by loading ezunits and then defining the problem data.

In [1]:

load (ezunits);

Out[1]:

\[\tag{${\it \%o}_{1}$}\mbox{ /home/robert/maxima/maxima-code/share/ezunits/ezunits.mac }\]

In [2]:

a: 30*%pi/180;

Out[2]:

\[\tag{${\it \%o}_{2}$}\frac{\pi}{6}\]

Units are indicated with a back quote, e.g. 1000 ` kg means 1000 kilograms.

Note that units are displayed in Roman font -- this is a conventional way to present units.

In [3]:

masse: 1000 ` kg;

Out[3]:

\[\tag{${\it \%o}_{3}$}1000\;\mathrm{kg}\]

I'll write decimals as rational numbers because Maxima is much more comfortable with exact numbers (integers and rationals) than with inexact numbers (floats and bigfloats).

A problem here is that within the context of this problem, g is gravitational acceleration at the Earth's surface. However, in other contexts, g might mean grams. One could work around the problem by naming the gravitational acceleration something else, such as %g, but so long as g is used only one way or the other, we won't run into a problem.

In [4]:

g: 980665/100000 ` m/s^2;

Out[4]:

\[\tag{${\it \%o}_{4}$}\frac{196133}{20000}\;\frac{\mathrm{m}}{\mathrm{s}^2}\]

In [5]:

b: 3/10 ` m;

Out[5]:

\[\tag{${\it \%o}_{5}$}\frac{3}{10}\;\mathrm{m}\]

In [6]:

B: 5/10 ` m;

Out[6]:

\[\tag{${\it \%o}_{6}$}\frac{1}{2}\;\mathrm{m}\]

In [7]:

L: 1/10 ` m;

Out[7]:

\[\tag{${\it \%o}_{7}$}\frac{1}{10}\;\mathrm{m}\]

Arithmetic on dimensional quantities such as masse and g here is carried out according to the usual rules for handling units. In this case, when the quantities are multiplied, the units are multiplied too.

Multiplying mass times an acceleration yields a force. In the Systeme Internationale (metric system), force is measured in newtons, so the result of masse * g could be expressed in newtons. However, ezunits doesn't have a preferred set of units, so it won't convert to newtons automatically. We'll carry out the conversion by hand later on.

In [8]:

F_g: masse * g;

Out[8]:

\[\tag{${\it \%o}_{8}$}\frac{196133}{20}\;\frac{\mathrm{m}\,\mathrm{kg}}{\mathrm{s}^2}\]

In [9]:

F_H: masse * g;

Out[9]:

\[\tag{${\it \%o}_{9}$}\frac{196133}{20}\;\frac{\mathrm{m}\,\mathrm{kg}}{\mathrm{s}^2}\]

Now we're ready to start solving equations. I'll write out each equation so that we can take a look at it, and then solve it. I'll wrap solve inside of dimensionally to help solve handle the units.

In [10]:

eq1: F_H - 2 * S * sin(a) = 0;

Out[10]:

\[\tag{${\it \%o}_{10}$}\frac{196133}{20}\;\frac{\mathrm{m}\,\mathrm{kg}}{\mathrm{s}^2}-S=0\]

In [11]:

solution1: dimensionally (solve (eq1, S));

Out[11]:

\[\tag{${\it \%o}_{11}$}\left[ S=\frac{196133}{20}\;\frac{\mathrm{m}\,\mathrm{kg}}{\mathrm{s}^2} \right] \]

So far, so good -- the result for S looks right. (Of course, it's just a simple equation.) Let's try the second equation for H.

In [12]:

eq2: - F_g + 2 * H;

Out[12]:

\[\tag{${\it \%o}_{12}$}2\,H+\left(-\frac{196133}{20}\right)\;\frac{\mathrm{m}\,\mathrm{kg}}{\mathrm{s}^2}\]

In [13]:

solution2: dimensionally (solve (eq2, H));

Out[13]:

\[\tag{${\it \%o}_{13}$}\left[ H=\frac{196133}{40}\;\frac{\mathrm{m}\,\mathrm{kg}}{\mathrm{s}^2} \right] \]

Third equation for Ly.

In [14]:

eq3: tan(a) = Ly / (B/2);

Out[14]:

\[\tag{${\it \%o}_{14}$}\frac{1}{\sqrt{3}}=4\,{\it Ly}\;\frac{1}{\mathrm{m}}\]

In [15]:

solution3: dimensionally (solve (eq3, Ly));

Out[15]:

\[\tag{${\it \%o}_{15}$}\left[ {\it Ly}=\frac{1}{4\,\sqrt{3}}\;\mathrm{m} \right] \]

The fourth and final equation for FN. This one relates FN to other variables we have already solved for. I'll state the equation first, and then substitute the solutions we found already. (There are other ways of going about it, this is one way that works so I'll just keep going like this.)

In [16]:

eq4: H * B/2 - FN * (L + Ly) + S * sin(a) * B/2 + S * cos(a) * Ly = 0;

Out[16]:

\[\tag{${\it \%o}_{16}$}\frac{\sqrt{3}\,{\it Ly}\,S}{2}+\left(\frac{S}{8}+\frac{H}{4}\right)\;\mathrm{m}-{\it FN}\,\left({\it Ly}+\frac{1}{10}\;\mathrm{m}\right)=0\]

In [17]:

prev_solutions: append (solution1, solution2, solution3);

Out[17]:

\[\tag{${\it \%o}_{17}$}\left[ S=\frac{196133}{20}\;\frac{\mathrm{m}\,\mathrm{kg}}{\mathrm{s}^2} , H=\frac{196133}{40}\;\frac{\mathrm{m}\,\mathrm{kg}}{\mathrm{s}^2} , {\it Ly}=\frac{1}{4\,\sqrt{3}}\;\mathrm{m} \right] \]

In [18]:

eq4a: subst (prev_solutions, eq4);

Out[18]:

\[\tag{${\it \%o}_{18}$}\left(-\left(\frac{1}{4\,\sqrt{3}}+\frac{1}{10}\right)\,{\it FN}\right)\;\mathrm{m}+\frac{588399}{160}\;\frac{\mathrm{m}^2\,\mathrm{kg}}{\mathrm{s}^2}=0\]

Now we're ready to solve for FN.

In [19]:

solution4: dimensionally (solve (eq4a, FN));

Out[19]:

\[\tag{${\it \%o}_{19}$}\left[ {\it FN}=\frac{196133\,3^{\frac{3}{2}}}{16\,\sqrt{3}+40}\;\frac{\mathrm{m}\,\mathrm{kg}}{\mathrm{s}^2} \right] \]

As before, Maxima finds the solution in terms of exact numbers, including the square root of 3. At some point we probably have to go back to decimal numbers, for example if we wanted to tell someone else about the solution. I'll apply the function float to convert the exact numbers to inexact.

In [20]:

float (solution4);

Out[20]:

\[\tag{${\it \%o}_{20}$}\left[ {\it FN}=15050.87322705384\;\frac{\mathrm{m}\,\mathrm{kg}}{\mathrm{s}^2} \right] \]

Finally, let's convert FN to newtons. In ezunits, that's implemented by the double back quote `` operator. Note that the operator distributes over the equal sign -- on the left-hand side, we get FN `` N while on the right we get m*kg/s^2 converted to N.

In [21]:

solution4 `` N;

Out[21]:

\[\tag{${\it \%o}_{21}$}\left[ {\it FN} `` N=\frac{196133\,3^{\frac{3}{2}}}{16\,\sqrt{3}+40}\;\mathrm{N} \right] \]

I'll let that be enough for now. ezunits has other capabilities that I haven't touched on here. If you are interested in ezunits or other aspects of Maxima, I encourage you to bring it up on the Maxima mailing list. Here is the subscription page for the maxima-discuss mailing list.

Thanks for reading and I hope you have a great day.

Probability density of sum of gamma and uniform

The probability density of the sum of gamma and uniform variables is the convolution of their densities, as it was for the sum of Gaussian and uniform variables, which we considered in a previous post. As before I'll define a simplifying unit boxcar function and a general boxcar to represent the uniform density.

matchdeclare (xx, numberp) $

tellsimp (U1(xx), if 0 <= xx and xx < 1 then 1 else 0) $

U(x, a, b) := 1/(b - a) * U1((x - a)/(b - a)) $

Let's let h be the density of a gamma variable, with parameters a and b.

h (x, a, b) := b^a / gamma(a) * x^(a - 1) * exp (-b*x);

⇒

$h(x,a,b):=\left(\frac{b^a}{\Gamma \left(a\right)}{x}^{a-1}\exp \left(\left(-b\right)x\right)\right)$

Here's the convolution integral. I'll quote the density functions to make the structure more obvious and therefore make further manipulation easier.

I : 'integrate ('h (r - s, a[1], b[1]) * 'U (s, a[2], b[2]), s, minf, r);

⇒

${\int }_{-\infty }^{r}h(r-s,a_1,b_1)U(s,a_2,b_2)ds$

As in a previous blog post, I'll define some simplification rules in order to manipulate the integral to the point that Maxima's built-in integrate function can handle it. There's a slightly obscure strategem in the definition of the second rule. yy_minus_xx is assigned before turning off the simplification flag, since it is important that the rule match the simplified form of y - x, which is internally represented as y + (-1)*x.

matchdeclare (xx, symbolp, [aa, bb], all) $

matchdeclare (ee, all) $

matchdeclare ([ll1, ll2], all) $

yy_minus_xx : yy - xx $

simp : false $

tellsimp

('integrate (ee * 'U(xx, aa, bb), xx, ll1, ll2),

1/(bb - aa) * 'integrate (ee, xx, max (ll1, aa), min (ll2, bb))) $

tellsimp

('integrate ('U(''yy_minus_xx, aa, bb)*ee, xx, ll1, ll2),

'integrate ('U(xx, yy - bb, yy - aa)*ee, xx, ll1, ll2)) $

simp : true $

Now let's allow Maxima to apply the rules we just defined. The effect of quote-quote is just as if we just typed in the value of I.

I1 : ''I;

⇒

$\frac{{\int }_{a_2}^{\min (b_2,r)}h(r-s,a_1,b_1)ds}{b_2-a_2}$

I originally defined the integral with h quoted in the integrand. Let's verbify it.

I2 : I, h;

⇒

$\frac{b_1^{a_1}{\int }_{a_2}^{\min (b_2,r)}{\left(r-s\right)}^{a_1-1}{e}^{-b_1\left(r-s\right)}ds}{\left(b_2-a_2\right)\Gamma \left(a_1\right)}$

I'll make a change of variables from r - s to u which simplifies the integrand somewhat.

I3 : changevar (I2, u - (r - s), u, s);

⇒

$-\frac{b_1^{a_1}{\int }_{r-a_2}^{r-\min (b_2,r)}{u}^{a_1-1}{e}^{-b_{1}u}du}{\left(b_2-a_2\right)\Gamma \left(a_1\right)}$

Now we're ready to evaluate the integral. I'll assert some inequalities to forestall questions from 'integrate'. These are all necessary conditions for the definition of the gamma and uniform density functions.

assume (a[1] > 0, b[1] > 0, b[2] > a[2]) $

assume (r > a[2]) $

assume (min (b[2], r) > a[2]) $

I4 : I3, nouns;

⇒

$-\frac{b_1^{a_1}\left(\frac{\Gamma (a_1,b_{1}r-b_1{a}_2)}{b_1^{a_1}}-\frac{\Gamma (a_1,b_{1}r-b_1\min (b_2,r))}{b_1^{a_1}}\right)}{\left(b_2-a_2\right)\Gamma \left(a_1\right)}$

Terrific, now we have computed the convolution symbolically. Note the gamma function, Γ(a), and the upper incomplete gamma function, Γ(a, z), in the result. We'll take a look at a specific example, as given by specific values for the parameters.

I5 : subst ([a[1] = 1/4, b[1] = 1/2, a[2] = 5, b[2] = 12], I4);

⇒

$-\frac{2^{\frac{1}{4}}\Gamma (\frac{1}{4},\frac{r}{2}-\frac{5}{2})-2^{\frac{1}{4}}\Gamma (\frac{1}{4},\frac{r}{2}-\frac{\min (12,r)}{2})}{72^{\frac{1}{4}}\Gamma \left(\frac{1}{4}\right)}$

Let's define a function which has those parameter values baked in, and plot it.

conv_symbolic (r) := ''I5 $

plot2d (conv_symbolic, [r, 0, 20]) $

It's interesting that the convolution function is asymmetric. I guess that's not surprising, since the gamma density is asymmetric too. Well, that's all for now.

Probability density of sum of Gaussian and uniform

The probability density function of the sum of a Gaussian variable and a uniform variable is the convolution of their densities. Let's compute that convolution symbolically, and compare it with a numerical approximation.

Constructing the symbolic convolution

First let's define the Gaussian density g and the uniform density U. These depend on some parameters.

g(x, a, b) := 1/(b*sqrt(2*%pi))*exp(-(1/2)*(x - a)^2/b^2) $

U(x, a, b) := 1/(b - a) * U1((x - a)/(b - a)) $

I've defined the general uniform density in terms of a unit uniform density U1 defined on the interval [0, 1). (Maxima doesn't have a built-in notion of a unit step or boxcar function so I'll do this part myself.) I'll make U1 a simplifying function in order to prevent it from being evaluated except when the argument is a number.

matchdeclare (xx, numberp) $

tellsimp (U1(xx), if 0 <= xx and xx < 1 then 1 else 0) $

Now let's construct the convolution integral. In general the convolution conv(f, g) = integral f(r - s) g(s) ds. I'll quote g and U so that they are not evaluated; that makes it easier to work with the integral in this case.

I : integrate ('g(r - s, a[1], b[1]) * 'U(s, a[2], b[2]), s, minf, inf);

⇒

${\int }_{-\infty }^{\infty }g(r-s,a_1,b_1)U(s,a_2,b_2)ds$

Simplifying and evaluating the convolution integral

The convolution integrand contains a function which is nonzero only on an interval, and constant there, so we can pull out the constant factor and change the limits of integration. I'll define a simplification rule for that.

matchdeclare (xx, symbolp, [aa, bb], all) $

matchdeclare (ee, all) $

matchdeclare ([ll1, ll2], all) $

simp : false $

tellsimp ('integrate (ee * 'U(xx, aa, bb), xx, ll1, ll2), 1/(bb - aa) * 'integrate (ee, xx, max (ll1, aa), min (ll2, bb))) $

simp : true $

First step toward evaluation is to let Maxima apply the simplification rule I just defined.

I1 : ''I;

⇒

$\frac{{\int }_{a_2}^{b_2}g(r-s,a_1,b_1)ds}{b_2-a_2}$

Terrific, now it's just a Gaussian density in the integrand. It's a noun (i.e., quoted function), so verbify it. Then 'integrate' can compute its integral.

I2 : I1, nouns;

⇒

$\frac{\frac{\sqrt{\pi }{b}_1\mathrm{erf}\left(\frac{\sqrt{2}r-\sqrt{2}{a}_2-\sqrt{2}{a}_1}{2b_1}\right)}{\sqrt{2}}-\frac{\sqrt{\pi }{b}_1\mathrm{erf}\left(\frac{\sqrt{2}r-\sqrt{2}{b}_2-\sqrt{2}{a}_1}{2b_1}\right)}{\sqrt{2}}}{\sqrt{2}\sqrt{\pi }{b}_1\left(b_2-a_2\right)}$

Rearrange it somewhat, just to make it nicer to look at.

I3 : ratsimp (I2);

⇒

$-\frac{\mathrm{erf}\left(\frac{\sqrt{2}r-\sqrt{2}{b}_2-\sqrt{2}{a}_1}{2b_1}\right)-\mathrm{erf}\left(\frac{\sqrt{2}r-\sqrt{2}{a}_2-\sqrt{2}{a}_1}{2b_1}\right)}{2b_2-2a_2}$

I'll pick out some typical values for the parameters, and bake them into a function to evaluate the convolution.

I4 : subst ([a[1] = -6, b[1] = 1/2, a[2] = 2 , b[2] = 9], I3);

conv_symbolic (r) := ''I4 $

⇒

$-\frac{\mathrm{erf}\left(\sqrt{2}r-3\sqrt{2}\right)-\mathrm{erf}\left(\sqrt{2}r+2^{\frac{5}{2}}\right)}{14}$

Let's make a plot. The convolution of a Gaussian and a uniform looks like a softened boxcar. Each end of the boxcar is a shifted and rescaled erf.

plot2d (conv_symbolic, [r, -10, 10]) $

Numerical approximation to the convolution via FFT

Here I'll just paste in the discrete convolution via FFT which was the subject of a recent blog post.

load (fft) $

conv (a, b) :=

block

([n : length(a) + length(b) - 1],

inverse_fft (fft (zero_extend (a, next2(n))) * fft (zero_extend (b, next2(n)))),

realpart (next2(n) * %%),

rest (%%, n - next2(n))) $

zero_extend (a, n) := makelist (if i <= length(a) then a[i] else 0, i, 1, n) $

next2 (n) := block ([i : ilog2 (n)], if n = 2^i then n else 2^(i + 1)) $

ilog2 (n) := block ([i : 0], while n > 1 do (n : floor (n/2), i : i + 1), i) $

Comparing the symbolic and numerical results

I'll create a grid of values and evaluate the Gaussian and uniform densities on the grid. The parameters for the densities are the same as what were baked into conv_symbolic.

n : 100 $

dx : (10 - (-10))/n $

xx : makelist (-10 + (i - 1) * dx, i, 1, n), numer $

gx : dx * makelist (g (xi, a[1], b[1]), xi, xx), a[1] = -6, b[1] = 1/2, numer $

Ux : dx * makelist (U (xi, a[2], b[2]), xi, xx), a[2] = 2, b[2] = 9, numer $

Then I'll compute the convolution of those discretized densities.

conv_numerical : conv (gx, Ux) $

conv_numerical_density : conv_numerical / dx $

conv_support : makelist (-20 + (i - 1) * dx, i, 1, 2*n - 1) $

Finally I'll plot both the discretized approximation and the symbolic density. The agreement is pretty close. That's comforting.

plot2d ([[discrete, conv_support, conv_numerical_density], conv_symbolic], [r, -10, 10])$