C語言 例題5-2用修正Euler尤拉近似法 解一階常微分方程ODE y' = -y+x+1 0<= x <= 1 , y(0)=1 , h=0.2

C語言 例題5-2用修正Euler尤拉近似法 解一階常微分方程ODE y' = -y+x+1 0<= x <= 1 , y(0)=1 , h=0.2

源自於 https://www.geeksforgeeks.org/predictor-corrector-or-modified-euler-method-for-solving-differential-equation/

For a given differential equation  with initial condition 
find the approximate solution using Predictor-Corrector method.

Predictor-Corrector Method :
The predictor-corrector method is also known as Modified-Euler method.
In the Euler method, the tangent is drawn at a point and slope is calculated for a given step size. Thus this method works best with linear functions, but for other cases, there remains a truncation error. To solve this problem the Modified Euler method is introduced. In this method instead of a point, the arithmetic average of the slope over an interval  is used.

Thus in the Predictor-Corrector method for each step the predicted value of is calculated first using Euler’s method and then the slopes at the points  and  is calculated and the arithmetic average of these slopes are added to  to calculate the corrected value of .

/* C Program to find approximation of a ordinary
   differential equation using modified Euler's method.
 
    This method is similar to Euler's method ,
    except here we use the average of slopes instead of one.

    x is calculate like the earlier method  and
    y1 = y0 + h*f(x0,y0)
    y2 = y1 +  h*(f(x0,y0)+f(x1,y1))/2
    y3 = y2 +  h*(f(x1,y1)+f(x2,y2))/2
    .
    .
    .
    .
    yn+1 = yn +  h*(f(xn,yn)+f(xn+1,yn+1))/2
    and so on...
*/

#include <stdio.h>
#include <math.h>
#define  w(t)    (exp(-t)+t)
// consider the differential equation
// for a given x and y
float f(float(x),float(y))
{
    return (-y+x+1);
}

// 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)
{
 
    printf("\n\tx\t\t\ty1p\t\t\ty1c\t\t\treal\t\terr\n");
    printf("=================================================================\n");
    while (x < xn) {
        double x1 = x + h;
        double y1p = predict(x, y, h);
        double y1c = correct(x, y, x1, y1p, h);
       
        printf("\t%0.2lf\t\t%0.4lf\t\t%0.4lf\t\t%0.4lf\t\t%0.4lf\n",x1,y1p,y1c,w(x1),fabs(y1c-w(x1)) );
        x = x1;
        y = y1c;
    }
 
    // at every iteration first the value
    // of for next step is first predicted
    // and then corrected.
    printf("\n\tThe final value of y at x =%0.2lf is : %0.4lf" ,x , y);
}
 
int main()
{
    // here x and y are the initial
    // given condition, so x=0 and y=1
    double x = 0, y = 1.0;
    printf("The value of x0=%0.2lf , ",x);
    printf("The Value of y0=%0.2lf\n",y);
    // final value of x for which y is needed
    double xn = 1;
 
    // step size
    double h = 0.2;
 
    printFinalValues(x, xn, y, h);
 
    return 0;
}


The value of x0=0.00 , The Value of y0=1.00

x y1p y1c real err
=========================================================
0.20 1.0000 1.0182 1.0187 0.0005
0.40 1.0545 1.0694 1.0703 0.0009
0.60 1.1355 1.1477 1.1488 0.0011
0.80 1.2382 1.2481 1.2493 0.0012
1.00 1.3585 1.3666 1.3679 0.0012

The final value of y at x =1.00 is : 1.3666