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>

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

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

2. 改进的欧拉法程序：

#include <stdio.h>

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

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

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

#include <stdio.h>

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

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

4. 求解 dy/dx=2x/3y^2 ,x属于[0,1] y(0)=1，比较三种算法的精度：

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

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

int main() {
    double x0 = 0, y0 = 1, h = 0.1, xn = 1;
    double x = x0, y = y0;
    while (x < xn) {
        y += h * f(x, y);
        x += h;
    }
    double y_euler = y;
    printf("Euler method: y(%lf) = %lf, error = %lf\n", xn, y_euler, pow(y_euler * y_euler - 1, 0.5));
    
    x = x0, y = y0;
    while (x < xn) {
        double k1 = f(x, y);
        double k2 = f(x + h, y + h * k1);
        y += h * (k1 + k2) / 2;
        x += h;
    }
    double y_improved_euler = y;
    printf("Improved Euler method: y(%lf) = %lf, error = %lf\n", xn, y_improved_euler, pow(y_improved_euler * y_improved_euler - 1, 0.5));
    
    x = x0, y = y0;
    while (x < xn) {
        double k1 = f(x, y);
        double k2 = f(x + h / 2, y + h * k1 / 2);
        double k3 = f(x + h / 2, y + h * k2 / 2);
        double k4 = f(x + h, y + h * k3);
        y += h * (k1 + 2 * k2 + 2 * k3 + k4) / 6;
        x += h;
    }
    double y_runge_kutta = y;
    printf("Runge-Kutta method: y(%lf) = %lf, error = %lf\n", xn, y_runge_kutta, pow(y_runge_kutta * y_runge_kutta - 1, 0.5));
    
    return 0;
}