A Volume Problem for EMT

People often ask me how to get started with Euler Math Toolbox (EMT). My answer is always to skim through the basic tutorials and then to solve a problem with the software.

Here is one for you to try:

In the sketch a tank is filled with a fluid up to height h, and the question is simply how much fluid is that? The tank is a horizontal cylinder with two half balls attached.

Of course, you can now go through the internet or a book of geometry to find the two formulas you need to compute this. But if you understand the principle of Fubini you can solve it with paper and pen or EMT. No look-ups are necessary. This principle tells us that we have to integrate the areas of the horizontal cuts through a solid from bottom to top and get the volume. The reasoning for this is that we can imagine the volume as a composition of small shelves stacked upon each other. Integrating is just the limit of taking the sums of the volumes of these shelves.

But how big is the area of the cut in height h? This is not too difficult to compute, if we know the radius of the half circles at the left and right in height h. Let’s make another sketch.

We get

$$x^2 = r^2 – (h-r)^2 = 2rh – h^2$$

With this, it is not very difficult to setup a formula in EMT for the surface at height h.

>function x(h) &= sqrt(2*r*h-h^2)

2
sqrt(2 h r - h )

>function S(h,r,w) &= 2*w*x(h) + pi*x(h)^2

2                   2
2 sqrt(2 h r - h ) w + pi (2 h r - h )

We can get a Latex form of this and paste it into the blog. Let us try.

>tex("S(h,r,w)")
2\,\sqrt{2\,h\,r-h^2}\,w+\pi\,\left(2\,h\,r-h^2\right)

$$S(h,r,w) = 2\,\sqrt{2\,h\,r-h^2}\,w+\pi\,\left(2\,h\,r-h^2\right)$$

The symbolic integral of this turns out to be a lot more verbose.

>function V(h,r,w) &= integrate(S(x,r,w),x,0,h)
Answering "Is  h  positive, negative, or zero?" with "positive"
Answering "Is  r  zero or nonzero?" with "nonzero"

2       r - h        2         r
(- (3 r  asin(-------) - 3 r  asin(-------)) w
mabs(r)              mabs(r)
2            2         3
- (3 r - 3 h) sqrt(2 h r - h ) w + 3 pi h  r - pi h )/3

Maxima is asking questions about the signs of h and r. We could avoid this by making assumptions for h and r, but the automatic answers of EMT look okay. As a test, we can plot the function now for a tank with r=1 and w=2.

>plot2d("V(x,1,2)",0,2):

To make sure, we compare the volume of a full tank with the known geometric formulas for the volume of a cylinder and a sphere. This time, we assume h>0 to get rid of the absolute value of h in the result. Try without this assumption first!

>&assume(h>0); &expand(V(2*h,h,w))

3
2     4 pi h
pi h  w + -------
3

The formula for the Volume at height h is by no means nice.

$$\frac{-\left(3\,r^2\,\arcsin \left(\frac{r-h}{\left| r\right| } \right)-3\,r^2\,\arcsin \left(\frac{r}{\left| r\right| }\right) \right)\,w-\left(3\,r-3\,h\right)\,\sqrt{2\,h\,r-h^2}\,w+3\,\pi\,h^2 \,r-\pi\,h^3}{3}$$

For comparison, we can also numerically integrate the function S(h,r,w). This is necessary anyway in most real world examples. We take the function integrate() of EMT, which uses an adaptive Gauß integration. A simple Gauß scheme would suffice for this function, but we do not care about speed.

>function V1(h,r,w) := integrate("S",0,h;r,w);

To test this, we ask for the necessary height of the fluid if we want to fill the tank with a volume of 7 (whatever units we have). And, indeed, the solution works well when checked with the exact formula.

>h1=solve("V1",1;1,2,y=7)
1.25080577552
>V(h1,1,2)
7

The question remains: How do we get the details of these command right? Indeed, the walk through above seems too complicated to understand at first glance.

First of all, I showed you an exemplary solution. You can make things much easier, if you define r=1 and w=2 globally, and omit the parameters r and w from the functions V, S and the calls to solve() and integrate(). This would look as follows.

>function x(h) &= sqrt(2*r*h-h^2);
>function S(h) &= 2*x(h)*w+pi*x(h)^2;
>&assume(h>0); &assume(w>0); &assume(r>0);
>function V(h) &= integrate(S(x),x,0,h);
>&expand(V(2*r))

3
2     4 pi r
pi r  w + -------
3

>r=1; w=2;
>V(2)
10.471975512
>function V1(h) := integrate("S",0,h);
>V1(2)
10.471975512
>solve("V1",1,y=7), V(%)
1.25080577552
7


Then, you can find the documentation for these functions including examples by double clicking on the functions or pressing F1. You can find numerous examples for the user in the tutorials too.

And there is still the community in the forums which is eager to help.

Diese Website verwendet Akismet, um Spam zu reduzieren. Erfahre mehr darüber, wie deine Kommentardaten verarbeitet werden.