数学研发论坛

 找回密码
 欢迎注册
查看: 312|回复: 6

[原创] 一道相当困难的三重积分以及MMA代码

[复制链接]
发表于 2019-5-20 22:14:15 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有帐号?欢迎注册

x
本帖最后由 葡萄糖 于 2019-5-20 23:22 编辑

\[ \Large\iiint\limits_{(x^2+y^2)^2+z^4\,\leqslant\,x+y+z}\mathrm{d}x\mathrm{d}y\mathrm{d}z \]
\begin{align*}
\Large{\left(\sqrt{\left(\sqrt{\color{red}5}-{\color{red}1}\right){\color{red}2}}+\arctan\left(\sqrt{\frac{\sqrt{\color{red}5}+{\color{red}1}}{\color{red}2}}\,\,\right)\,\right)\frac{\pi}{3}}
\end{align*}
对于封闭曲面\(\,\Gamma\colon(x^2+y^2)^2+z^4=x+y+z\,\,\)围成的封闭区域\(\,J\colon(x^2+y^2)^2+z^4\leqslant\,x+y+z\,\,\)
笔者利用Wolfram Mathematica 10.3的Integrate+ImplicitRegion命令还是利用Integrate+Boole命令计算封闭区域\(\,J\,\)的体积,程序都没有给出解析结果
  1. RegionJ = ImplicitRegion[(x^2 + y^2)^2 + z^4 \[LessSlantEqual] x + y + z, {x, y, z}]
  2. Integrate[1, {x, y, z} \[Element] RegionJ]
  3. Integrate[Boole[(x^2 + y^2)^2 + z^4 \[LessSlantEqual] x + y + z], {x, -1, 2}, {y, -1, 2}, {z, -1, 2}]
复制代码

而,利用Wolfram Mathematica 10.3的NIntegrate+ImplicitRegion命令还是利用NIntegrate+Boole命令计算封闭区域\(\,J\,\)的体积,程序给出的数值结果与解析值一致
  1. RegionJ = ImplicitRegion[(x^2 + y^2)^2 + z^4 \[LessSlantEqual] x + y + z, {x, y, z}]
  2. NIntegrate[1, {x, y, z} \[Element] RegionJ]
  3. NIntegrate[Boole[(x^2 + y^2)^2 + z^4 \[LessSlantEqual] x + y + z], {x, -1, 2}, {y, -1, 2}, {z, -1, 2}]
复制代码


************************************************************


另外,将斜着的封闭曲面\(\,\Gamma\colon(x^2+y^2)^2+z^4=x+y+z\,\,\)
绕向量\(\,\overset{\small\rightharpoonup}{l}(-1,1,0)\,\)逆时针旋转\(\,\arccos{\frac{1}{\sqrt{\,3\,}}}\,\),
得到正着的封闭曲面\(\,\Sigma\colon(x+y-z)^4+4(x^2+y^2+z^2-xy+xz+yz)^2=9\sqrt{3}\,z\,\)
  1. Clear["Global`*"];(*Clear all variables*)
  2. (*z轴的单位向量为(0,0,1),待转轴上的一个向量P是(1,1,1),
  3. 计算向量积得到旋转轴上的一个向量q(-1,1,0),计算旋转角度*)
  4. (*绕着向量q(-1,1,0)旋转ArcCos[1/Sqrt[3]],
  5. 所求曲面上的点(a,b,c)旋转后就在斜着的曲面上*)
  6. out = RotationMatrix[ArcCos[1/Sqrt[3]], {-1, 1, 0}].{a, b, c}
  7. (*化简出来的点在方程上,逐项替换,c轴就是z轴*)
  8. xyz = (x^2 + y^2)^2 + z^4 == x + y + z /. Thread[{x, y, z} -> out]
  9. xyz = FullSimplify@xyz
复制代码

其中,封闭曲面\(\,\Gamma\,\)旋转的代码参考该贴:https://bbs.emath.ac.cn/thread-15672-1-1.html

************************************************************

