油罐问题——用mathematica计算闭合曲面体积

我们知道,Mathematica是一款无比强大的符号计算软件,有着极为简单粗暴的计算定积分方式:

Integrate[f,{x,xmin,xmax}]

其中xmin,xmax 积分的上下限,那么,用这个函数就可以方便的求得各种规则的体积。

但是在实际建模和应用中,我们常常遇到一些非常复杂的情况,比如油罐容积问题就因为倾斜时候油面形状在不同阶段形状不断变化而需要分段求得,计算量就非常大了。

其实求体积的时候主要困难就是积分表达式无法直接求得,或者积分表达式太过复杂导致运算量太大,在需要求得大量关于油面高度的体积时需要耗费大量时间,甚至有可能根本算不出结果等情况。这个时候就需要用一些比较奇特的思路或者函数来解决问题了。

下面,按照这一次的题目中的不同情况依次说明三种不同的体积求解方式:

1)直接积分法

油罐容积问题中水平放置的时候可以通过直接积分法快速的计算出结果,但是在油罐倾斜情况下,由于油面变化的方式过于复杂,如果继续使用直接积分法,一般需要分为三个阶段,求得各个阶段油面面积关于油面高度h的表达式后进行积分,三个阶段示意图如下:

额,以上。表达式太过复杂,当时求得的结果是什么我忘了,不过想做肯定是可以做出来的。其中,在求油面面积的方程中有一个小Tips,就是倾斜斜面切椭圆柱时得到的切面肯定是一个一个椭圆或者正圆,油面在不同高度的时候其实就是截取这个椭圆不同部分得到的面积。这样求油面面积的表达式比较好写。

可以发现,这样的积分已经是非常复杂,而且计算机也用长时间的沉默(不给出结果)来向我们抗议这种积分方式的计算量,你TM实在是太大啦!那么这次建模的第二问如果继续用这种积分方式,绝对是可行而不可取的方式了。下面我来换一种思路继续求这个油罐的容油量。

2)通过Boole函数来求体积

在此之前先扯点有用没用的知识,Boole是历史上一名著名的数学家(https://en.wikipedia.org/wiki/George_Boole),1864年,布尔死于肺炎。肺炎是他在暴风雨淋湿了之后坚持给学生上课引发的。在这之后,她的妻子坚信以毒攻毒的方式,在布尔病倒之后将他安置在床上倾倒大量凉水后裹好被子,安送了他的丈夫去了西天(这说明有一位无比”疼爱”自己的妻子是多么荣幸的一件事情(●˘˘●))。由于其在符号逻辑运算中的特殊贡献,很多计算机语言中将逻辑运算称为布尔运算,将其结果称为布尔值。突出贡献之一就是将自己的名字赋予了c++中的一个变量(bool)。

Mathematica中也有这个函数,主要功能就继承于c++:判断表达式的正误。

举个栗子:

Boole[True]

给出的结果是1。

通过这个例子,大家想到了什么呢?对啦!就是通过积分Boole函数中的不等式判断某个点是否在区域内,然后对所有满足要求的点进行积分,我们就可以求出一个区域的面积或者体积啦~

栗子2:

In:Integrate[Boole[x^2+y^2≤1],{x,-1,1},{y,-1,1}]

out:Pi

这个式子是怎么求到半径为一的圆的面积呢?看一下Boole函数中的表达式

只有x,y满足上述不等式的时候Boole函数的值是才为1。这个时候我们对这个区域里对Boole函数进行积分,得到的结果就是圆的内部的面积。

类似的,我们可以求由一系列不等式确定的区域体积。以这次题目中油罐为例,油罐内部分的不等式表达式可以表示为:

  1. x^2+z^2<1
  2. y>0&&y<3

额,石怡图如上。下面计算它的体积:

in:vv=x^2+z^2<1&&y>0&&y<3

Integrate[Boole[vv],{x,-Infinity,Infinity},{y,-Infinity,Infinity},{z,-Infinity,Infinity}]

         out:3Pi

事实证明这种思路的正确性,下面我们直接计算本次模型最复杂的情况:

额,示意图如上。下面计算这块的体积:

in:vv=x^2+z^2<1&&y>0&&y<3&&z+0.2y<0.2
Integrate[Boole[vv],{x,-Infinity,Infinity},{y,-Infinity,Infinity},{z,-Infinity,Infinity}]

out:4.122564758856962

计算得到这个结果的时间在1s以内,所以大量取点的情况也不用在意时间的问题~

啊,这个时候可能就有童鞋要问啦,老湿你胡说八道,这怎么能是最复杂的情况呢!最复杂的情况两端是圆的!如下所示:

 

in:in:∫_(-5)^10∫_(-5)^5∫_(-5)^5[Boole[y^2+z^2<2.25&&x>5/8-√(169/64-y^2-z^2 )&&x<59/8+√(169/64-y^2-z^2 )]ⅆzⅆyⅆx]

out:64.66444878638991

呐,结果如上。

其实这种情况也是类似的,那就来示范一下怎么做这个蛋形的体积。

首先得到描述这个蛋的不等式约束条件:

  1. y^2+z^2<2.25(蛋的侧面)
  2. x>5/8-Sqrt[169/64-y^2-z^2](蛋的左底面)
  3. x<59/8+Sqrt[169/64-y^2-z^2](蛋的右底面)
  4. z+xtana+ytanb<h(油面方程)(可选)

额,示意图如上。下面是体积计算代码:

in:vv=y^2+z^2<2.25&&x>5/8-Sqrt[169/64-y^2-z^2]&&x<59/8+Sqrt[169/64-y^2-z^2]&&z+0.2x<1

Integrate[Boole[vv],{x,-Infinity,Infinity},{y,-Infinity,Infinity},{z,-Infinity,Infinity}]

Out:

嗯,你以为我算出来上面这个体积了吗!哈哈哈哈,太年轻啦骚年,

等了足足五分钟啊,尼玛直接给了我一个标准形式的表达式啊,(っ °Д °;)っ

下面是本次体积计算的终极大招,仅限Mathematica 10以上的用户使用。

3)Volume函数(Mathematica 10版本以上可用)

Volume函数顾名思义,其实就是计算体积的意思。注意区分一下内置指示器Volume。

继续举栗子:

in:Volume[Ball[]]

out:4/3Pi

有没有觉得很方便啊!是不是超方便!为湿跑了一天的代码Mathematica 9版本出了很多原来没有出现过的问题啊!坑货啊,听说有中文版本才换上9的啊!最后实在受不了抱着试试看的态度装回了10啊,原来的报错的情况没有了啊!

这个时候只要把我们第二种思路中Boole函数中的表达式用隐函数表达一下放到Volume函数里面就可以求到结果了~

以上。

PS:查资料的时候发现有一篇论文翻译了mathematica中的帮助文档,举了几个计算体积的栗子就水了一篇论文,强烈鄙视这种行为。

发表评论

电子邮件地址不会被公开。 必填项已用*标注