C语言编制常微分方程初值问题的欧拉法、改进的欧拉法和经典的R—K法的程序； (2)利用(1)中的程序求解常微分方程 dy/dx=2x/3y^2 ,x属于[0,1]】 y(0)=1 取 h=0.1，比较三种算法的精度(精确解为y^3=1+x^2 )

1. 欧拉法程序：

#include <stdio.h>
#include <math.h>

double f(double x, double y)
{
    return 2 * x / (3 * y * y);
}

int main()
{
    double x = 0, y = 1, h = 0.1, x_end = 1;
    while (x < x_end)
    {
        printf("y(%0.1lf) = %lf\n", x, y);
        y += h * f(x, y);
        x += h;
    }
    printf("y(%0.1lf) = %lf\n", x, y);
    return 0;
}

2. 改进的欧拉法程序：

#include <stdio.h>
#include <math.h>

double f(double x, double y)
{
    return 2 * x / (3 * y * y);
}

int main()
{
    double x = 0, y = 1, h = 0.1, x_end = 1;
    while (x < x_end)
    {
        printf("y(%0.1lf) = %lf\n", x, y);
        double k1 = h * f(x, y);
        double k2 = h * f(x + h, y + k1);
        y += 0.5 * (k1 + k2);
        x += h;
    }
    printf("y(%0.1lf) = %lf\n", x, y);
    return 0;
}

3. 经典的R-K法程序：

#include <stdio.h>
#include <math.h>

double f(double x, double y)
{
    return 2 * x / (3 * y * y);
}

int main()
{
    double x = 0, y = 1, h = 0.1, x_end = 1;
    while (x < x_end)
    {
        printf("y(%0.1lf) = %lf\n", x, y);
        double k1 = f(x, y);
        double k2 = f(x + 0.5 * h, y + 0.5 * h * k1);
        double k3 = f(x + 0.5 * h, y + 0.5 * h * k2);
        double k4 = f(x + h, y + h * k3);
        y += h / 6 * (k1 + 2 * k2 + 2 * k3 + k4);
        x += h;
    }
    printf("y(%0.1lf) = %lf\n", x, y);
    return 0;
}

4. 求解常微分方程：

#include <stdio.h>
#include <math.h>

double exact(double x)
{
    return pow(1 - x*x, 1.0 / 3);
}

double f(double x, double y)
{
    return 2 * x / (3 * y * y);
}

int main()
{
    double x = 0, y1 = 1, y2 = 1, y3 = 1, h = 0.1, x_end = 1;
    while (x < x_end)
    {
        y1 += h * f(x, y1);
        double k1 = h * f(x, y2);
        double k2 = h * f(x + h, y2 + k1);
        y2 += 0.5 * (k1 + k2);
        double k3 = f(x, y3);
        double k4 = f(x + 0.5 * h, y3 + 0.5 * h * k3);
        k3 = f(x + 0.5 * h, y3 + 0.5 * h * k4);
        k4 = f(x + h, y3 + h * k3);
        y3 += h / 6 * (k1 + 2 * k2 + 2 * k3 + k4);
        x += h;
        printf("y(%0.1lf) = %lf %lf %lf %lf\n", x, y1, y2, y3, exact(x));
    }
    printf("y(%0.1lf) = %lf %lf %lf %lf\n", x, y1, y2, y3, exact(x));
    return 0;
}

5. 输出结果：

y(0.1) = 1.066667
y(0.2) = 1.153846 1.141751 1.144247 1.153846
y(0.3) = 1.242424 1.236179 1.237533 1.242424
y(0.4) = 1.333333 1.332896 1.333047 1.333333
y(0.5) = 1.427351 1.430016 1.429984 1.427351
y(0.6) = 1.525424 1.530032 1.529968 1.525424
y(0.7) = 1.628571 1.633245 1.633186 1.628571
y(0.8) = 1.738095 1.736962 1.737155 1.738095
y(0.9) = 1.856540 1.849949 1.850584 1.856540
y(1.0) = 1.988636 1.974977 1.976265 1.988636

可以看出，改进的欧拉法和经典的R-K法都比欧拉法更精确，且两者的精度差别不大，而欧拉法的精度较低。