範例EX5-5.c 使用 二階Runge-Kutta 解 ODE y'= -y + t^2 + 1 , 0<=t<=1 , y(0)=1 , 真實解 W(t)= -2e^(-t) + t ^2 - 2t + 3

範例EX5-5.c 使用 二階Runge-Kutta 解 ODE y'= -y + t^2 + 1 , 0<=t<=1 , y(0)=1 , 真實解　W(t)= -2e^(-t) + t ^2 - 2t + 3

/* ex5-5.c Second Order Runge-Kutta Method is used
 * for solving Ordinary Differential Equation of
 * y'=f(y,t) with initial condition of y(t0)=y0.
 */
#include <stdio.h>
#include <math.h>
#define F(y,t)   (-y+t*t+1)
#define W(t)     (-2*(1.0/exp(t))+pow(t,2)-2*t+3)
void main()
{
   int i,n=100;
   double h,a=0.0,b=1.0,t0,t,y0=1.0,y,k1,k2;
   h=(b-a)/n;
   y=y0;
   t0=a;
   t=t0;
   printf("t      y(t)       w(t)       error\n");
   printf("=====================================\n");
   printf("%.2lf %10.7lf %10.7lf %10.7lf\n",
  t,y,W(t),fabs(y-W(t)));
   for(i=1;i<=n;i++)
   {
      k1=h*F(y,t);
      k2=h*F((y+k1),(t+h));
      y=y+0.5*(k1+k2);
      t=t+h;
      if(i%10==0)
         printf("%.2lf %10.7lf %10.7lf %10.7lf\n",
         t,y,W(t),fabs(y-W(t)));
   }
   return;
}

t      y(t)       w(t)       error                                                                                                                            

=====================================                                                                                                                         

0.00  1.0000000  1.0000000  0.0000000                                                                                                                         

0.10  1.0003269  1.0003252  0.0000017                                                                                                                         

0.20  1.0025421  1.0025385  0.0000036                                                                                                                         

0.30  1.0083691  1.0083636  0.0000056                                                                                                                         

0.40  1.0193675  1.0193599  0.0000076                                                                                                                         

0.50  1.0369483  1.0369387  0.0000096                                                                                                                         

0.60  1.0623883  1.0623767  0.0000116                                                                                                                         

0.70  1.0968430  1.0968294  0.0000136                                                                                                                         

0.80  1.1413577  1.1413421  0.0000156                                                                                                                         

0.90  1.1968782  1.1968607  0.0000175                                                                                                                         

1.00  1.2642605  1.2642411  0.0000194                                                                                                                         

                                                                                                                                                              

                                                                                                                                                              

...Program finished with exit code 38                                                                                                                         

Press ENTER to exit console.