2009年1月20日星期二

Simulating N-Body Problem


理解多体问题的数值模拟

多体问题
物理上,多体问题(N-Body Problem)是指:3维空间中给n个质点在某一时刻的位置和速度,已知他们的质量,假定他们彼此之间只受万有引力作用,求他们后续任意时刻的位置和速度。数学上,多体问题是一个高阶的常微分方程初值问题,具体说来:
设位置函数向量x=(x1,...,xn),,质量mi
满足方程

其分量形式为

和初值条件

求x(t)
当n=2时(又称二体问题),有解析解,解曲线是平面上的二次曲线(《常微分方程教程》P132)。
n>2时目前没有解析解,虽然理论上得到了很多很漂亮的定性分析结果,但是实践中还要归诸数值模拟,也就是用差分代替微分。

多体问题的数值模拟
Part1:差分格式
参考http://en.wikipedia.org/wiki/Verlet_integration
数值模拟最朴素的想法是欧拉迭代,在每一个时刻,根据当前位置计算加速度,更新下一时刻的速度,根据当前速度更新下一时刻位置:(设时间步长为h)
a(t) = f(x(t))
v(t + h) = v(t) + a(t)h
x(t + h) = x(t) + v(t)h
这个方法误差很大,即使假设速度v没有产生误差,也很容易分析出(通过Taylor展开到二阶项),迭代一步产生的误差为O(h2)。

经典力学中的数值模拟常用的是Verlet方法,它利用两侧Taylor展开消去速度项,构造误差阶数更小的差分格式:
首先Taylor展开到4阶

注意到

并将-h带入上面的h,得

两式相加,消去1阶项,顺带把3阶项也消掉了(赚到了,呵呵)。得到:

定义差分格式为

显然一步误差O(h4)
然后来分析误差扩散。
(注:通常Big O Notation是不能做减法的,这里我们不加说明地假设它可以做减法,使用Lagrange余项的Taylor展开也许能解释清楚,我们有深究它。)
来归纳证明,k步误差有如下二次扩散关系

k=1、2时不难直接验证。假设直到k-1都成立,来看k时

其中假设f连续可微,x取值在一个紧集上,就有f是Lipschitz连续的,从而两个h2f项做差是高阶小,直接扔掉。从而上式成立。
因此,若模拟总时间T=kh固定,则总误差O(h2)


Part2:受力计算的分治算法
参考http://www.cs.berkeley.edu/~demmel/cs267/lecture26/lecture26.html
计算每个质点受的合力,按照定义,是其他质点对它的万有引力的矢量和。n个点,每个点计算(n-1)个分力,这样迭代一步的复杂度是O(n2),不能忍。
Barnes-Hut算法使用分治策略,将复杂度降低到O(nlogn),它的思想是,将较远的若干个点当做一个点,一起计算来自他们的分力。举例来说,假设x,y到点集A的质心的距离很远,A包含a个点,则按照定义计算x,y受A的力需要计算2a个分力,而若把A当做一个对象处理,预先计算它的质心和总质量,则完成同样的工作,仅需要a次计算。
算法用到了四叉树的结构(这是对于2维情形,3维使用八叉树)。
递归地将矩形分成四等分,落入第i象限的点属于当前节点的第i子节点,直到当前节点只有一个节点为止,依此从平面上的点集构造四叉树,如图——

通过适当选取划分轴,可以保证四叉树均衡,也称自适应四叉树,这一步复杂度O(nlogn)
得到四叉树后,后续遍历一次得到每个节点的质心和总质量。最后,计算每个质点的受力:
从根节点开始考虑每个子树,如果质点到该子树的根节点的质心足够远,则累加质点受到该子树(作为整体对应的质心和总质量)的万有引力,否则,质点受到该子树的合力等于该子树的每个子树的引力矢量和,而计算子树的引力又递归调用当前算法。由四叉树对应的几何空间划分,很容易看出,计算一个质点的受力,要遍历树上的节点个数不超过树的高度。若四叉树均衡,则复杂度O(nlogn),从而总复杂度也如此。

总结:
1、利用Taylor展开构造有效差分格式可以减少计算误差
2、利用分治思想改进算法可以降低计算复杂度

题外话:
计算数学看算法就看三件事:收敛性、复杂性、稳定性。这个例子中讨论了前两个,第三个稳定性主要是针对良好的问题的,也就是说,如果真实解对初值的依赖是连续的,则要求算法的数值解对初值的依赖也是连续的。但是多体问题是一个典型的混沌系统,问题本身是病态的,讨论它的算法稳定性也就意义不大了。
一直以来,我就特别看好物理的sense好的人,像上面这个分治算法的设计,物理背景就很强(berkeley的讲义上用仙女座星系做四叉树构造思想的导言),或者说,对问题域有透彻了解的人,往往能发明出新的、好的方法。

没有评论: