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()