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)
 #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;
}


輸出畫面
t      y(t)       w(t)      error
=====================================
0.00  1.0000000  1.0000000  0.0000000
0.10  1.0003252  1.0003252  0.0000000
0.20  1.0025385  1.0025385  0.0000000
0.30  1.0083636  1.0083636  0.0000000
0.40  1.0193599  1.0193599  0.0000000
0.50  1.0369387  1.0369387  0.0000000
0.60  1.0623767  1.0623767  0.0000000
0.70  1.0968294  1.0968294  0.0000000
0.80  1.1413421  1.1413421  0.0000000
0.90  1.1968607  1.1968607  0.0000000
1.00  1.2642411  1.2642411  0.0000000

Command exited with non-zero status 101