C語言 例題5-6 四階 Runge-Kutta 解 ODE y'= -y + t^2 + 1 , 0<=t<=1 , y(0)=1 , 真實解 W(t)= -2e^(-t) + t ^2 - 2t + 3

C語言 例題5-6  四階 Runge-Kutta 解 ODE y'= -y + t^2 + 1 , 0<=t<=1 , y(0)=1 , 真實解　W(t)= -2e^(-t) + t ^2 - 2t + 3

/* ex5-6.c based on Four-Order Runge-Kutta
 * Method to approximate the solution of the
 * initial-value problem
 *   y'=f(y,t), a<=t<=b, y(a)=y0
 * at (n+1) equally spaced numbers in the interval
 * [a,b]: input a,b,n,and initial condition y0.
 */
 #include <stdio.h>
 #include <math.h>
 #define  F(y,t)   (-y+t*t+1)   // dy/dt= -y + t^1 +1
 #define  W(t)     (-2*(1/exp(t))+pow(t,2)-2*t+3)
 void main()
 {
    int i,n=100;
    double a=0.0,b=1.0,y0=1.0,k1,k2,k3,k4,h,t,y,err;
    h=(b-a)/n;
    t=a;
    y=y0;
    err=fabs(y-W(t));
    printf("t      y(t)       w(t)      error\n");
    printf("=====================================\n");
    printf("%.2lf %10.7lf %10.7lf %10.7lf\n",t,y,W(t),err);
    for(i=1;i<=n;i++)
    {
       k1=h*F(y,t);
       k2=h*F((y+k1/2.0),(t+h/2.0));
       k3=h*F((y+k2/2.0),(t+h/2.0));
       k4=h*F((y+k3),(t+h));
       y=y+(k1+2*k2+2*k3+k4)/6.0;
       t=a+i*h;
       err=fabs(y-W(t));
       if(i%10==0)
          printf("%.2lf %10.7lf %10.7lf %10.7lf\n",t,y,W(t),err);
    }
    return;
}