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

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


// C code for solving the differential equation
// using Predictor-Corrector or Modified-Euler method
// with the given conditions, y(0) = 0.5, step size(h) = 0.2
// to find y(1)
 
#include <stdio.h>
#include <math.h>

// consider the differential equation
// for a given x and y, return v
double f(double x, double y)
{
    double v = y + 2 * exp(4*x) ;
    return v;
}
 
// predicts the next value for a given (x, y)
// and step size h using Euler method
double predict(double x, double y, double h)
{
    // value of next y(predicted) is returned
    double y1p = y + h * f(x, y);
    return y1p;
}
 
// corrects the predicted value
// 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;
}
 
void printFinalValues(double x, double xn, double y, double h)
{
    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(" predict----The value of y , at x = %2.4lf is :%2.4lf\n " ,x1 ,y1p);
                printf(" correct----The value of y , at x = %2.4lf is :%2.4lf\n " ,x1 ,y1c);
            }
        x = x1;
        y = y1c;
        i++;
    }
 
    // at every iteration first the value
    // of for next step is first predicted
    // and then corrected.
    printf("\n\nThe final value of y  , at x = %2.4lf is :%2.4lf " ,x ,y);
}
 
int main()
{
    // here x and y are the initial
    // given condition, so x=0 and y=0.5
    double x = 0, y = -3;
 
    // final value of x for which y is needed
    double xn = 1;
 
    // step size
    double h = 0.01; 
    printFinalValues(x, xn, y, h);
 
    return 0;
}


   輸出畫面
  predict----The value of y , at x = 0.1000 is :-3.0583 ---修正前
  correct----The value of y , at x = 0.1000 is :-3.0577 ---修正後

  predict----The value of y , at x = 0.2000 is :-2.9956
  correct----The value of y , at x = 0.2000 is :-2.9947

  predict----The value of y , at x = 0.3000 is :-2.7373
  correct----The value of y , at x = 0.3000 is :-2.7359

  predict----The value of y , at x = 0.4000 is :-2.1699
  correct----The value of y , at x = 0.4000 is :-2.1676

  predict----The value of y , at x = 0.5000 is :-1.1222
  correct----The value of y , at x = 0.5000 is :-1.1186

  predict----The value of y , at x = 0.6000 is :0.6633
  correct----The value of y , at x = 0.6000 is :0.6687

  predict----The value of y , at x = 0.7000 is :3.5728
  correct----The value of y , at x = 0.7000 is :3.5810

  predict----The value of y , at x = 0.8000 is :8.1849
  correct----The value of y , at x = 0.8000 is :8.1973

  predict----The value of y , at x = 0.9000 is :15.3656
  correct----The value of y , at x = 0.9000 is :15.3843

  predict----The value of y , at x = 1.0000 is :26.4097
  correct----The value of y , at x = 1.0000 is :26.4378


  The final value of y  , at x = 1.0000 is :26.4378


The final value of y  , at x = 1.0000 is :26.4378