lab4_1

from copy import deepcopy
import numpy as np
import matplotlib.pyplot as plt

# Вторая производная, выраженная через x, f, f'
def ddf(x, f, df):
    return 4 * x * df - (4 * x ** 2 - 3) * f - np.e ** (x ** 2)

# Точное решение
def f(x):
    return (np.e ** x + np.e ** (-x) - 1) * np.e ** (x ** 2)

# Явный метод Эйлера
def Euler(ddy, borders, y0, z0, h):
    x = np.arange(borders[0], borders[1]+h, h)
    N = np.shape(x)[0]
    y = np.zeros(N)
    z = np.zeros(N)
    y[0] = y0; z[0] = z0
    for i in range(0, N-1):
        z[i+1] = z[i]+h*ddy(x[i], y[i], z[i])
        y[i+1] = y[i]+h*z[i]
    return y

# Явный метод Эйлера-Коши
def EulerCauchy(ddy, borders, y0, z0, h):
    x = np.arange(borders[0], borders[1]+h, h)
    N = np.shape(x)[0]
    y = np.zeros(N)
    z = np.zeros(N)
    y[0] = y0; z[0] = z0
    for i in range(0, N-1):
        pre_z = z[i]+h*ddy(x[i], y[i], z[i])
        pre_y = y[i]+h*z[i]
        z[i+1] = z[i]+h*(ddy(x[i], y[i], z[i])+ddy(x[i+1], pre_y, pre_z))/2
        y[i+1] = y[i]+h*z[i]
    return y

def RungeKutta(ddy, borders, y0, z0, h):
    x = np.arange(borders[0], borders[1] + h, h)
    N = np.shape(x)[0]
    y = np.zeros(N)
    z = np.zeros(N)
    y[0] = y0; z[0] = z0
    for i in range(N-1):
        K1 = h * z[i]
        L1 = h * ddy(x[i], y[i], z[i])
        K2 = h * (z[i] + 0.5 * L1)
        L2 = h * ddy(x[i] + 0.5 * h, y[i] + 0.5 * K1, z[i] + 0.5 * L1)
        K3 = h * (z[i] + 0.5 * L2)
        L3 = h * ddy(x[i] + 0.5 * h, y[i] + 0.5 * K2, z[i] + 0.5 * L2)
        K4 = h * (z[i] + L3)
        L4 = h * ddy(x[i] + h, y[i] + K3, z[i] + L3)
        delta_y = (K1 + 2 * K2 + 2 * K3 + K4) / 6
        delta_z = (L1 + 2 * L2 + 2 * L3 + L4) / 6
        y[i+1] = y[i] + delta_y
        z[i+1] = z[i] + delta_z
    return y, z

def Adams(ddy, borders, calc_y, calc_z, h):
    x = np.arange(borders[0], borders[1] + h, h)
    N = np.shape(x)[0]
    y = calc_y[:4]
    z = calc_z[:4]
    for i in range(3, N-1):
        # Этап: предиктор
        z_next = z[i] + (h/24) * (55 * ddy(x[i], y[i], z[i]) -
                          59 * ddy(x[i-1], y[i-1], z[i-1]) +
                          37 * ddy(x[i-2], y[i-2], z[i-2]) -
                          9 * ddy(x[i-3], y[i-3], z[i-3]))
        y_next = y[i] + (h/24)*(55*z[i]-59*z[i-1]+37*z[i-2]-9*z[i-3])

        # Этап: корректор
        z_i = z[i] + (h/24) * (9 * ddy(x[i+1], y_next, z_next) +
                          19 * ddy(x[i], y[i], z[i]) -
                          5 * ddy(x[i - 1], y[i - 1], z[i - 1]) +
                          1 * ddy(x[i - 2], y[i - 2], z[i - 2]))
        z = np.append(z, z_i)

        y_i = y[i] + (h/24)*(9*z_next + 19*z[i] - 5*z[i-1] + 1*z[i-2])
        y = np.append(y, y_i)
    return y

def sqr_error(y, y_correct):
    return np.sqrt(np.sum((y-y_correct)**2))

def RungeRombert(y1, y2, h1, h2, p):
    if h1 > h2:
        k = int(h1 / h2)
        y = np.zeros(np.shape(y1)[0])
        for i in range(np.shape(y1)[0]):
            y[i] = y2[i*k]+(y2[i*k]-y1[i])/(k**p-1)
        return y
    else:
        k = int(h2 / h1)
        y = np.zeros(np.shape(y2)[0])
        for i in range(np.shape(y2)[0]):
            y[i] = y1[i * k] + (y1[i * k] - y2[i]) / (k ** p - 1)
        return y

borders = (0, 1)
y0 = 1
z0 = 1
h = 0.1
x = np.arange(borders[0], borders[1]+h, h)
y = f(x)
y1 = Euler(ddf, borders, y0, z0, h)
y2 = EulerCauchy(ddf, borders, y0, z0, h)
y3, z = RungeKutta(ddf, borders, y0, z0, h)
y4 = Adams(ddf, borders, deepcopy(y3), z, h)
# Получаем решения методов для шага h/2
h2 = h/2
y1_2 = Euler(ddf, borders, y0, z0, h2)
y2_2 = EulerCauchy(ddf, borders, y0, z0, h2)
y3_2, z2 = RungeKutta(ddf, borders, y0, z0, h2)
y4_2 = Adams(ddf, borders, deepcopy(y3_2), z2, h2)
print("Оценка погрешности методом Рунге-Ромберта:")
print("Для явного метода Эйлера:", sqr_error(y1, RungeRombert(y1, y1_2, h, h2, 1)))
print("Для явного метода Эйлера-Коши:", sqr_error(y2, RungeRombert(y2, y2_2, h, h2, 1)))
print("Для метода Рунге-Кутты:", sqr_error(y3, RungeRombert(y3, y3_2, h, h2, 4)))
print("Для метода Адамса:", sqr_error(y4, RungeRombert(y4, y4_2, h, h2, 4)))
print("\nСравнение с точным решением:")
print("Для явного метода Эйлера:", sqr_error(y1, y))
print("Для явного метода Эйлера-Коши:", sqr_error(y1, y))
print("Для метода Рунге-Кутты:", sqr_error(y3, y))
print("Для метода Адамса:", sqr_error(y4, y))
plt.figure(figsize=(12, 7))
plt.plot(x, y, label='Точное решение', linewidth=1, color="green")
plt.plot(x, y1, label='Явный метод Эйлера', linewidth=1, color="yellow")
plt.plot(x, y2, label='Явный метод Эйлера-Коши', linewidth=1, color="orange")
plt.plot(x, y3, label='Метод Рунге-Кутты', linewidth=1, color="red")
plt.plot(x, y4, label='Метод Адамса', linewidth=1, color="blue")
plt.grid()
plt.legend()
plt.show()