C/C++语言实现改进欧拉格式

这里使用改进欧拉格式的平均化形式：

y p = y n + h f ( x n , y n ) y c = y n + h f ( x n + 1 , y p ) y n + 1 = 1 2 ( y p + y c ) y_p = y_n + h f(x_n,y_n) \ y_c = y_n + h f(x_{n+1},y_p) \ y_{n+1} = frac{1}{2}(y_p + y_c) yp​=yn​+hf(xn​,yn​)yc​=yn​+hf(xn+1​,yp​)yn+1​=21​(yp​+yc​)

得到的欧拉格式如下：

// 改进欧拉格式 
void ModifiedEuler(double x0,double y0,double h,int N)
{
	
    double yp = 0,yc = 0,y1 = 0,x1 = 0 ; // x1,y1为新值,即x_n+1,y_n+1 
    // 逐步计算结果 
    for(int i = 0;i < N; i++,x0=x1,y0=y1)
    {
        x1 = x0 + h;
        yp = y0 + h * f(x0,y0);
        yc = y0 + h * f(x1,yp);
        y1 = 0.5 * (yp + yc);
        printf("%lf  %lfn",x1,y1); 
    }
}

给定初值问题：
{ d y d x = y x + x e 2 , 1 < x ≤ 2 y ( 1 ) = 0 , begin{cases} frac{dy}{dx} = frac{y}{x}+xe^2,1dxdy​=xy​+xe2,1 其精确解为 y = x ( e x − e ) y =x(e^x-e) y=x(ex−e)，步长取0.1

#include 
#include 

// 右函数f(x,y)
double f(double x,double y)
{
    return (y/x + x * exp(2));
}

// 精确解 
double exact(double x)
{
    return (x  *(exp(x) - exp(1) ));
}

// 改进欧拉格式 
void ModifiedEuler(double x0,double y0,double h,int N)
{
	//  x0、y0为老值,h 为步长, N 为步数
    double yp = 0,yc = 0,y1 = 0,x1 = 0 ; // x1,y1为新值 
    printf(" xn%8syn%4sy(xn)%4s误差n"," "," "," ");
    // 逐步计算结果 
    for(int i = 0;i < N; i++,x0=x1,y0=y1)
    {
        x1 = x0 + h;
        yp = y0 + h * f(x0,y0);
        yc = y0 + h * f(x1,yp);
        y1 = 0.5 * (yp + yc);
        printf(" %.2lf %8.4lf %8.4lf %8.4lfn",x1,y1,exact(x1),fabs(exact(x1)-y1)); 
    }
}

int main()
{
    printf("步长为0.1时,改进欧拉格式:n");
    ModifiedEuler(1,0,0.1,10);
    return 0;
}

结果如下：