xzcyr xzcyr
关注数: 3 粉丝数: 743 发帖数: 13,773 关注贴吧数: 18
【交流】这个二维N-S方程边值问题在其他软件里是怎么解的? 标题又写不下了……其实原本我拟的标题是 【交(gǎo)流(shì)】这个经典(大概)的纳维斯托克斯方程的边值问题(二维稳态不可压缩流)在其他软件里是怎么解的? 大家都知道,NDSolve自版本10开始加入了有限元方法("FiniteElement"),终于可以求解 1. 不含时的偏微分方程定解问题(或者说得更“数学”一点:不管在哪个方向上都不是 适定的初值问题的偏微分方程定解问题) 2. 不规则区域上的的偏微分方程定解问题 了,尽管目前NDSolve的有限元方法展现出了一些吸引人的特性(比如简洁易懂的语法),但是,它还很年轻, 截止版本11.2,它还不能求解非线性问题。比如这一帖 mathematica.stackexchange.com/a/96579/1871 里面提及的问题。如大家所见,user21 使用了一些较低级的有限元函数来求解这个问题,老实说,我一个用了几年Mathematica的人,看了他的代码都觉得脑仁疼……显然,在当前,如果需要在不规则区域上求解非线性的偏微分方程,Mathematica恐怕不是一个很好的选择。 那什么才称得上是“很好的选择”呢? 开这个帖就是想让各位来秀一下。(我知道吧里存在复数个其他有限元软件的用户。) 事先声明下规则:所有讨论围绕上面所举的例子进行,谢绝空谈,违者……可能删。(毕竟这是个涉及软件对比的帖子,我得最大限度地压制住出现无谓争端的可能性。) 为了方便大家阅读,这里重述一下上面帖子中出现的问题。首先是方程组:参数: \[Mu] = 10^-3; \[Rho] = 1; 区域(就直接上代码了,应该都能看懂吧……): l = 2.2; h = 0.41; \[CapitalOmega] = RegionDifference[Rectangle[{0, 0}, {l, h}], ImplicitRegion[(x - 1/5)^2 + (y - 1/5)^2 < (1/20)^2, {x, y}]]; RegionPlot[\[CapitalOmega], AspectRatio -> Automatic]边界条件 \[CapitalGamma] = { DirichletCondition[p[x, y] == 0., x == l], DirichletCondition[{u[x, y] == 4*0.3*y*(h - y)/h^2, v[x, y] == 0}, x == 0], DirichletCondition[{u[x, y] == 0., v[x, y] == 0.}, y == 0 || y == h || (x - 1/5)^2 + (y - 1/5)^2 <= (1/20)^2]}; ……其实就是左边界速度分布为4*0.3*y*(h - y)/h^2,压力的右边界为0。 好了,大家来展示下自己所用的软件是怎么解这个问题的吧。
3分钟复现系数为正负1的全体n次幂方程的根的复平面图形 原本是想直接续在这帖:http://tieba.baidu.com/p/2134509239的后面直接回复的,但是考虑到那帖楼层数已经较多,默默发在那边无法满足我的虚荣心(殴)所以决定重开一帖。(其实@隨意超 似乎已经在那帖11楼找对了方向,但却似乎没有引起后续的参与者的重视。) 先上代码和图再讲解吧。2GB内存2GHz老爷机,3 min出图: modifiedroots[c_List] := Block[{a = DiagonalMatrix[ConstantArray[1., Length@c - 1], -1]}, a[[1]] = -c; Eigenvalues[a]]; s = 1000; l = ConstantArray[0., {s + 1, s + 1}]; l[[#, #2]] += 1. & @@@ Round[1 + s Rescale@ Function[z, {Im@z, Re@z}, Listable]@ Flatten[modifiedroots /@ Tuples[{-1., 1.}, 18]]]; // AbsoluteTiming ArrayPlot[UnitStep[80 - l] l, ColorFunction -> "AvocadoColors"]当然,我说“复现”其实是标题党了,因为原图算到了24次多项式,我这幅图只算到了19次。(这是2G内存的极限,次数每加1,内存占用大约要翻一番。)但是如大家所见,这幅图的效果已经和那幅Sam Derbyshire画了3天3夜的图相差无几了。 下面简单讲解一下这段代码。要做到在多数人可以 接受的时间内使用多数人 可以接受的内存获得 和那幅图相似的效果,需要: 1 省略对这个特别的问题而言其实不必要的(主要是符号)计算。 2 昨晚对像素点十分有限的屏幕而言其实不必要的仅有着细微差别的数据。 3 正确地注意到原图的颜色代表的是什么。 先说第一点。因为这幅图表示的是幂方程的根,所以我们自然需要求解幂方程。如何解方程?我们的第一反应是Solve。是的,Solve可以求解这个问题,但是,作为一个通用求解器,它对这个问题而言太慢了。有没有更为专用的函数可以求解这个问题?在本文开头的帖子里@草红样 找到了NRoots,一个专解多项式方程的函数,从而大大提速了求解。可是,这是否是终点?不,NRoots依旧需要进行多项式结构的分析,我们能不能把这一部也给省掉,直接使用(以确定次序排列的)多项式系数表作为输入来求根?再一次,答案是肯定的,那就是我们上面使用的modifiedroots。(对于NRoots的提速问题的讨论其实尚未结束,有兴趣的参考这帖:http://tieba.baidu.com/p/3575212088) 另,不难注意到“全体系数为正负1的全体n次幂方程”其实有一半是重复的。上面的modifiedroots已经利用了这个特性。 再说后两点。(这两点其实是相联系的,所以一起讲。)对方程的求解将会产生大量的数据,这倒是其次,真正够呛的是如何将这些数据绘成图形。(具体参看本文开头帖子中大伙儿的哀嚎。)正常情况下,要画这么多的点很花时间,还吃内存,更麻烦的是,还怎么画也不像Sam Derbyshire那幅“震撼”的图!要跨过这个槛,需要另一个技巧,好吧,其实这个技巧我是从这帖学来的:http://tieba.baidu.com/mo/q/checkurl?url=http%3A%2F%2Fmathematica.stackexchange.com%2Fquestions%2F50839%2Fsmooth-peter-de-jong-attractor%23comment149926_50839&urlrefer=826a9b85381fccc015cc5edd4575e535 (注意下面的答案最好也看看),这里的重点是,Sam那幅图中的颜色明暗到底表现的是什么?复平面上同一点处重合的点量?我们不难想像完全重合的点其实是很少的,是的, 颜色的明暗其实表现的只是某一个点附近小区域内的点量。如何决定这个小区域?对于这点,我们的计算机屏幕其实已经帮我们限制好了: 屏幕上的像素点是很有限的。于是这又反过来启发了我们:既然最终在屏幕上呈现的点数有限,我们又何必把所有的点都拿去画图?只要事先统计好每个像素内有多少点,再把这个统计结果拿去画图就好了嘛。哪个函数可以经济地画出这种原始数据点数极多的明暗图?我们有ArrayPlot。ArrayPlot可以调配色方案嗯。 另,最后画图时使用UnitStep剪掉了实轴上的点,因为影响视觉效果。 感觉其实还有些地方地方可以解说下不过现在也累了所以算了。 最后,其实即便抛开方程阶数,单论配色方案,我也没有完全“复现”人家的图。(小声:不过我自以为我选的方案更好看。)有谁有兴趣可以研究下怎么用Blend把原图的颜色混合出来。(应该不是某个Mathematica里现成的配色方案吧……)
关于NRoots函数的(可能)提速,以及某matlab函数roots的测试请 为了避免一会儿答昏了头又忘了先把这帖发出来。会折腾起这个问题是因为前两天在知乎上看到了一个问题(那个问题现在貌似被删了……)。该问题问的是如何提升NRoots的效率,并声称NRoots函数比某matlab的roots()函数(定义在这里http://tieba.baidu.com/mo/q/checkurl?url=http%3A%2F%2Fradio.feld.cvut.cz%2Fmatlab%2Ftoolbox%2Fmathlib%2Fcppmathug%2Fmatlab11.html&urlrefer=c718dc5954add794b75faa368c11b46c)效率低一个数量级。 虽然我的电脑上没有装matlab,但我认为这一“声称”还是可信的,因为NRoots是使用多项式(而非列表)作为输入的,那么几乎可以肯定,对符号的分析是会占用时间的。类似的情况我们可以在LinearSolve上看到:在解一次方程组时它的速度大致比Solve快一个数量级。 然后我就在想,NRoots会不会像Solve之于LinearSolve一样,拥有一个抽象的(我的意思是,直接以列表为输入的)对应函数呢?在帮助里找了一下,找不到,想了想就把问题扔SE去了,然后,被踩了……路过的Daniel Lichtblau猜是不是因为我没给出具体的测试结果,我也觉得有可能,于是……有人能帮忙测一下然后截个图什么的吗? 问题在这里:http://tieba.baidu.com/mo/q/checkurl?url=http%3A%2F%2Fmathematica.stackexchange.com%2Fq%2F73312%2F1871&urlrefer=2e2bc58f5501205c6d2be4a3e2a39fba 要是能帮我直接编辑进问题里那就更感谢了! 要是能顺便再Upvote一下那就更更感谢了! 直接在这帖下面讨论这个问题也是欢迎的,不过要是能帮忙到那边把问题给答了那就更更更感谢了! 好了,先这样。我继续屠版去了。
我是怎样学起Mathematica来 本文系自此帖(http://tieba.baidu.com/p/2736346260)5楼的回复修订增补而来。吧里的老人们阅读本文时可能会觉得各种眼熟:毕竟重要的经验我都是不断强调的。 我最初会得知Mathematica的存在,是因为《数学建模》老师在课堂上的推荐。我第一眼就被它那极度接近传统数学式书写的使用方法给吸引了——相信这个吧里的许多人也是一样。《数学建模》结课后不久校内就搞了数学建模大赛(目的是为了选拔参加全国大赛的选手),我就参加了,用的当然是Mathematica。 注意这时的时间点是快三年前,也就是2012年的上半年,那个时候具有全中文帮助文档的Mathematica 8已经出世,但是,版本8的体积比之前的版本要大,于是我这个吝啬鬼想了半天最后 错误地选择了版本7,并从图书馆借了一本《Mathematica7实用教程》,觉得这样应该也够用了。 比赛时间是三天,我所属的组在最后一天算是弄出了模型,然后,解不出来。 站在现在看来,这一局面是理所当然的,因为我那时对Mathematica的使用,几乎是胡来,因为Mathematica的代码接近传统数学式,我便 想当然的在Mathematica里写起了传统数学式。我还能记得的在这一时期犯过的错误——实际犯的肯定不止这些——有: 为了获得纯正的传统数学表达,大量依赖TraditionalForm和自由格式输入。 以为在代码的开头写上 a ∈ Integers 就可以使 a 被全局声明为整数。 不会清变量,对语法着色也毫无概念,只好通过反复重开软件来使软件“恢复正常”。 ………… 总而言之,那时的我对 Mathematica其实是一种严谨的编程语言这点可以说是完全没有认识,而那本 万恶的《Mathematica7实用教程》也基本上只是把Mathematica当成一个计算器在介绍,这便进一步封死了我认知自己错误的路。参赛的结果,也就不消我多说了吧。这次失败让我陷入了强烈的沮丧之中,从此一厥不振,Mathematica成为了我心中又一段不堪回首的记忆——要真是这样我现在就不会在这儿了,故事当然还有后续。 这次经历固然让我受到了打击,但也让我感到了不甘心,结果,Mathematica还是留在了我的电脑里没有删掉。不久,我所在的教研室迎来了一个组会日,我那时还没进实验室,手上没什么东西可讲,搜肠刮肚,我忽地想起自己本科毕设时的一个问题可以用Mathematica来解决,就拿这个当材料做PPT了。 注意,这个问题是连当时的我也能自力解决的问题,不消说,它是十分简单的。(实际上就是用Reduce求解整数不定方程。)但是我却忘记了自己当初是怎么被形如Solve[x^3 + 2 x^2 + 3 x + 4 == 0, x]的代码迷得神魂颠倒,不知道这一切在旁人眼里是多么的神奇,没想到自己已经被选为了目标……暑假临近时,导师忽然找到了我,让我帮忙解一个数学模型。我在这个时候才得知,教研室内部不知何时起已经流传起了一个奇怪的传闻,说是我参加建模大赛获奖了囧。 “既然导师信任我,那我就不能辜负人家的期望”——这种假话我当然说不出口,但是我拼命想要快点解完这个模型那是千真万确的,因为,我想回家。毕竟这是我大学以来第一次暑假没家回,虽然理性上清楚这对研究生来说是正常的,但是感性上这就好像一直拥有的权利突然被剥夺了一样浑身难受。我所要求解的模型实际上是一个偏微分方程组,我一看到这个问题自然想起了NDSolve(《实用教程》介绍了这个函数)。“那就赶快把方程塞进这个指令里解出来然后回家!”我第一时间 把电脑里的Mathematica换成了版本8——我在这个时候已经意识到了纸质教程的内容恐怕会不足,而且在回家的诱惑面前硬盘空间当然是浮云——输入方程,Shift+Enter,我得到的是——一大片的警告。 但是此刻我已经不是几个月前的我了!几个月前的失败也就是浪费几天时间,这次可是关乎回不回得了家的大事!……而且,或许是不断的失败使我在潜移默化中有了反省,最终我竟基本冷静了下来。我开始意识到我面临的其实是复数个问题的集合, 贪婪地想要一口气把问题彻底解决是不现实的。于是我开始读着自带帮助和网上看到的各色杂七杂八的资料开始了跌跌撞撞的探索。基本掌握Mathematica的最基础语法并没花掉我很多时间,可是,还是不够……我最终开始在百度知道上提问(这算是Mathematica学习再开后的第一问:http://tieba.baidu.com/mo/q/checkurl?url=http%3A%2F%2Fzhidao.baidu.com%2Fquestion%2F447269535.html&urlrefer=4efec66f5a57e52c755730b0ed6bf5fd)。 读到这里,或许会有读者对我此时的选择嗤之以鼻或是感觉不理解:一个学术问题,怎么想到跑百度知道上来问了?为什么不问问你的同学和老师呢?比如你的《数学建模》老师?关于这点,一个原因是我的身边没有什么人会Mathematica,而《数学建模》与我甚至不同院,当时还是暑假,我怎么好意思去麻烦别人——好吧,这其实是借口,说白了我就是从小怕老师怕惯了,所以对找老师这事有抵触——不过,现在看来,即使我当时选择这条路,等在我前方的恐怕依旧是死胡同,说句听起来有些不知天高地厚的话:我丝毫不认为在课堂上教我们破解版本4(还是版本3?记不起来了)Mathematica的老师能够给我帮什么忙——当然了,也可能他会因此为我介绍一位新的老师,最终闯出一条新路……遐想就到此为止,历史不容假设,而且我丝毫不后悔我当时的选择。 虽然困难,但我开始将我要求解的问题一点一点的分解开,放到知道上提问,虽然得到的答案并不总是如意,但是还是有收获的,至少我开始对我所面临的问题的难度以及其中重大的难点看得越来越清楚了—— 百度知道并不像某些人想的那么糟糕,还是有人在认真的回答问题的。这里要特别感谢 @tryytr ,尽管他对我所要解的问题并不熟悉,却依旧在这段我最为艰难的时期积极地与我进行了多次探讨。(比如这里:http://tieba.baidu.com/mo/q/checkurl?url=http%3A%2F%2Fzhidao.baidu.com%2Fquestion%2F452380013.html&urlrefer=e9800c58323aab1db01a1b2ee478212d 和这里:http://tieba.baidu.com/mo/q/checkurl?url=http%3A%2F%2Fzhidao.baidu.com%2Fquestion%2F455225843.html&urlrefer=b7e6b068687046b7a06b2bd2662a9303 )但是,依旧不够……有的问题不解决我还是回不了家……就在我再度陷入僵局,甚至在漫天飞的“Matlab更擅长数值计算”的 流言影响下装上了Matlab的时候,神( @yang_bigarm )出现了!:(这个提问因为最终含了外链被百度和谐了,所以这里只贴个图。) 得知 Stackexchange的存在(由stackoverflow找到mathematica的分站去并没花掉我多少时间),可以说是我学习Mathematica过程中最重要的 转折点——今天看来,我所要求解的模型中所蕴含的问题,其实是超乎想象的,其中的许多内容至今没有彻底解决。(鉴于其复杂度以及吧里关心偏微分方程问题的人员稀少度,此处按下不表。)但是,mathematica.stackexchange的成员们的实力,同样超乎想象。通过在Stackexchange的提问,我所面临的主要问题均得到了解决。(这是其中最重要的一个:http://tieba.baidu.com/mo/q/checkurl?url=http%3A%2F%2Fmathematica.stackexchange.com%2Fq%2F9277%2F1871&urlrefer=cb49ff4e68e2a11a79b0a665418d1039 ,我至今认为这是我在Stackexchange的最有价值的提问。)——而这,依旧是次要的,更重要的是,在我在Stackexchange阅读帖子的过程中, 我对Mathematica的热情再度迸发了出来。 诚然,我在首次接触Mathematica时也曾为之着迷,但那更多的是被它的光鲜一面给吸引,没法长久,而这一次,我是在幻想破灭之后, 被Mathematica真正的强大吸引的。同时,对于一些基本原则的掌握也使我有了进一步学习的资本。(比如 勤查自带帮助这一点——这也是stackexchange的各位强人们在不断强调的。)我再次迎来了一个学习Mathematica的高潮。 我一面 在Stackexchange大量阅读帖子,一面开始在百度知道疯狂答题(之所以那个时候没在Stackexchange答,一来是因为本身实力不济,二来是因为那时候Stackexchange开站不久,原始会员们答题远比现在积极,我根本抓不到机会),那时我的答题近乎饥不择食(话说那段日子里抢了 @yang_bigarm 不少“生意”真是不好意思……) ,甚至连要注册码的问题我也来者不拒(http://tieba.baidu.com/mo/q/checkurl?url=http%3A%2F%2Fzhidao.baidu.com%2Fquestion%2F492396412.html&urlrefer=82fc9008bf83e97dae09e607ce261125),现在看来,这种做法无疑造成了一定的 负面影响,那就是这种无差别答题给我造成了巨大的精神压力,以致我在某个时期回答他人问题时的态度很不好(啥,现在也不咋的?……总之比那时候好多了吧。这都是靠着我如今对答题时间的限制),但是,不可否认的是尝试 回答他人的问题确实是一个让自己快速进步的方法。 我在这一时期的进步是非常快的,毫不夸张地说,我每阅读一个Stackexchange的帖子,每回答一个百度知道的问题,我都能切实地感觉到自己比上一刻更懂Mathematica。到了大约十月的时候,我已经打遍百度知道无敌手了(我不否认这和百度知道上的提问总体容易有关,其实在我在知道上开始答题后的很长一段时间里,甚至连Table都不能熟练使用:http://tieba.baidu.com/mo/q/checkurl?url=http%3A%2F%2Fzhidao.baidu.com%2Fquestion%2F475944231.html&urlrefer=4fedc0aa0e4f13347436644fe195e9b3)。同时百度知道有的时候还是会冒出一些困难得光靠自查帮助解决不了的问题的,这种时候我就会把问题搬到Stackexchange去——我就是靠着这个到手了我在stackexchange的第一个nice question:http://tieba.baidu.com/mo/q/checkurl?url=http%3A%2F%2Fzhidao.baidu.com%2Fquestion%2F499080092.html&urlrefer=6d164ecb644233c5e9392c8661637db7 http://tieba.baidu.com/mo/q/checkurl?url=http%3A%2F%2Fmathematica.stackexchange.com%2Fq%2F15021%2F1871&urlrefer=62029c8278b3e49483380533c9c7a4af 。顺便这一时期我最重要的回答是:http://tieba.baidu.com/mo/q/checkurl?url=http%3A%2F%2Fzhidao.baidu.com%2Fquestion%2F505315229.html&urlrefer=270afee2707bc9201f907f56ceef926c 。 再往后的事,就波澜不惊了。先是把答题的战场扩张到了贴吧。(我曾因为一些历史原因,犹豫了很长时间是否要进驻贴吧,但是知道上的问题已日益无法满足我的需求,最后,“管他呢!”顺便这大概是吧内第一帖:http://tieba.baidu.com/p/1950670750)然后又开始 阅读Leonid Shifrin写的那本书。(我已经记不起来最初是因为什么才开始读这本书的了,不过,长时间在Stackexchange摸爬滚打,接触到这本书绝对是时间问题。)大概到了年底,也就是正式开始学习Mathematica约半年的时候,我发觉自己的进步速度似乎慢了下来——其实可学的东西还有非常非常多(http://tieba.baidu.com/p/2025917205),但是对软件的学习毕竟受限于自己的知识和兴趣——再往后,我的Mathematica技能的进步进入了一个比较平缓的时期,经常很长时间都没有什么显著的进步,当然,绝不是没有进步,实际上,包括编译、模式匹配在内的许多重要内容我都是在这一时期学会的。(关于模式匹配,必须要感谢 @33珊珊2013 于2013年3月前后在知道的一系列问题的促进作用。) 以上,就是我学习Mathematica的主要过程的一个啰啰嗦嗦的回顾。如今的我虽然不再是个新手,但是和Stackexchange的众神相比还是有着光年级的差距,即使是在吧里,我也绝不是最强的。(不知道全吧打排位赛的话能不能进前十啊……)我在学习过程中走了不少弯路,某些做法也不见得对,不过,若是各位能从这个帖子里吸取到些许经验或教训,我会非常高兴。 凡是在本帖下面仅回复“太长不看”或类似词句者,删。 最后的最后,依旧用从前在 @青衣瓦屋 的帖子下面胡诌的小诗做结吧: 假如StandardForm欺骗了你, 不要悲伤,不要心急! 出错的日子里需要镇静, 相信吧,习惯InputForm的日子将会来临。
【吧务】关于精品帖门槛的放低 嗯……细心的同学可能发现最近一个月吧里精品帖的数量多了起来,很多的老帖突然被加了精。 会产生这个想法的直接的原因是,过去的一个月里的某几天,突然又出现了数位求问诸如“你们的Mathematica是怎么学的的”同学……为什么又有人问?这个我们不是说过很多次么?教程的地址也贴了很多次,为什么没有人看?难道……于是我试着新建了个“教程”的精品分类,然后把@mm_酱 那个翻译Leonid的书的帖子给加了精,结果……这个沉了不知多少年的帖子果然被立刻挖出来了!明明过去我们不管怎么自顶,都会迅速沉掉的! 看样子过去吧里加精的门槛设的太高了,没能活用这个制度,发挥出它的力量。 大家应该也有感觉,那就是百度贴吧的搜索功能越来越渣了,有时候,有的内容,明明过去讨论过,却怎么搜也搜不到。如果是在我自己发的主题里那还好,如果是别人的帖,甚至只是在我自己过去的回复里,那找起来简直就是大海捞针,类似的事我已经做了多次,不想再做了。我相信我绝不是唯一一个这么想的人。 而且连我们自己也搜不到的东西,别人怎么能搜到呢? 如果别人搜不到、答完了就废,那我坚持不加Mathematica相关的QQ和QQ群、一直在公开场合答题的意义又在哪呢? 为什么不把这种需要多次引用的帖子放到精品里呢? 就因为它们中有价值的往往只有某一层楼? 就因为帖子本身没有整理的四平八稳和个论文似的? 就因为问题本身其实很简单,答题者没有很麻烦很累就把问题解决了? 没必要弄这种形式主义。 今后,只要是感觉值得一看的帖子都会考虑加精,不论帖子中实际蕴含的工作量是高还是低,也不论帖子内容是否“整齐”,更不论帖子的可看之处是在主楼还是在回复或是在楼中楼。 我自己能记忆起来的(以及能找到的)帖子毕竟有限,大家看到或想到有意思的帖子时大可以挖上来,我会加精的。 对于这个“政策”有什么意见或建议的话,欢迎回帖讨论。
首页 1 2 3 下一页