罗兰谈MATLAB的艺术

将想法转化为MATLAB

请注意

罗兰谈MATLAB的艺术已退役,不会更新。

利用符号数学工具箱计算面积惯性矩

我很高兴再一次介绍客座博主Kai Gehrs。在过去的五年里,Kai一直是MathWorks的软件工程师,致力于符号数学工具箱.他有数学和计算机科学的背景。他已经为我的博客贡献了过去的写作在MATLAB中使用符号方程和符号函数以及方法简化符号结果

内容

简而言之:这篇文章是关于什么的?

如果你有兴趣使用MATLAB和符号数学工具箱来教授一些机械工程的基础知识,这可能会对你感兴趣。

计算面积惯性矩是力学中的一项重要任务。例如,面积惯性矩在结构的应力、动力和稳定性分析中起着至关重要的作用。在本文中,我们使用符号数学工具箱计算椭圆管截面的面积弯矩。

我们从只涉及数字参数的基本情况开始,然后通过引入符号参数使计算更一般。

本文中使用的所有图都是在MuPAD笔记本应用程序.您可以在本文的末尾找到这些图的源代码。

基本示例:椭圆管的横截面

下图显示了椭圆管的横截面。

下面的椭圆定义了部分的外部和内部轮廓:

  • $ y = \下午2 \ \√6{1 - \压裂{x ^ 2}{9}} $描述等高线$x \in[-3,3]$。
  • $y = \pm \√{1 - \frac{x^2}{4}}$描述内心的等高线$x \in[-2,2]$。

我们将计算这部分关于x轴的面积转动惯量。

面积惯性矩

面积A相对于x轴的转动惯量用二重积分来定义

$$I_x = {\int \int}_A y^2\, \math {d}A.$$

你可以在维基百科

这个例子背后的数学

在我们的例子中,区域$A$是前面图中显示的孵化区域。利用关于$x$-和$y$-轴的截面的对称性,我们可以将计算限制在$x-y$-平面的第一象限。

为了在孵出的区域上计算必要的二重积分,我们将这个区域分为两个独立的区域$A1$和$A2$。

为了计算关于x轴的最终面积转动惯量,我们用$A_1$和$A_2$的二重积分之和乘以4。

$ $ I_x = {\ int \ int} _A y ^ 2 \ \ mathrm dA = 4 \ cdot \大({\ int \ int} _ {A_1}
y ^ 2 \ \ mathrm dA_1 + {\ int \ int} _ {a₂}y ^ 2 \ \ mathrm dA_2 \大),$ $

现在,我们只需要计算这两个二重积分。多维积分背后的数学理论告诉我们,我们可以用两个积分来重写这些二重积分:

$ $ I_1 = {\ int \ int} _ {A_1} y ^ 2 \ \ mathrm dA_1 = \ int_0 ^ 2 \ int_{\√{1 -
\压裂{x ^ 2}{4}}} ^{2 \ \√6{1 - \压裂{x ^ 2} {9}}} y ^ 2 \ \ mathrm dy \,
\ mathrm dx, $ $

$ $ I_2 = {\ int \ int} _ {a₂}y ^ 2 \ \ mathrm dA_2 = \ int_2 ^ 3 \ int_{0} ^{2 \ \√6{1 - \压裂{x ^ 2} {9}}}
Y ^2\, \mathrm dy\, \mathrm dx $$

int函数从符号数学工具箱可以计算这些积分。

象征性的集成

我们将$x$和$y$定义为符号变量:

信谊xy

我们开始计算$I_1$的内积分$\int_{\sqrt{1 - \frac{x^2}{4}}}^{2 \, \sqrt{1 - \frac{x^2}{9}}} y^2 \, \mathrm dy$:

innerI1 = int (y ^ 2, y, sqrt (1 - x ^ 2/4), 2 *√(1 - x ^ 2/9))
innerI1 = (8 * (9 - x ^ 2) ^ (3/2)) / 81 - (4 - x ^ 2) ^ (3/2) / 24

接下来,我们进行积分innerI1关于x,从0到2。这提供了$I_1$:

I1 = int(innerI1, x, 0,2)
I1 = 3*asin(2/3) - /8 + (74*5^(1/2))/81

我们用同样的策略来计算$I_2$。首先,计算内积分$\int_{0}^{2 \, \sqrt{1 - \frac{x^2}{9}}} y^2\, \mathrm dy$。然后,我们将得到的表达式对$x$从$2$到$3$进行积分:

