lab4_2

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

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

# Функция перед y'
def p(x):
    return -np.tan(x)

# Функция перед y
def q(x):
    return 2

# Функция в правой части уравнения
def right_f(x):
    return 0

def diff_right(BCondition, h, y):
    return (BCondition['c']+(BCondition['b']/h)*y[-2]) / (BCondition['a']+(BCondition['b']/h))

def ShotingMethod(ddy, borders, eta1, eta2, BCondition1, BCondition2, h):
    y0 = BCondition1['c']/BCondition1['a']
    eps = 0.1*h
    resolve1 = RungeKutta(ddy, borders, y0, eta1, h)[0]
    resolve2 = RungeKutta(ddy, borders, y0, eta2, h)[0]
    Phi1 = resolve1[-1] - diff_right(BCondition2, h, resolve1)
    Phi2 = resolve2[-1] - diff_right(BCondition2, h, resolve2)
    x = np.arange(borders[0], borders[1] + h, h)
    plt.figure(figsize=(12, 7))
    plt.plot(x, resolve1, '--', label='1 приближение', linewidth=1)
    plt.plot(x, resolve2, '--', label='2 приближение', linewidth=1)
    iter=2;
    while abs(Phi2 - Phi1) > eps:
        iter = iter + 1
        temp = eta2
        eta2 = eta2 - (eta2 - eta1) / (Phi2 - Phi1) * Phi2
        eta1 = temp
        resolve1 = RungeKutta(ddy, borders, y0, eta1, h)[0]
        resolve2 = RungeKutta(ddy, borders, y0, eta2, h)[0]
        Phi1 = resolve1[-1] - diff_right(BCondition2, h, resolve1)
        Phi2 = resolve2[-1] - diff_right(BCondition2, h, resolve2)
        plt.plot(x, resolve2, '--', label=f'{iter} приближение', linewidth=1)

    plt.grid()
    plt.legend()
    plt.show()
    return RungeKutta(ddy, borders, y0, eta2, h)[0]

def FiniteDifferenceMethod(BCondition1, BCondition2, equation, borders, h):
    x = np.arange(borders[0], borders[1] + h, h)
    N = np.shape(x)[0]

    # Составляем СЛАУ
    A = np.zeros((N, N))
    b = np.zeros(N)
    A[0][0] = BCondition1['a'] - BCondition1['b']/h
    A[0][1] = BCondition1['b']/h
    b[0] = BCondition1['c']
    for i in range(1, N-1):
        A[i][i-1] = 1/h**2 - equation['p'](x[i])/(2*h)
        A[i][i] = -2/h**2 + equation['q'](x[i])
        A[i][i+1] = 1/h**2 + equation['p'](x[i])/(2*h)
        b[i] = equation['f'](x[i])
    A[N-1][N-2] = -BCondition2['b']/h
    A[N-1][N-1] = BCondition2['a'] + BCondition2['b']/h
    b[N-1] = BCondition2['c']

    return Solve(A, b)

# Метод Рунге-Кутты решения задачи Коши
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 Solve(A, b):
    p = np.zeros(len(b))
    q = np.zeros(len(b))

    # Прямой ход: поиск прогоночных коэффициентов P и Q
    p[0] = -A[0][1]/A[0][0]
    q[0] = b[0]/A[0][0]
    for i in range(1, len(p)-1):
        p[i] = -A[i][i+1]/(A[i][i] + A[i][i-1]*p[i-1])
        q[i] = (b[i] - A[i][i-1]*q[i-1])/(A[i][i] + A[i][i-1]*p[i-1])

    p[-1] = 0
    q[-1] = (b[-1] - A[-1][-2]*q[-2])/(A[-1][-1] + A[-1][-2]*p[-2])

    # Обратный ход: поиск x
    x = np.zeros(len(b))
    x[-1] = q[-1]
    for i in reversed(range(len(b)-1)):
        x[i] = p[i]*x[i+1] + q[i]
    return x

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

equation = {'p': p, 'q': q, 'f': right_f}
BCondition1 = {'a': 1, 'b': 0, 'c': 2}
BCondition2 = {'a': 1, 'b': 0, 'c': 2.5 - 0.5*np.log(3)}
borders = (0, np.pi/6)
h = np.pi/120
x = np.arange(borders[0], borders[1]+h, h)
y = f(x)
eta1 = abs(y[-1]-y[1])/abs(x[-1]-x[1])
eta2 = eta1/2
y1 = ShotingMethod(ddf, borders, eta1, eta2, BCondition1, BCondition2, h)
y2 = FiniteDifferenceMethod(BCondition1, BCondition2, equation, borders, h)
h2 = h / 2
y1_2 = ShotingMethod(ddf, borders, eta1, eta2, BCondition1, BCondition2, h2)
y2_2 = FiniteDifferenceMethod(BCondition1, BCondition2, equation, borders, 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("\nСравнение с точным решением:")
print("Для метода стрельбы:", sqr_error(y1, y))
print("Для метода конечных разностей:", sqr_error(y2, y))
plt.figure(figsize=(12, 7))
plt.plot(x, y, label='Точное решение', linewidth=1, color="blue")
plt.plot(x, y1, label='Метод стрельбы', linewidth=1, color="red")
plt.plot(x, y2, label='Метод конечных разностей', linewidth=1, color="green")
plt.grid()
plt.legend()
plt.show()