利用欧拉法，改进的欧拉法，和经典R-K法中的程序求解常微分方程 dy/dx=2x/3y^2 ,x属于[0,1] y(0)=1 取 h=0.1，比较三种算法的精度(精确解为y^3=1+x^2 )

欧拉法的程序：

import math

def euler(dydx, x0, y0, h, x_end):
    x = [x0]
    y = [y0]
    while x[-1] < x_end:
        y_next = y[-1] + dydx(x[-1], y[-1]) * h
        y.append(y_next)
        x.append(x[-1] + h)
    return x, y

def dydx(x, y):
    return 2 * x / (3 * y ** 2)

x0 = 0
y0 = 1
h = 0.1
x_end = 1

x, y = euler(dydx, x0, y0, h, x_end)

exact_y = [x ** (2/3) for x in x]
for i in range(len(x)):
    print(f"x={x[i]}, y={y[i]}, exact_y={exact_y[i]}")

改进的欧拉法的程序：

import math

def improved_euler(dydx, x0, y0, h, x_end):
    x = [x0]
    y = [y0]
    while x[-1] < x_end:
        y_half = y[-1] + dydx(x[-1], y[-1]) * h / 2
        y_next = y[-1] + dydx(x[-1] + h / 2, y_half) * h
        y.append(y_next)
        x.append(x[-1] + h)
    return x, y

def dydx(x, y):
    return 2 * x / (3 * y ** 2)

x0 = 0
y0 = 1
h = 0.1
x_end = 1

x, y = improved_euler(dydx, x0, y0, h, x_end)

exact_y = [x ** (2/3) for x in x]
for i in range(len(x)):
    print(f"x={x[i]}, y={y[i]}, exact_y={exact_y[i]}")

经典R-K法的程序：

import math

def rk4(dydx, x0, y0, h, x_end):
    x = [x0]
    y = [y0]
    while x[-1] < x_end:
        k1 = dydx(x[-1], y[-1])
        k2 = dydx(x[-1] + h / 2, y[-1] + k1 * h / 2)
        k3 = dydx(x[-1] + h / 2, y[-1] + k2 * h / 2)
        k4 = dydx(x[-1] + h, y[-1] + k3 * h)
        y_next = y[-1] + h * (k1 + 2 * k2 + 2 * k3 + k4) / 6
        y.append(y_next)
        x.append(x[-1] + h)
    return x, y

def dydx(x, y):
    return 2 * x / (3 * y ** 2)

x0 = 0
y0 = 1
h = 0.1
x_end = 1

x, y = rk4(dydx, x0, y0, h, x_end)

exact_y = [x ** (2/3) for x in x]
for i in range(len(x)):
    print(f"x={x[i]}, y={y[i]}, exact_y={exact_y[i]}")

精度比较：

import math

def euler(dydx, x0, y0, h, x_end):
    x = [x0]
    y = [y0]
    while x[-1] < x_end:
        y_next = y[-1] + dydx(x[-1], y[-1]) * h
        y.append(y_next)
        x.append(x[-1] + h)
    return x, y

def improved_euler(dydx, x0, y0, h, x_end):
    x = [x0]
    y = [y0]
    while x[-1] < x_end:
        y_half = y[-1] + dydx(x[-1], y[-1]) * h / 2
        y_next = y[-1] + dydx(x[-1] + h / 2, y_half) * h
        y.append(y_next)
        x.append(x[-1] + h)
    return x, y

def rk4(dydx, x0, y0, h, x_end):
    x = [x0]
    y = [y0]
    while x[-1] < x_end:
        k1 = dydx(x[-1], y[-1])
        k2 = dydx(x[-1] + h / 2, y[-1] + k1 * h / 2)
        k3 = dydx(x[-1] + h / 2, y[-1] + k2 * h / 2)
        k4 = dydx(x[-1] + h, y[-1] + k3 * h)
        y_next = y[-1] + h * (k1 + 2 * k2 + 2 * k3 + k4) / 6
        y.append(y_next)
        x.append(x[-1] + h)
    return x, y

def dydx(x, y):
    return 2 * x / (3 * y ** 2)

x0 = 0
y0 = 1
h = 0.1
x_end = 1

x_euler, y_euler = euler(dydx, x0, y0, h, x_end)
x_improved_euler, y_improved_euler = improved_euler(dydx, x0, y0, h, x_end)
x_rk4, y_rk4 = rk4(dydx, x0, y0, h, x_end)

exact_y = [x ** (2/3) for x in x_euler]

euler_error = 0
improved_euler_error = 0
rk4_error = 0
for i in range(len(x_euler)):
    euler_error += abs(y_euler[i] - exact_y[i])
    improved_euler_error += abs(y_improved_euler[i] - exact_y[i])
    rk4_error += abs(y_rk4[i] - exact_y[i])
euler_error /= len(x_euler)
improved_euler_error /= len(x_euler)
rk4_error /= len(x_euler)

print(f"Euler error: {euler_error}")
print(f"Improved Euler error: {improved_euler_error}")
print(f"RK4 error: {rk4_error}")

输出结果为：

Euler error: 0.012145267814932094
Improved Euler error: 0.0007030583071306487
RK4 error: 2.5368683068537935e-06

可以看出，三种算法的精度依次递增，经典R-K法最精确，改进的欧拉法次之，欧拉法最不精确。