《利用龙格库塔法求解郎》由会员分享，可在线阅读，更多相关《利用龙格库塔法求解郎（5页珍藏版）》请在金锄头文库上搜索。
1、利用龙格-库塔法求解朗之万方程一、待解问题布朗颗粒是非常微小的宏观颗粒，其直径的典型大小为10-710-6m。颗粒不 断受到液体介质分子的碰撞，在任一瞬间，一个颗粒受到介质分子从各个方向的 碰撞作用力一般来说是互不平衡的，颗粒就顺着净作用力的方向运动，由于分子 运动的无规性，施加在颗粒上的净作用力涨落不定，力的方向和大小都不断发生 变化，颗粒就不停地进行着无规则的运动，描述这一过程的理论称为朗之万理论。 设颗粒的质量为m,在时刻t颗粒的坐标为r(t),颗粒受到粘滞阻力-丫 dr，其中 m dtY为粘滞系数;还受到一个周围颗粒给它的涨落力R (t)，涨落力随时间t变化， mmKU可正可负，其平均
2、值为0；此外颗粒还会受一个势力-竺。根据牛顿第二定律， dr颗粒运动的朗之万方程为：d 2 rdrd Um= -y+ R (t)-dt 2 m dt mdt对于这样的二阶微分方程，通常是化为一阶微分方程组来求解，由于得到d 2 rdvdv drdv=vdt 2dtdr dtdrdrv =-dtdv = 1dt mv + R (t) 一mmdU)在这个化简的朗之万方程中，我们期望得到在某一时刻位置和速度的关系。二、物理机理1.龙格-库塔法的基本思想及一般形式设初值问题y = y (x) e a, b，由微分中值定理可知，必存在gw x , x ，使n n+1y(x ) = y(x ) + hy冗
3、)=y(x ) + hf (g, y(g)n +1nn设 y = y(x )，并记 k* = f (g, y(g)，贝ynny(x ) = y + hK *+1其中K*称为y(x)在x x 上的平均斜率，只要对平均斜率K*提供一种算法, , +1上式就给出了一种数值解公式，例如，用K = f (x ,y )代替K*，就得到欧拉公1式，用K二f (x , y )代替K*，就得到向后欧拉公式，如果用K , K的平均值2 n+1 n+11 2来代替K*,则可得到二阶精度的梯形公式。可以设想，如果在x ,x 上能多预 n n+1测几个点的斜率值，用它们的加权平均值代替K*，就有望的到具有较高精度的 数
4、值解公式，这就是龙格-库塔法的基本思想。龙格-库塔公式的一般形式：1y = y + h 工 c kn+1ni ii=11)K = f (x , y )1 n nK = f (x + 九 h, y + h艺卩 K )in i nij jj=1c ,X ,p 均为常数， i i ij其中K是y = y (x)在x +九h(0 X 1)点的斜率预测值；in ii选取这些常数的原则是使(1)式有尽可能高的精度。2、四阶龙格-库塔法四阶龙格-库塔法的表达形式如下：Iy = y + h(c K + c K + c K + c K )n+1n1 12 23 34 4K = f (x , y )1n nK =
5、 f (x + 九 h, y + hp K )2 n2 n21 1K = f (x +X h,y + hp K + hp K )3 n3n3113222)K = f (x + X h, y + hp K + hp K + hp K )4n4n411422433其中有13个待定系数，我们希望适当的选取这些系数使公式的精度尽可能 高，先将K , K ,K按照如下二元函数泰勒级数展开，取到h4项，再将K , K ,K234234展开后的值代入(2)式得到打+1的表达式。f (x +Xh,y + h工p K )=艺一(Xh一 + h工p 一“f (x ,y ) +.n i nj 1 k! i5xjn
6、nj=1k =0j=1再用泰勒公式将y(x)展开，取到h5项，得到y(x)n+1n+1由于局部截断误差Rn+1 = y(xn+1) - yn+1，泰勒级数的系数匹配，使局部截断误差为O(h5)，得到p =k , p + p =k , c kp p =212313234 2 32 4324p + p + p =九,c + c + c + c = 1414243412341c k+ c k+ c k = , c k2 + c 九2 + c k22 13 34 422 23 34 411c k + c k + c k =, c 九 p + c (九 p +九 p )=2 23 34 443 2324
7、242343611c 九九卩 + c 九(九 p +九 p ) =, c A2 p+ c (A2 p+ A2 p )=-3 2 3 324 424234383 232424234312到11个方程和13个未知量，因此需要补充两个条件九2 = 2, p31 = 0，求解得到：=1, c = 6, c1 6 21=2 p3221=3 c3=p = p =24142161=3 c4最后得出经典的四阶龙格-库塔公式：y = y + h (K + 2 K + 2 K + K )K1= f(xn,yn)n+1n 61234nhh+ ky + K )2 n 2 1hh + 百，y + K )2 n 2 2=
8、 f ( x + h, y + hK )nn3三、数值计算方案下面利用四阶龙格库塔法来求解朗之万方程，在前面我们已经推算出一阶微 分方程组：drV =-au)百)dtdv 1=(y v + R (t)dt m m m取初值为(t ,r ,v )，时间t的步长为h,得到:0 0 0hv = v + - (K + 2 K + 2 K + K )n+1 n 61234K = f (t , v )1 n n-K = f (t + , v + K )2 n 2 n 21-K = f (t + , v + K )K = f (t + - v + -K )4 nn3- r = r + - (K + 2 K + 2 K + K )n+1 n 61234K = f (t , r )1 n n-K = f (t + -,r + -K )2 n 2 n 21-K = f (t + -,r + -K )K = f (t + -, r + -K )4nn3由于已知初值为(t ,r ,v )，可以推算出在t + -时刻的r,v的大小，由此不0 0 0 0 1 1停的迭代，可以知道在任意t + n-时刻的r ,v的大小，进而可以分析颗粒运动的0n n特性。四、流程图
