利用龙格库塔法求解郎

利用龙格-库塔法求解朗之万方程一、待解问题布朗颗粒是非常微小的宏观颗粒，其直径的典型大小为10-7~10-6m。颗粒不断受到液体介质分子的碰撞，在任一瞬间，一个颗粒受到介质分子从各个方向的碰撞作用力一般来说是互不平衡的，颗粒就顺着净作用力的方向运动，由于分子运动的无规性，施加在颗粒上的净作用力涨落不定，力的方向和大小都不断发生变化，颗粒就不停地进行着无规则的运动，描述这一过程的理论称为朗之万理论。设颗粒的质量为m,在时刻t颗粒的坐标为r(t)，颗粒受到粘滞阻力-ym*,其中Ym为粘滞系数;还受到一个周围颗粒给它的涨落力Rm(t)，涨落力随时间t变化，可正可负，其平均值为0;此外颗粒还会受一个势力-竺。根据牛顿第二定律，dr颗粒运动的朗之万方程为：TOC\o"1-5"\h\z\o"CurrentDocument"d2r dr dUF—md+Rm(t)—诙对于这样的二阶微分方程，通常是化为一阶微分方程组来求解，由于\o"CurrentDocument"d2r dvdvdr dv= = - =v，得到\o"CurrentDocument"dt2 dtdrdt drdrv dtdv1 au、遍=m("V+Rm(t)F)在这个化简的朗之万方程中，我们期望得到在某一时刻位置和速度的关系。二、物理机理1.龙格•库塔法的基本思想及一般形式设初值问题y=y(x)e［a,b］，由微分中值定理可知，必存在&e［x,x「，使y(x)=y(x)+hy&)=y(x)+hf(&,y(&))n+1 n n设yn=y(xn)，并记k*=f(&,y(&))，则其中K*称为y(x)在［xx］上的平均斜率，只要对平均斜率K*提供一种算法，n,n+1上式就给出了一种数值解公式，例如，用%=f(xn,yn)代替K*，就得到欧拉公

式，用K2=f(x打yn1)代替K*,就得到向后欧拉公式，如果用K1,K2的平均值来代替K*，则可得到二阶精度的梯形公式。可以设想，如果在［七,七」上能多预测几个点的斜率值，用它们的加权平均值代替K•，就有望的到具有较高精度的数值解公式，这就是龙格-库塔法的基本思想。龙格-库塔公式的一般形式：1y=y+h'^Lckn+1 n iii=1(1)K1=f(xn,yn)(1)K=f(x+人h,y+h尤日K)i nin ijjj=1其中K是y=y(x)在x+人h(0电＜1)点的斜率预测值；i ni i选取这些常数的原则是使(1)式有尽可能高的精度。2、四阶龙格-库塔法四阶龙格-库塔法的表达形式如下：Iy=y+h(cK+cK+cK+cK)TOC\o"1-5"\h\zn+1 n1122 33 44K1=f(x^,yn)K=f(x+Xh,y+h日K)K=f(x+Xh,y+h日K+h口K)3 n 3 n 31 1 32 2(2)K=f(x+Xh,y+h口K+h口K+h口K)(2)4 n 4 n 41 1 42 2 43 3其中有13个待定系数，我们希望适当的选取这些系数使公式的精度尽可能高，先将K2,K3,K4按照如下二元函数泰勒级数展开，取到h4项，再将K2,K3,K4展开后的值代入(2)式得到七+1的表达式。f(xf(x+Xh,y+h切RK)=£!(Xh_i+h如R£)kf(x,y)+...

k!k=0nin ij1j=1再用泰勒公式将y(x)展开，取到h5项，

n+1由于局部截断误差Rn+1=变+1)一七+18、，、.8 1R,, )kf(x，y)+…j=1得到y(x)n+1泰勒级数的系数匹配，使局部截断误，差为O(h5)，得到日=N,日+日=N,c人日日TOC\o"1-5"\h\z21 2 31 32 3 423243 24日+日+日=人,c+c+c+c=141 42 43 4 1 2 3 4c人+c人+c人=—,cN+c人2+cN21 33 44 222 33 44 3cN+cN+cN=L,c人日+c(人日+人日)=-22 33 44 4 3232 4 242 343 61112c人人日+c人(人日 +人日)=,c如日+c(如日+如日)=一1232332 44 242 343 8 3232 4 242 343 12=0=0，求解得到11个方程和13个未知量，因此需要补充两个条件气=2,四31到：人=!，人=1，R=!，R=!，R=°，R=°,TOC\o"1-5"\h\z3 2 4 21 2 32 2 41 421 1111U=1,c=—,c=—,c=—,c=—43 1 62 33 34 6最后得出经典的四阶龙格-库塔公式：h, - 、y=y+_(K+2K+2K+K)

n+1n61 2 3 4%=f(%,yn)hhK=f(x+,y+/)n2n2 1hhK=f(X+-,y+-K)n2n2 2K=f(x+h,y+hK)三、数值计算方案下面利用四阶龙格库塔法来求解朗之万方程，在前面我们已经推算出一阶微分方程组：dr[v=—dtdv1, ,d,八如、石=m㈠卜+R(t)—石)取初值为(t°,r,v°)，时间t的步长为h，得到:

h,h, 、v1=v+^CK1+2K2+2K3+K4)K1=fCt,,vn)hh、K=f(t+,v+K)n2， 21hhK=f(t+,v+K)」n2，22K4=f(t+h,v+hK3)hr=r+-(K+2K+2K+K)

n+1n6 1 2 3 4K1=f(\,r「hhK=f(t+—,r+-K)」n2n2 1hhK=f(t+—,r+-