範例 EX5-3 使用 Euler修正法 ,求ODE y'=y+2 e^(4t) y(0)=-3 , 0<= t <=1 ---C語言

C語言 例題5-3 使用 Euler修正法　,求ODE y'=y+2 e^(4t) y(0)=-3 , 0<= t <=1 y(1)= ??

double h=0.01
n=(b-a)/h;

結果：　h=0.01　比　　h=0.0001　誤差還小
===========================================
Euler修正法　h=0.01
　t        y(t)        　  w(t)    　　   error
1.00   26.433458   26.431733    0.001725
===========================================
Euler法　h=0.0001
　t                y(t)           　  w(t)    　　         error
1.00   26.423010   26.431733    0.008724  


/* ex5-3.c based on modified Euler
 * 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+2*exp(4*t))
#define W(t)   ((2.0/3.0)*exp(4.0*t)-(11.0/3)*exp(t))
 void main()
 {
    int i,n;
    double h=0.01,a=0.0,b=1.0,y0=-3.0,k1,k2,y,t;
    n=(b-a)/h;
    t=a;
    y=y0;
    printf("t        y(t)          w(t)         error\n");
    printf("=========================================\n");
    printf("%.2lf  %10.6lf  %10.6lf  %10.6lf\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+(k1+k2)/2.0;
       t=a+i*h;
       if(i%10==0)
          printf("%.2lf  %10.6lf  %10.6lf  %10.6lf\n",
         t,y,W(t),fabs(y-W(t)));
    }
    return;
}


＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
輸出畫面
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

t        y(t)          w(t)         error

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

0.00   -3.000000   -3.000000    0.000000

0.10   -3.057725   -3.057744    0.000018

0.20   -2.994738   -2.994783    0.000045

0.30   -2.735987   -2.736071    0.000084

0.40   -2.167862   -2.168002    0.000140

0.50   -1.119051   -1.119274    0.000223

0.60    0.668025    0.667682    0.000343

0.70    3.579857    3.579338    0.000519

0.80    8.195482    8.194703    0.000779

0.90   15.381439   15.380278    0.001161

1.00   26.433458   26.431733    0.001725