Загрузка данных


import math
import matplotlib.pyplot as plt

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

def exact(x):
    return x**2 + math.sqrt(x)
def euler(x, y, h):
    return y + h * f(x, y)

def rk4(x, y, h):
    k1 = f(x, y)
    k2 = f(x + h/2, y + h*k1/2)
    k3 = f(x + h/2, y + h*k2/2)
    k4 = f(x + h, y + h*k3)
    return y + h * (k1 + 2*k2 + 2*k3 + k4) / 6

h = 0.1
x = [1.0 + i*h for i in range(6)]
y_rk = [2.0]
y_euler = [2.0]

for i in range(3):
    y_rk.append(
        rk4(x[i], y_rk[i], h)
    )
    y_euler.append(
        euler(x[i], y_euler[i], h)
    )


for i in range(3, 5):

    y_new_rk = y_rk[i] + h / 24 * (
        55 * f(x[i], y_rk[i])
        - 59 * f(x[i-1], y_rk[i-1])
        + 37 * f(x[i-2], y_rk[i-2])
        - 9 * f(x[i-3], y_rk[i-3])

    )

    y_new_euler = y_euler[i] + h / 24 * (
        55 * f(x[i], y_euler[i])
        - 59 * f(x[i-1], y_euler[i-1])
        + 37 * f(x[i-2], y_euler[i-2])
        - 9 * f(x[i-3], y_euler[i-3])

    )

    y_rk.append(y_new_rk)
    y_euler.append(y_new_euler)



y_exact = [exact(xi) for xi in x]



error_rk = [

    abs(y_rk[i] - y_exact[i])

    for i in range(len(x))

]

error_euler = [

    abs(y_euler[i] - y_exact[i])

    for i in range(len(x))

]




print(

    "k   x     Адамс+РК    Адамс+Эйлер   "

    "точное y    ошибка РК    ошибка Эйлера"

)

for i in range(len(x)):

    print(

        i,

        round(x[i], 1),

        round(y_rk[i], 6),

        round(y_euler[i], 6),

        round(y_exact[i], 6),

        round(error_rk[i], 8),

        round(error_euler[i], 8)

    )

# График решений

plt.figure()

plt.plot(

    x, y_rk,

    "o-",

    label="Адамс + Рунге-Кутта"

)

plt.plot(

    x, y_euler,

    "o-",

    label="Адамс + Эйлер"

)

plt.plot(

    x, y_exact,

    "s--",

    label="Точное решение"

)

plt.xlabel("x")

plt.ylabel("y")

plt.title("Сравнение методов")

plt.grid()

plt.legend()

plt.show()

# График погрешностей

plt.figure()

plt.plot(

    x, error_rk,

    "o-",

    label="Ошибка: Рунге-Кутта"

)

plt.plot(

    x, error_euler,

    "o-",

    label="Ошибка: Эйлер"

)

plt.xlabel("x")

plt.ylabel("Погрешность")

plt.title("Сравнение погрешностей")

plt.grid()

plt.legend()

plt.show()