I2 = int (int y ^ 2, y, 0, 2 *√(1 - x ^ 2/9)), x, 2, 3);

因此,椭圆管相对于$x$轴的面积惯性矩为

Ix = 4 * (I1 + I2)
Ix = (11*pi)/2

我们可以在数值上近似地使用vpa函数。例如,我们用5位有效(非零)数字来近似计算结果:

vpa(第九,5)
Ans = 17.279

高级示例:由符号参数定义的椭圆管截面

我们可以使用数值积分器,比如MATLABintegral2,计算上例中的面积惯性矩。数值积分器可能会返回稍微不太准确的结果,但除此之外,使用符号积分器并没有太大的好处。

但是,如果我们考虑一个形状由任意符号参数定义的更一般的椭圆管的截面呢?

例如,我们要计算这个椭圆管的面积转动惯量。

它的等高线仍然是椭圆,但它们由符号参数$r_1$, $r_2$, $r_1$和$r_2$定义:

  • $y = \pm R_2 \, \sqrt{1 - \frac{x^2}{R_1^2}}$描述$x \in [-R_1,R_1]$的外轮廓线。
  • $y = \pm r_2 \, \sqrt{1 - \frac{x^2}{r_1^2}}$描述$x \in [-r_1,r_1]$的内等高线。

现在我们需要计算这些积分:

$ $ I_1 = {\ int \ int} _ {A_1} y ^ 2 \ \ mathrm dA_1 = \ int_0 ^ {r_1} \ int_ {r_2 \ \√6 {1 -
\压裂{x ^ 2} {r_1 ^ 2}}} ^ {R_2, \ \√6{1 - \压裂{x ^ 2} {r_1 ^ 2}}} y ^ 2 \ \ mathrm dy \,
\ mathrm dx, $ $

$ $ I_2 = {\ int \ int} _ {a₂}y ^ 2 \ \ mathrm dA_2 = \ int_ {r_1} ^ {r_1}
\ int_ {0} ^ {R_2 \ \√6{1 - \压裂{x ^ 2} {R_1 ^ 2}}} y ^ 2 \ \ mathrm dy \ \ mathrm dx。$ $

虽然现在的情况更加复杂,但我们仍然可以使用与上面相同的策略,使用符号变量而不是数字。尽管如此,我们还需要增加一个步骤:指定符号参数之间的关系。这是因为变量$r_1$, $r_2$, $r_1$和$r_2$可以是任何复数,除非我们明确限制它们的值。例如,系统不知道这些变量是正还是负,是$r_1 < r_1 $还是其他。在这个例子中,我们想要指定$r_1>0$, $r_2>0$, $r_1> r_1$和$r_2> r_2$。在符号数学工具箱中,这样的限制被称为关于变量的假设.属性可以设置假设而且assumeAlso功能:

信谊xyr1r2R1R2;假设(r1 > 0);假设(r2 > 0);假设(r1 < r1);假设(r2 < r2);I1 = int (int y ^ 2, y, r2 *值(1 - x ^ 2 / r1 ^ 2), r2 *值(1 - x ^ 2 / r1 ^ 2)), x, 0, r1);I2 = int (int y ^ 2, y, 0, R2 *值(1 - x ^ 2 / R1 ^ 2)), x, R1, R1);Ixgeneral = 4 * (I1 + I2);漂亮的(Ixgeneral)
/ 4 / r1 \ \ / / r1 \ \ | 3 r1最佳 | -- | | | 4 3个R1 asin | - - - | | 3 | \ R1 / | 3 | 3πR1、R2 R1 / | 4 | # 1  + ---------------- | 4 R2 | # 1  - -------- + ---------------- | \ 8 / \ 16 8  / ------------------------------- - ------------------------------------------ 3 3 3 R1 3 R1 3πR1, r2  - --------- 4, / 2 3 \ | 5 R1 R1 R1 | 2 2 1/2 # 1  == | -------- - --- | ( R1 - R1) \ 8 4 /

现在,让我们用第一个例子中相同的值替换符号参数。我们得到了相同的结果吗?让我们再检查一下:当替换$r_1=2,r_2,=1, r_1= 3, r_2 =2$ in时Ixgeneral,我们是否得到了相同的结果9获得了最初考虑的部分?

