C語言 例題5-2 利用 Forward Euler Method求ODE y'=y+2 e^(4t) y(0)=-3 , 0<= t <=1 取 h=0.01 , 0.001 , 0.0001

C語言 例題5-2 利用 Forward Euler Method求ODE y'=y+2 e^(4t)  y(0)=-3 , 0<= t <=1 取 h=0.01 , 0.001 , 0.0001 真實值 = #define  w(t)   ((2.0/3.0)*exp(4.0*t)-(11.0/3)*exp(t))

/* ex5-1.c Forward Euler Method is used for
 * solving y'=f(y,t) of first order
 * ordinary differential equation
 * with initial condition y(0)=y0
 * known.

 #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))
要加小數點 否則當成整數 2與2.0 不同
 */
#include <stdio.h>
#include <math.h>
//===========================
double w(double t)
{
    double w1,w2;
    w1=(2.0/3.0)*exp(4.0*t);
    w2=(11.0/3)*exp(t);
    return w1-w2;
}
//===========================
double  F(double y, double t)
{
    return (y+2*exp(4*t));
}   
//===========================
void main()
{
    double h[]={0.01,0.001,0.0001},hj;
    int i,j,n;
    double y,t,y0=-3,t0=0,t1=1;
   
    for(j=0;j<=2;j++)
    {
        hj=h[j];
        printf("The h=%.4lf\n\n",hj);
        n=(t1-t0)/hj ;

        y=y0;
        t=t0;
        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++)
        {
            //printf("j=%1d",j);
            y=y+hj*F(y,t);
            t=t+hj;
            switch (j)
            {
                case 0:
                     if(i%10==0) 
                        printf("%.2lf  %10.6lf  %10.6lf  %10.6lf\n",t,y,w(t),fabs(y-w(t)));
                case 1:
                     if(i%100==0) 
                        printf("%.2lf  %10.6lf  %10.6lf  %10.6lf\n",t,y,w(t),fabs(y-w(t)));
                case 2:
                     if(i%1000==0) 
                        printf("%.2lf  %10.6lf  %10.6lf  %10.6lf\n",t,y,w(t),fabs(y-w(t)));
            }
           
       
        }
        printf("\n\n");
    }   
    return;
}

輸出畫面
The h=0.0100

t         y(t)         w(t)        error
========================================
0.00   -3.000000   -3.000000    0.000000
0.10   -3.062524   -3.057744    0.004780
0.20   -3.007973   -2.994783    0.013190
0.30   -2.763301   -2.736071    0.027230
0.40   -2.217918   -2.168002    0.049915
0.50   -1.205054   -1.119274    0.085780
0.60    0.526053    0.667682    0.141629
0.70    3.351680    3.579338    0.227657
0.80    7.835576    8.194703    0.359127
0.90   14.821414   15.380278    0.558865
1.00   25.570747   26.431733    0.860987
1.00   25.570747   26.431733    0.860987


The h=0.0010

t         y(t)         w(t)        error
========================================
0.00   -3.000000   -3.000000    0.000000
0.10   -3.058228   -3.057744    0.000484
0.20   -2.996118   -2.994783    0.001335
0.30   -2.738828   -2.736071    0.002757
0.40   -2.173055   -2.168002    0.005052
0.50   -1.127956   -1.119274    0.008682
0.60    0.653348    0.667682    0.014334
0.70    3.556298    3.579338    0.023040
0.80    8.158359    8.194703    0.036344
0.90   15.323721   15.380278    0.056557
1.00   26.344602   26.431733    0.087131
1.00   26.344602   26.431733    0.087131


The h=0.0001

t         y(t)         w(t)        error
========================================
0.00   -3.000000   -3.000000    0.000000
0.10   -3.057792   -3.057744    0.000048
0.20   -2.994917   -2.994783    0.000134
0.30   -2.736347   -2.736071    0.000276
0.40   -2.168508   -2.168002    0.000506
0.50   -1.120143   -1.119274    0.000869
0.60    0.666247    0.667682    0.001435
0.70    3.577031    3.579338    0.002307
0.80    8.191065    8.194703    0.003639
0.90   15.374616   15.380278    0.005663
1.00   26.423010   26.431733    0.008724