龙格-库塔(Runge-Kutta)方法数学原理及实现

我就是我 2022-08-04 16:42 145阅读 0赞

# 龙格-库塔(Runge-Kutta)方法 #

龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法。由于此算法精度高，采取措施对误差进行抑制，所以其实现原理也较复杂。该算法是构建在数学支持的基础之上的。

--------------------

对于一阶精度的欧拉公式有：

yi\+1=yi\+hki

其中  h 为步长，则  yi\+1 的表达式与  y(xi\+1) 的Taylor展开式的前两项完全相同，即 **局部截断误差**为  O(h2) 。  
当用点  xi 处的斜率近似值  k1 与右端点  xi\+1 处的斜率  k2 的算术平均值作为平均斜率  k∗ 的近似值，那么就会得到二阶精度的改进欧拉公式：

yi\+1=yi\+h(k1\+k2)

其中  k1=f(xi,yi) ，  k2=f(xi\+h,yi\+hk1)  
依次类推，如果在区间  \[xi,xi\+1\] 内多预估几个店上的斜率值  k1,k2,…,km ，并用他们的加权平均数作为平均斜率  k∗ 的近似值，显然能够构造出具有很高精度的高阶计数公式。  
上述两组公式在形式删过的共同点：都是用  f(x,y) 在某些点上值得线性组合得出  y(xi\+1) 的近似值  yi\+1 ，且增加计算的次数，可以提高截断误差的阶，他们的误差估计可以用  f(x,y) 在  xi 处的Taylor展开来估计。

--------------------

于是可考虑用函数f(x,y)在若干点上的函数值的线性组合老构造金斯公式，构造时要求近似公式在f(xi,yi)处的Taylor展开式与解y(x)在xi处的Taylor展开式的前面几项重合，从而使金斯公式达到所需要的阶数。既避免求高阶导数，又提高了计算方法精度的阶数。或者说，在\[xi,xi\+1\]这一步内计算多个点的斜率值，若够将其进行加权平均作为平均斜率，则可构造出更高精度的计算格式，这就是龙格-库塔(Runge-Kutta)方法。  
一般的龙格-库塔法的形式为

![这里写图片描述][20150522223232202]  
称为P阶龙格-库塔方法。  
其中  ai,bij,cj 为待定参数，要求上式  yi\+1 在点  (xi,yi) 处作Taylor展开，通过相同项的系数确定参数。

--------------------

当然，经典的龙格-库塔方法是四阶的。也就是在\[xi,xi\+1\]上用四个点处的斜率加权平均作为平均斜率k∗的近似值，构成一系列四阶龙格-库塔公式。具有四阶精度，即局部截断误差是O(h5)。  
下面介绍最常用的一种四阶龙格-库塔方法。  
设

yi\+1=yi\+c1K1\+c2K2\+c3K3\+c4K4

这里  K1,K2,K3,K4 为四个不同点上的函数值，分别设其为

![这里写图片描述][20150522224118450]  
其中  c1,c2,c3,c4,a2,a3,a4,b21,b31,b32,b41,b42,b43 均为待定系数。  
把  K2,K3,K4 分别在  xi 点占城h的幂级数，带入线性组合式中，将得到的公式与  y(xi\+1) 在  xi 点上的泰勒展开式比较，使其两式右端知道  h4 的系数相等，经过较复杂的解方程过程便可得到关于  ai,bij,cj 的一组特解。

a2=a3=b21=b32=12

b31=b41=b42=0

a4=b43=1

c1=c4=16

c2=c3=13

带入之后得到  
![这里写图片描述][20150522225041443]

龙格-库塔方法的推导基于Taylor展开方法，因而它要求所求的解具有较好的光滑性。如果解的光滑性差，那么，使用四阶龙格-库塔方法求得的数值解，其精度可能反而不如改进的欧拉方法。在实际计算时，应正对问题的具体特点选择适合的算法。对于光滑性不太好的解，最好采用**低阶算法**而将步长h<script type="math/tex" id="MathJax-Element-47">h</script>取小。

--------------------

龙格-库塔法的C语言实现

#include "stdio.h"
    #include "stdlib.h"
    
    void RKT(t,y,n,h,k,z)
    int n;              /*微分方程组中方程的个数，也是未知函数的个数*/
    int k;              /*积分的步数(包括起始点这一步)*/
    double t;           /*积分的起始点t0*/
    double h;           /*积分的步长*/
    double y[];         /*存放n个未知函数在起始点t处的函数值,返回时,其初值在二维数组z的第零列中*/
    double z[];         /*二维数组,体积为n x k.返回k个积分点上的n个未知函数值*/
    {
        extern void Func();             /*声明要求解的微分方程组*/
        int i,j,l;
        double a[4],*b,*d;
        b=malloc(n*sizeof(double));     /*分配存储空间*/
        if(b == NULL)
        {
            printf("内存分配失败\n");
            exit(1);
        }
        d=malloc(n*sizeof(double));     /*分配存储空间*/
        if(d == NULL)
        {
            printf("内存分配失败\n");
            exit(1);
        }
        /*后面应用RK4公式中用到的系数*/
        a[0]=h/2.0;                     
        a[1]=h/2.0;
        a[2]=h; 
        a[3]=h;
        for(i=0; i<=n-1; i++) 
            z[i*k]=y[i];                /*将初值赋给数组z的相应位置*/
        for(l=1; l<=k-1; l++)
        {
            Func(y,d);
            for (i=0; i<=n-1; i++)
                b[i]=y[i];
            for (j=0; j<=2; j++)
            {
                for (i=0; i<=n-1; i++)
                {
                    y[i]=z[i*k+l-1]+a[j]*d[i];
                    b[i]=b[i]+a[j+1]*d[i]/3.0;
                }
                Func(y,d);
            }
            for(i=0; i<=n-1; i++)
              y[i]=b[i]+h*d[i]/6.0;
            for(i=0; i<=n-1; i++)
              z[i*k+l]=y[i];
            t=t+h;
        }
        free(b);            /*释放存储空间*/
        free(d);            /*释放存储空间*/
        return;
    }
    main()
    {
        int i,j;
        double t,h,y[3],z[3][11];
        y[0]=-1.0; 
        y[1]=0.0; 
        y[2]=1.0;
        t=0.0; 
        h=0.01;
        RKT(t,y,3,h,11,z);
        printf("\n");
        for (i=0; i<=10; i++)           /*打印输出结果*/
        {
            t=i*h;
            printf("t=%5.2f\t ",t);
            for (j=0; j<=2; j++)
              printf("y(%d)=%e ",j,z[j][i]);
            printf("\n");
        }
    }
    
    void Func(y,d)
    double y[],d[];
    {
        d[0]=y[1];      /*y0'=y1*/
        d[1]=-y[0];     /*y1'=y0*/
        d[2]=-y[2];     /*y2'=y2*/
        return;
    }

ps:如果有时间的话，可能会回过头来加一分解方程的推到吧…

[20150522223232202]: /images/20220805/1b102ca8f74945299aa128d5751c96f1.png
[20150522224118450]: /images/20220805/0c832cee9ae14b90b050749a9146702d.png
[20150522225041443]: /images/20220805/c06027152c15410e95854b5172b49cb7.png