Untitled

mail@pastecode.io avatar
unknown
python
2 years ago
2.4 kB
6
Indexable
Never
import numpy as np
import matplotlib.pyplot as plt


def f_x(x):
    return (np.cos(x) - 1) / x


def g_x(x):
    return np.tan(x)


def find_integral(f, begin, end):
    n = 1000
    x = np.linspace(begin, end, n)
    # if begin == 0:
    #     x = x[1:]
    if end == np.pi/2:
        x = x[:len(x)-1]
    return sum(abs(f(x))/n)


def next_element(integr, x, n):
    return - integr * x * x * n / ((n + 1) * (2 * n + 1) * (2 * n + 2))


def find_value(x, eps):
    n = 1
    next_m = - x * x / 4
    summ = next_m
    j = 1
    while abs(next_m) >= eps:
        next_m = next_element(next_m, x, n)
        n += 1
        j += 1
        summ += next_m
    return summ, j


def first_task(x, eps):
    res_s, j_s, int_s = [], [], []
    for i in x:
        res, j = find_value(i, eps)
        res_s.append(res)
        j_s.append(j)
    # plt.plot(x, res_s)
    # plt.suptitle("Ci(x)")
    # plt.xlabel("x")
    # plt.ylabel("Ci(x)")
    # # plt.plot(x, int_s)
    # plt.show()
    return res_s


def lagranj(x, y, t):
    z = np.zeros_like(t)
    for j in range(len(y)):
        p1 = np.ones_like(t)
        p2 = np.ones_like(t)
        for i in range(len(x)):
            if i == j:
                p1 = p1*1
                p2 = p2*1
            else:
                p1 = p1*(t-x[i])
                p2 = p2*(x[j]-x[i])
        z = z + y[j] * p1 / p2
    return z

def task_2(bounds, eps, n, m):
    x = np.linspace(bounds[0], bounds[1], m)
    y = first_task(x, eps)
    t = np.linspace(bounds[0], bounds[1], n)
    res = lagranj(x, y, t)
    acc_res = first_task(t, eps)
    error = max(abs(res - acc_res))
    #print(error)
    #chebyshev_t = [((bounds[1]-bounds[0])/2)*np.cos((2*i+1)*np.pi/(2*n+2))+(bounds[1]+bounds[0])/2 for i in range(n)]
    #chebyshev_res = lagranj(x, y, chebyshev_t)
    #chebyshev_res = np.flip(chebyshev_res)
    #chebyshev_error = max(abs(chebyshev_res - acc_res))
    #plt.plot(chebyshev_t, chebyshev_res)
    #plt.plot(chebyshev_t,chebyshev_res)
    #plt.show()
    #print(chebyshev_error)
    return error




if __name__ == "__main__":
    bounds = [0.4, 4]
    eps = 1e-6
    n = 12
    m = 10
    err_s = []
    ch_err_s = []
    m_s = []
    while m < 1000:
        err = task_2(bounds, eps, n, m)
        err_s.append(err)
        m_s.append(m)
        m += 50

    plt.plot(m_s, err_s)
    plt.show()