利用龙格库塔法求解郎之万方程

一、待解问题

布朗颗粒是非常微小的宏观颗粒，其直径的典型大小为10-7~10-6m。颗粒不断受到液体介质分子的碰撞，在任一瞬间，一个颗粒受到介质分子从各个方向的碰撞作用力一般来说是互不平衡的，颗粒就顺着净作用力的方向运动，由于分子运动的无规性,施加在颗粒上的净作用力涨落不定，力的方向和大小都不断发生变化，颗粒就不停地进行着无规则的运动，描述这一过程的理论称为朗之万理论.

dr设颗粒的质量为m，在时刻t颗粒的坐标为r(t）,颗粒受到粘滞阻力m，其

dt中m为粘滞系数;还受到一个周围颗粒给它的涨落力Rm(t),涨落力随时间t变化，可正可负,其平均值为0；此外颗粒还会受一个势力颗粒运动的朗之万方程为：

U。根据牛顿第二定律，rd2rdrUm2mRm(t) dtdtt对于这样的二阶微分方程，通常是化为一阶微分方程组来求解，由于

d2rdvdvdrdvv，得到 dt2dtdrdtdrdrdtdv1U(mvRm(t))dtmt

v在这个化简的朗之万方程中，我们期望得到在某一时刻位置和速度的关系。

二、物理机理

1.龙格-库塔法的基本思想及一般形式

设初值问题yy(x)[a,b]，由微分中值定理可知，必存在[xn,xn1]，使

y(xn1)y(xn)hy()y(xn)hf(,y())设

yny(xn),并记K*f(,y()),则

y(xn1)ynhK*

其中K*称为y(x)在[xn,xn1]上的平均斜率，只要对平均斜率K*提供一种算法,上式就给出了一种数值解公式，例如，用K1f(xn,yn)代替K*，就得到欧拉公式，用K2f(xn1,yn1)代替K*，就得到向后欧拉公式，如果用K1,K2的平均值来代替K*，则可得到二阶精度的梯形公式。可以设想，如果在[xn,xn1]上能多预测几个点的斜率值，用它们的加权平均值代替K*，就有望的到具有较高精度的数值解公式，这就是龙格—库塔法的基本思想。 龙格-库塔公式的一般形式:

yn1ynhcikii1r

K1f(xn,yn)Kif(xnih,ynhijKj)j1i1 （1）

其中Ki是yy(x)在xnih(0i1)点的斜率预测值；ci,i,ij均为常数，选取这些常数的原则是使(1）式有尽可能高的精度。

2、四阶龙格—库塔法

四阶龙格—库塔法的表达形式如下：

yn1ynh(c1K1c2K2c3K3c4K4)K1f(xn,yn)K2f(xn2h,ynh21K1)K3f(xn3h,ynh31K1h32K2)K4f(xn4h,ynh41K1h42K2h43K3) （2）

其中有13个待定系数，我们希望适当的选取这些系数使公式的精度尽可能

高，先将K2,K3,K4按照如下二元函数泰勒级数展开,取到h4项，再将K2,K3,K4展开后的值代入（2）式得到yn1的表达式.

i11kf(xnih,ynhijK1)(ihhij)f(xn,yn)...xyj1k0k!j1

i1n再用泰勒公式将y(xn1)展开，取到h5项，得到y(xn1)

由于局部截断误差Rn1y(xn1)yn1，泰勒级数的系数匹配,使局部截断误差为O(h5),得到

,,c12123132342324324 4142434,c1c2c3c4111222ccc,ccc 2133442233442311 33c23cc,cc()233443232424234346 1122c32332c44(242343),c32c()2324242343812

得到11个方程和13个未知量，因此需要补充两个条件2,310，求解得到：

123,41,21,32,410,420,11221111431,c1,c2,c3,c4633612

最后得出经典的四阶龙格-库塔公式：

hyn1yn(K12K22K3K4)6K1f(xn,yn)hhK2f(xn,ynK1)22hhK3f(xn,ynK2)22K4f(xnh,ynhK3)

三、数值计算方案

下面利用四阶龙格库塔法来求解朗之万方程，在前面我们已经推算出一阶微分方程组:

vdrdt dv1U(mvRm(t))dtmt取初值为(t0,r0,v0)，时间t的步长为h,得到：

hvn1vn(K12K22K3K4)6K1f(tn,vn)hhK2f(tn,vnK1)22hhK3f(tn,vnK2)22K4f(tnh,vnhK3)

hrn1rn(K12K22K3K4)6K1f(tn,rn)hhK2f(tn,rnK1)22hhK3f(tn,rnK2)22K4f(tnh,rnhK3)

由于已知初值为(t0,r0,v0)，可以推算出在t0h时刻的r1,v1的大小,由此不停的迭代，可以知道在任意t0nh时刻的rn,vn的大小，进而可以分析颗粒运动的特性。

四、流程图

输入t0,r0,v0,取t的步长为h tnt0nhhrn1rn(K12K22K3K4)6vn1vnh(K12K22K3K4)6tn,rn,vn