龙格库塔的计算写法

参考《导弹飞行力学》

对部分参数的解释:

dx/dt=f(t,x):   之前一直看不懂f（t,x）到底指的哪个式子，其实在开头就提出来了，f是dy 

 K2=△t*f(tk+△t/2,xk+1/2*K1)：t在导弹飞行力学，是y(0)(M中是1)，所以可忽略。然后xk+1/2*K1即dx/dt的变量进行变换，在飞行力学没有直接加，需要通过y来改变dy

而y的变换则是借助y+h*dy去改变，进而代入右函数改变dy.

龙格库塔公式：

xk+1=xk+(K1+2K2+2K3+K4)/6         dx=f(t,x)               对应的  应该是   dy=dery(y)           y自带t，所以不用管原式t.

K1=h*f(tk,xk)                          K1=h*dy                           K1=h*d(y_old)            

K2=h*f(tk+h/2,xk+K1/2)       K2=h*dy/2                         dy2=dery(y_old+K1/2)        K2=h*dy2

K3=h*f(tk+h/2,xk+K2/2)       K3=h*dy/2                         dy3=dery(y_old+K2/2)       K3=h*dy3

K4=h*f(tk+h,xk+K3)             K4=h*dy/2                          dy4=dery(y_old+K3)          k4=h*dy4          y_new=y_old+1/6(K1+2K2+2K3+K4);

所以主要的步骤

C语言版本  和导弹飞行力学后面的一模一样

1）迭代

void rk(int n, double h)                                 //传递量为数组数量和步长
{
int i, j;
double a[4], old_y[dim], y1[dim], *dy;                  

dy = (double *)malloc(n*sizeof(double));             //赋内存

a[0] = h / 2;
a[1] = a[0];
a[2] = h;
a[3] = h;                    //可以单纯记一下，0,1h,23,二分之一h

grkt2f(dy, y);              //先运行一次

for (i = 0; i<n; i++)
{
old_y[i] = y[i];             //留旧
}

for (j = 0; j<3; j++)                                                   //迭代3次，之后补上
{
for (i = 0; i<n; i++)                              
{
y1[i] = old_y[i] + a[j] * dy[i];              //涉及到迭代，正好是1/2 1/2 1
y[i] = y[i] + a[j + 1] * dy[i] / 3.;            //此句不涉及到迭代，只是相加的式子.  1/2 1 1刚好
}

grkt2f(dy, y1);                 //更新
}

for (i = 0; i<n; i++)
{
y[i] = y[i] + h*dy[i] / 6.;                      //补上最后一次  1/2 即1/6
}

free(dy);                   //释放内存，在C/C+语言很重要
return;
}

2）终止

do{}while(y(r)>=0)等等

Matlab 版本

问题：如何针对步长做出调整（大步长可能会在末尾出错，小步长计算时间太长）

1）迭代：可省n,内存分配，n的内循环

function yl=rk(h,y)
dy=dif_missle_trace1(y);
old_y=y;
K1=h*dy;
y=y+K1/2;
dy=dif_missle_trace1(y);
K2=h*dy;
y=y+K2/2;
dy=dif_missle_trace1(y);
K3=h*dy;
y=y+K3;
dy=dif_missle_trace1(y);
K4=h*dy;
yl=old_y+(K1+2*K2+2*K3+K4)/6;

2）终止：没有do while ，可以先做一次，然后补上。