numlab2

import numpy as np, matplotlib.pyplot as plt

def analytic_func(x, t):
    return np.exp(-t - x) * np.sin(x) * np.sin(2 * t)

def first_border(x, t):
    return 0

def second_border(x, t):
    return 0

def tridiagonal_alg(a, b, c, d, s):
    P = np.zeros(s + 1)
    Q = np.zeros(s + 1)

    P[0] = -c[0] / b[0]
    Q[0] = d[0] / b[0]

    k = s - 1
    for i in range(1, s):
        P[i] = -c[i] / (b[i] + a[i] * P[i - 1])
        Q[i] = (d[i] - a[i] * Q[i - 1]) / (b[i] + a[i] * P[i - 1])
    P[k] = 0
    Q[k] = (d[k] - a[k] * Q[k - 1]) / (b[k] + a[k] * P[k - 1])

    x = np.zeros(s)
    x[k] = Q[k]

    for i in range(s - 2, -1, -1):
        x[i] = P[i] * x[i + 1] + Q[i]

    return x

def implicit(K, t, tau, h, x):
    N = len(x)
    U = np.zeros((K, N))
    t += tau
    sigma = (tau/h) ** 2
    for j in range(N):
        U[0, j] = 0
        U[1, j] = tau * 2 * np.exp(-x[j]) * np.sin(x[j])

    for k in range(1, K - 1):
        a = np.zeros(N)
        b = np.zeros(N)
        c = np.zeros(N)
        d = np.zeros(N)
        t += tau

        for j in range(1, N - 1):
            a[j] = (1-h) * tau**2
            b[j] = -(2 + 3 * h**2) * tau**2 - (1 + tau) * h**2
            c[j] = (1+h) * tau**2
            d[j] = (U[k-1, j] - 2 * U[k, j] - tau * U[k-1, j]) * h**2

            b[0] = -(2 + 3 * h**2) * tau**2 - (1 + tau) * h**2
            c[0] = (1+h) * tau**2
            d[0] = first_border(x[j], t)

            a[N - 1] = (1-h) * tau**2
            b[N - 1] = -(2 + 3 * h**2) * tau**2 - (1 + tau) * h**2
            d[N - 1] = second_border(x[j], t)

        u_new = tridiagonal_alg(a, b, c, d, N)
        for i in range(N):
            U[k + 1, i] = u_new[i]

    return U

def explicit(K, t, tau, h, x, a, b, c, d):
    N = len(x)
    U = np.zeros((K, N))
    t += tau
    for j in range(N):
        U[0, j] = 0
        U[1][j] = tau * 2 * np.exp(-x[j]) * np.sin(x[j])

    for k in range(1, K - 1):
        t += tau
        for j in range(1, N - 1):
            U[k + 1, j] = ((2 + a * tau) * U[k, j] - U[k - 1, j] +
                           tau ** 2 * (b * (U[k, j + 1] - 2 * U[k, j] + U[k, j - 1]) / h ** 2 + c * (
                            U[k, j + 1] - U[k, j - 1]) / (2 * h) + d * U[k, j])) \
                          / (1 + a * tau)

            U[k + 1, 0] = first_border(x[0], t)
            U[k + 1, N - 1] = second_border(x[N - 1], t)

    return U

def main(N, K, time):
    a = 2
    b = 1
    c = 2
    d = -3
    h = np.pi / N
    tau = time / K
    x = np.arange(0, np.pi + h / 2 - 1e-4, h)
    T = np.arange(0, time, tau)
    t = 0
    U = explicit(K, t, tau, h, x, a, b, c, d)
    #U = implicit(K, t, tau, h, x)
    dt = int(input("Введите момент времени:"))
    U_analytic = analytic_func(x, T[dt])
    plt.title("График точного и численного решения задачи")
    plt.plot(x, U_analytic, label="Точное решение", color="red")
    plt.scatter(x, U[dt, :], label="Численное решение")
    error = abs(U_analytic - U[dt, :])
    plt.xlabel("x")
    plt.ylabel("U")
    plt.text(0.2, -0.25, "Максимальная ошибка метода: " + str(max(error)))
    plt.axis([-0.2, 3.3, -1, 1])
    plt.grid()
    plt.legend()
    plt.show()

    plt.title("График ошибки по шагам")
    error_time = np.zeros(len(T))
    for i in range(len(T)):
        error_time[i] = max(abs(analytic_func(x, T[i]) - U[i, :]))
    plt.plot(T, error_time, label="По времени")
    plt.legend()
    plt.grid()
    plt.show()

    h = []
    for i in range(0, 40):
        h.append(np.pi / (50 - i))
    er_ex = np.zeros(len(h))
    tau = time / K

    for i in range(len(h)):
        x = np.arange(0, np.pi + h[i] / 2 - 1e-4, h[i])

        if a * tau / h[i] ** 2 <= 0.5:
            Uex_1 = explicit(K, t, tau, h[i], x, a, b, c, d)

        er_ex[i] = max(abs(analytic_func(x, T[dt]) - Uex_1[dt, :]))

    print(er_ex)

    plt.title('Ошибка явной схемы')
    plt.xlabel('h')
    plt.ylabel('error')
    plt.plot(h, er_ex)
    plt.show()
    return 0

if __name__ == '__main__':
    N = 50
    K = 7000
    time = 3
    main(N, K, time)