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

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

/***********************************************
****EULER METHOD FOR DIFFERENTIAL EQUATIONS*****
* C語言 例題5-3 利用 Modofied Euler Method
* 求ODE y'=y+2 e^(4t) y(0)=-3 , 0<= t <=1
* 取 h=0.01 
*
***********************************************/
#include<stdio.h>
#include<math.h>
#define W(x)  ( (2.0/3.0)*exp(4.0*x) - (11.0/3)*exp(x))

/*Define the RHS of the first order differential equation here(Ex: dy/dx=f(x,y))  */
double f(double x, double y){
return y+2*exp(4*x);
}

double predict(double x, double y, double h)
{
    // value of next y(predicted) is returned
    double y1p = y + h * f(x, y);
    return y1p;
}


// using Modified Euler method
double correct(double x, double y, double x1, double y1,double h)
{
    // (x, y) are of previous step
    // and x1 is the increased x for next step
    // and y1 is predicted y for next step
    double e = 0.00001;
    double y1c = y1;

    do {
        y1 = y1c;
        y1c = y + 0.5 * h * (f(x, y) + f(x1, y1));
    } while (fabs(y1c - y1) > e);

    // every iteration is correcting the value
    // of y using average slope
    return y1c;
}



int main()
{
    // here x and y are the initial
    // given condition, so x=0 and y=15
    double x = 0, y = -3;

    // final value of x for which y is needed
    double xn = 1;

    // step size
    double h = 0.01;
    int i=1;
        while (x < xn) {
            double x1 = x + h;
            double y1p = predict(x, y, h);
            double y1c = correct(x, y, x1, y1p, h);
            if (i%10==0)
                printf("Euler Modified method---at x= %2.2lf y=%2.6lf ,真值解=%2.6lf , 差值=%2.6lf \n" ,x1 ,y1c,W(x1) ,fabs(W(x1)-y1c));
            x = x1;
            y = y1c;
            i++;
    }

return 0;
}


輸出畫面
Euler Modified method---at x= 0.10 y=-3.057701 ,真值解=-3.057744 , 差值=0.000042
Euler Modified method---at x= 0.20 y=-2.994672 ,真值解=-2.994783 , 差值=0.000111
Euler Modified method---at x= 0.30 y=-2.735851 ,真值解=-2.736071 , 差值=0.000220
Euler Modified method---at x= 0.40 y=-2.167611 ,真值解=-2.168002 , 差值=0.000391
Euler Modified method---at x= 0.50 y=-1.118620 ,真值解=-1.119274 , 差值=0.000654
Euler Modified method---at x= 0.60 y=0.668738 ,真值解=0.667682 , 差值=0.001056
Euler Modified method---at x= 0.70 y=3.581004 ,真值解=3.579338 , 差值=0.001666
Euler Modified method---at x= 0.80 y=8.197291 ,真值解=8.194703 , 差值=0.002588
Euler Modified method---at x= 0.90 y=15.384255 ,真值解=15.380278 , 差值=0.003977
Euler Modified method---at x= 1.00 y=26.437797 ,真值解=26.431733 , 差值=0.006064