四阶龙格库塔法-常微分方程数值计算

龙格-库塔法是数值求解常微分方程经典方法，
下面是百度百科的四阶龙格库塔法的公式截图

h是t的步长
设n为0时，t0=0.0，y0也已知，可以计算k1，k2，k3，k4，然后y1=y0+h*(k1+2.0* k2+2.0 *k3+k4)/6.0，t1=t0+h
这样由t0算出了t1，y1，
然后再带公式推算出t2，y2
写成程序就是

double f(t, y)
{
    2.0*t;
}
double step_rk4(double tn，double yn, double h)
{
   double h2=0.5*h;
   double h6=h/6.0;
   double k1=f(tn,fn);
   double k2=f(tn+h2,yn+h2*k1);
   double k3=f(tn+h2,yn+h2*k2);
   double k4=f(tn+h,yn+h*k3);
   double y_nadd1=yn+h6*(k1+2.0* k2+2.0 *k3+k4);
   return y_nadd1;
}

这样可以从t0开始一步步推算

#include<iostream>

double f(t, y)
{
    2.0*t;
}
double step_rk4(double tn，double yn, double h)
{
   double h2=0.5*h;
   double h6=h/6.0;
   double k1=f(tn,fn);
   double k2=f(tn+h2,yn+h2*k1);
   double k3=f(tn+h2,yn+h2*k2);
   double k4=f(tn+h,yn+h*k3);
   double y_nadd1=yn+h6*(k1+2.0* k2+2.0 *k3+k4);
   return y_nadd1;
}

int main()
{
    double t=0.0;
    double y=0.0;
    double y_next=0.0;
    double h=0.1;
    
    for(int i=0;i<1000;i++)
    {
       std::cout << "(" << t << ", " << y << ")\n";
       
       y_next=step_rk4(t,y,h);
       t+=h;
       y=y_next;
    }
    return 0;
}