对于封闭曲面\(\,\Sigma\colon(x+y-z)^4+4(x^2+y^2+z^2-xy+xz+yz)^2=9\sqrt{3}\,z\,\)围成的封闭区域
\(\,K\colon(x+y-z)^4+4(x^2+y^2+z^2-xy+xz+yz)^2\leqslant9\sqrt{3}\,z\,\)
由于Wolfram Mathematica 10.3无法解出隐式区域,于是
笔者利用Wolfram Mathematica 10.3的Integrate+ImplicitRegion命令以及Integrate+Boole命令计算封闭区域\(\,K\,\)的体积,程序依旧没有给出解析结果
  1. RegionK = ImplicitRegion[(x + y - z)^4 + 4 (x^2 + y^2 + z^2 - x*y + x*z + y*z)^2 \[LessSlantEqual] 9 Sqrt[3] z, {x, y, z}]
  2. Integrate[1, {x, y, z} \[Element] RegionK]
  3. Integrate[Boole[(x + y - z)^4 + 4 (x^2 + y^2 + z^2 - x*y + x*z + y*z)^2 \[LessSlantEqual] 9 Sqrt[3] z], {x, -1, 2}, {y, -1, 2}, {z, -1, 2}]
复制代码

更令人诧异的是,利用Wolfram Mathematica 10.3的NIntegrate+ImplicitRegion命令以及NIntegrate+Boole命令计算封闭区域\(\,K\,\)的体积,程序给出错误的数值结果
  1. RegionK = ImplicitRegion[(x + y - z)^4 + 4 (x^2 + y^2 + z^2 - x*y + x*z + y*z)^2 \[LessSlantEqual] 9 Sqrt[3] z, {x, y, z}]
  2. NIntegrate[1, {x, y, z} \[Element] RegionK]
  3. NIntegrate[Boole[(x + y - z)^4 + 4 (x^2 + y^2 + z^2 - x*y + x*z + y*z)^2 \[LessSlantEqual] 9 Sqrt[3] z], {x, -1, 2}, {y, -1, 2}, {z, -1, 2}]
复制代码


[AoPS]The volume of the solid enclosed by surface
https://artofproblemsolving.com/community/c7h1835622
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2019-5-21 10:12:38 | 显示全部楼层
解析解不行,那就数值解

点评

你有没有读完这个帖子呀?我的意思是MMA处理不了这个的解析值,而手算可以得出解析值。  发表于 2019-5-21 10:37
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2019-5-21 11:00:22 | 显示全部楼层
写了这么多代码,你应该把代码的允许结果也贴出来,
然后别人看了结果后,再决定是否看你的,
没结果,准备让别人运行代码后再看对你帖子是否感兴趣?
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2019-5-23 13:31:56 | 显示全部楼层
旋转后的曲面好像算的不太对【不应该还是ArcCos[1/Sqrt[3]]】。所以才得出错误的数值结果。
参考https://bbs.emath.ac.cn/thread-15842-3-1.html 的29楼代码。

  1. from={0,0,1};to={1,1,1};
  2. f=Function[{x,y,z},(x^2+y^2)^2+z^4-(x+y+z)];
  3. t=FullSimplify[f@@(RotationTransform[VectorAngle[from,to],Cross[from,to]][{x,y,z}])]
复制代码



===
另外,再怎么旋转,都应该不会比当前的表达式 更简洁的吧
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2019-5-23 16:02:26 | 显示全部楼层
本帖最后由 葡萄糖 于 2019-5-23 16:18 编辑
wayne 发表于 2019-5-23 13:31
旋转后的曲面好像算的不太对【不应该还是ArcCos[1/Sqrt[3]]】。所以才得出错误的数值结果。
参考https://bbs.emath.ac.cn/thread-15842-3-1.html 的29楼代码。 ...


嘻嘻,你算的也是对的!这两个结果是等价的!
不信你看看这两个:
  1. (x + y - z)^4 + 4 (x^2 + y^2 + z^2 - x*y + x*z + y*z)^2 - 9 Sqrt[3] z // Expand

  2.    (5 x^4 - 4 x^3 y + 18 x^2 y^2 - 4 x y^3 + 5 y^4 - 9 Sqrt[3] z +
  3.    4 (x + y) (x^2 - 4 x y + y^2) z + 6 (3 x^2 + 2 x y + 3 y^2) z^2 +
  4.    4 (x + y) z^3 + 5 z^4) // Expand
复制代码

点评

恩,也对。 变换不改变 体积  发表于 2019-5-23 17:47
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

小黑屋|手机版|数学研发网 ( 苏ICP备07505100号 )

GMT+8, 2019-6-26 10:42 , Processed in 0.065851 second(s), 17 queries .

Powered by Discuz! X3.4

© 2001-2017 Comsenz Inc.

快速回复 返回顶部 返回列表