Blog List View

冲击摆的力学模型,让我来做力学板块第一篇入门教程吧 

2015-04-08 17:32:48
3502 次阅读


    力学里面有很多神奇的模型,可能随便给出的条件和特别简单的构件就能有一个神奇的运动,画出很多漂亮的运动轨迹,今天便借新科学网新开的机会发一篇mma在力学中的简单应用。

先来一个很简单的模型,单摆。

对于单摆我们都熟知这样的一个模型,有这样的微分方程 [math] heta ''[t] + frac{g}{L} sin heta[t]= 0[/math],这样的微分方程是没有解析解的。所以而我们通常用小量近似[math]sin( heta)sim heta[/math]来将其转变成可解析的微分方程 [math] heta ''[t] + frac{g}{L} heta[t]= 0[/math],然后很容易就得出了[math] heta[t]= heta _{m}sin(omega t+varphi )[/math]

但是我们明确的知道这个近似要求角度小于5°,这个条件是近似的依据,所以我们先来看看近似和不近似的差别:因为mma强大的解微分方程能力,我们来做这样的一个例子,也是对mma的求微分方程功能做一个简单介绍:

首先定义g和L两个基本参数,然后sol1是真实情况下的摆动,即未作近似的解,sol2是小量近似后的结果,二者很快就显示出了差异。这个例子就很好的反应了mma 的几个特性:1.支持输入各种字符。2.变量不用声明,直接就可以用。3.图片质量高。4.代码简单。

在这个例子中我们可以看到,sol1是用了NDSolve这个函数,其中N代表数值运算,D代表微分,Solve代表解方程,合起来就是数值求解微分方程,他的调用形式就是这样:NDSolve[{微分方程及初始条件},{待求解函数列表},{变量及其范围}]运行这个返回的就是一个插值函数的替换规则。而sol2用的是DSolve,这个函数是用来求微分方程的解析解,调用方式和NDSolve很类似,只是变量不用指定范围。这里值得注意的是,在写等式的时候一定要用双等号即==表示等式,不然就成为了赋值语句,很容易出错。还有一点:数值求解只在指定范围内成立,超出范围的解是不可信的。然后我们再用Plot函数把这两个图画出来,其中的Evaluate函数是用来在plot之前计算具体的函数的,即将sol1和sol2这两个替代规则应用到函数上,而替换正是mma的精华。怎么样,是不是很简单,也很神奇?

然后你就可能问了,mma是不是能解决更加复杂的问题呢?答案是肯定的,mma具有强大而自动化程度极高的微分求解系统,下面举一个复杂的例子来展示这一点:

冲击摆,和上述的简单单摆不同,这个模型摆的上端不是固定的,而是可以自由滑动的,下端垂直,然后突然给下端一个冲击,使其获得一个初始速度,这种题目很显然用普通的力学方法求解会很麻烦,这里采用拉格朗日方法求解。条件是这样的:g=9.8m/s^2,绳长L=1m,小球质量为1KG视为质心,上方小球静止,下方小球有初速度(水平方向)1m/s。显然这个系统有两个广义坐标,一个是x一个是[math] heta[/math],即上端小球的横坐标和之间连线与竖直方向的角度。由拉格朗日方程我们有:[math]frac{mathrm{d} }{mathrm{d} t}(frac{partial L}{partial dot{q_{i}}})=frac{partial L}{partial q_{i}}[/math] 只要写出拉格朗日函数L,就可以对该系统进行求解,下面给出mma解法:

先输入条件,g,L,m,其中取x1和[math] heta[/math]为广义坐标,所以用他们表示出来下面小球的位置(x,y),mma的符号计算引擎极为强大,完全可以直接用坐标的导数来写动能表达式,从例子中也可以看出这一点,然后很容易的写出拉格朗日函数,即图中的L。关于前面这一部分的写法要特别提一下:mma里面表示函数的方法有4种,第一是最常见的 f[x],第二种是 f@x ,第三种是后置式 x//f,第四种是针对二元函数的,形式为a~f~b,这里采用了第二种形式。求导也很简单,D[函数,变量 ,次数]。

点击下载附件

然后就到了拉格朗日方程的部分了,我们知道对第i个广义坐标,拉格朗日方程的形式都是一样的,我们可以写一个通式[math]frac{mathrm{d} }{mathrm{d} t}(frac{partial L}{partial dot{q_{i}}})=frac{partial L}{partial q_{i}}[/math] ,只要把其中的[math]q_{i}[/math]替换成x1和[math] heta[/math]就行,所以就有了图中的写法,即mma中的匿名函数,其中#代表第一个变量,#2代表第二个变量,#n代表第n个变量,然后在函数最后加上&便可以了,例如要定义一个f(x)=x^2,我们可以写成 f[x_]:=x^2,也可以直接写成 #^2&,明白了这一点以后,我们就可以像上面那样写一个匿名函数D[D[#, #2'@t], t] - D[#, #2@t] == 0& ,只要我们将这个函数中的变量分别替换为L和广义坐标就行了,于是就有神奇的@@@,帮助里是这么解释的 f @@@ {{a, b, c}, {d, e}}={f[a, b, c], f[d, e]},所以[D[D[#, #2'@t], t] - D[#, #2@t] == 0] & @@@ {{L, x1}, {L, [math] heta[/math]}}这句话的意思就是两个拉格朗日方程啦。之后就按部就班,利用NDSolve对其求解,结果保存在S中。

有了S也就有了我们所求的x1和[math] heta[/math],只要用替换规则就可以得到相应的函数,我们可以看看底端小球的运行轨迹是什么样子,这里x和y都是t的函数,应该会想到用参数方程画图,当然,mma有这样的函数,ParametricPlot,具体用法查帮助吧,mma的帮助是最好的。然后我们有了这样的图,是不是很像一本翻开的书?

点击下载附件

有人在意我为什么刚好要取1.44为结束点吗?这可不是随便取的哦,看起来这就是一个周期,我们可以对它做一个验证,我们可以画出x1[t]和x1’[t]的图像,理论上说,当x1又一次速度为0时,就应该是一个周期。先画x1[t]和x1‘[t]图:

点击下载附件

从图中我们大致可以知道一个周期在1.4左右,为了准确的求出这个周期。我们可以用Findroot来求。就像这样:

In[155]:= FindRoot[Evaluate[x1'@t/.S]==0,{t,1.4}]

Out[155]= {t->1.44163}

所以周期就是1.44163秒啦。。是不是很神奇~~

竟然巴拉巴拉讲了这么多,不管啦,这是我在新科学的第一篇文章哦~快来用mma来解决你们的问题吧~




发表评论