和预期的一样,通解Ixgeneral还原到具体的解9当我们代入原始维度时:

vpa(潜艇(Ixgeneral (r1, r1, r2, r2)[2、3、1、2]),5)
Ans = 17.279

请注意,必须对$r_1$, $r_2$, $r_1$和$r_2$设置假设。默认情况下,“符号数学工具箱”假定包括符号参数在内的所有变量都是任意复数。这使得计算积分更加复杂和耗时。例如,尝试计算Ixgeneral不需要对$r_1$, $r_2$, $r_1$和$r_2$进行假设。一般来说,如果没有任何额外的假设,甚至不可能得到结果。

如何创建本文中显示的MuPAD图形

要生成本文中显示的图形,请使用MuPAD笔记本应用程序.请注意,本节中的代码不会在MATLAB命令窗口中运行。

打开MuPAD Notebook应用程序:

  1. 在MATLAB工具条上,单击“应用程序”选项卡。
  2. 在“应用程序”选项卡上,单击“应用程序”部分末尾的向下箭头。
  3. 在数学,统计和优化下,单击MuPAD笔记本按钮。

另外,类型mupad在MATLAB命令窗口中。

将以下命令粘贴到MuPAD笔记本上,即可获得上述四种图形:

椭圆:= (x,a,b) ->平方根(b²*(1-x²/a²)):
f1: =情节::Function2d(椭圆(x 2 1), x = 3 . . 3): f2: =情节::Function2d(椭圆(x 3 2), x = 3 . . 3): f3: =情节::Function2d(椭圆(x 3 2), x = 3 . . 2): f4: =情节::Function2d(椭圆(x 3 2), x = 2…3):
f5: =情节::Function2d(椭圆(x 2 1), x = 3 . . 3): f6: =情节::Function2d(椭圆(x 3 2), x = 3 . . 3): f7: =情节::Function2d(椭圆(x 3 2), x = 3 . . 2): f8: =情节::Function2d(椭圆(x 3 2), x = 2…3):
f9:=图::函数2d(椭圆(x,2,1), x = 0..2,VerticalAsymptotesVisible = FALSE): f10 := plot::Function2d(Ellipse(x,3,2), x = 0..3, VerticalAsymptotesVisible = FALSE): f11 := plot::Function2d(Ellipse(x,3,2), x = 2..3, VerticalAsymptotesVisible = FALSE):
plot(plot::Hatch(f1,f2),plot::Hatch(f3),plot::Hatch(f4),f1,f2, plot::Hatch(f5,f6),plot::Hatch(f7),plot::Hatch(f8),f5,f6,缩放=受限,GridVisible = TRUE, VerticalAsymptotesVisible = FALSE,高度= 120,宽度= 200,头部= "椭圆管截面(x-y平面)"):
plot(plot::Hatch(f9,f10),plot::Hatch(f11),f9,f10,缩放=受限,GridVisible = TRUE,高度= 120,宽度= 200,标题= "限制到第一象限(x-y平面)"):
plot(plot::Hatch(f9,f10,Color=RGB::Magenta),plot::Hatch(f11,Color =RGB:: Green),f9,f10,缩放=受限,GridVisible = TRUE, plot::Text2d("A1",[1.0,1.25]), plot::Text2d("A2",[2.5,0.25]),高度= 120,宽度= 200,标题= "集成区域"):
情节(情节::孵化(f1、f2)剧情::孵化(f3)剧情::孵化(f4), f1, f2,情节::孵化(f5、f6)剧情::孵化(f7)剧情::孵化(f8), f5、f6,情节::Text2d (r1,[1.8, 0.05]),情节::Text2d (r2,[0.05, 0.85]),情节::Text2d (r1,[2.74, 0.05]),情节::Text2d (r2,[0.05, 1.84]),比例=约束,GridVisible = TRUE, VerticalAsymptotesVisible = FALSE, XTicksLabelsVisible = FALSE, YTicksLabelsVisible = FALSE, ViewingBox =(-3.5, 3.5, -2.5, 2.5),身高= 120,宽度= 200,标题=“椭圆管截面(更一般的情况)”):

你试过符号数学工具箱吗?

你在与力学相关的计算中使用过符号数学工具箱吗?让我知道在这里

MATLAB®R2013a发布

|
  • 打印
  • 发送电子邮件