ISC_421_Controle_4_Saroukhanian_Iliya.py

import numpy as np
import math
from matplotlib import pyplot as plt
from scipy.interpolate import lagrange, barycentric_interpolate
from numpy.polynomial.polynomial import Polynomial


#############################################################################
######################### Variables globales ################################
#############################################################################

Nmbre_derivees = 14
Nmbre_pts = 201


def Base_function(x): return np.sin(np.pi*np.cos(4*x*x))


Taylor_points = np.array([-0.8, 0.25, 0.8])

Borne_a = 1
Borne_b = 4


def sigma(a, b, x): return (b-a)/2*x + (b+a)/2


def amgis(a, b, x): return -1 + 2*(x-a)/(b-a)


#############################################################################
############ Fonctions necessaires a la production des donnees ##############
#############################################################################

def deltas2(lstx, lsty, depth=1):
    N = len(lstx)
    deltas_dic = {tuple([l]): y for l, y in enumerate(lsty)}
    for k in range(1, min(depth, N)):
        for l in range(N - k):
            indices_now = tuple(range(l, k+l+1))
            deltas_dic[indices_now] = (
                deltas_dic[indices_now[1:]] - deltas_dic[indices_now[:-1]])/(lstx[k+l]-lstx[l])
    return deltas_dic


def approx_der_lst(lstx, lsty, a, depth):
    D = deltas2(lstx, lsty, depth=depth)
    indices = [l for l, x in enumerate(lstx) if x >= a]
    index = indices[0]
    lst_return = []
    l = 0
    while (index - l >= 0) and (index + l < len(lstx)) and (l < depth):
        der_temp = []
        indices_now = range(index - l, index + 1)
        for m in indices_now:
            approx_der_temp = D[tuple(range(m, m+l+1))]*np.math.factorial(l)
            der_temp.append(approx_der_temp)
        approx_der = sum(der_temp)/(len(der_temp))
        lst_return.append(approx_der)
        l += 1

    return lst_return


#############################################################################
##### Classe regroupant toutes les donnees pertinentes d'une fonction #######
#############################################################################

class Function_data:

    def __init__(self, f=lambda x: x):
        self.function = f

        def approx_der_temp(pt, width=.1, number_points=Nmbre_derivees):
            t = np.linspace(pt - width, pt + width, 2*number_points + 1)
            return approx_der_lst(t, f(t), pt, number_points)

        def max_der(a, b, width=.1, number_points=12, total_number=Nmbre_derivees):
            ts = np.linspace(a, b, total_number)
            max_ret = [0]*(number_points-1)
            for p in ts:
                der = approx_der_temp(p)
                max_ret = [max(oldmax, abs(der[l]))
                           for l, oldmax in enumerate(max_ret)]
            return max_ret
        self.approx_derivatives = approx_der_temp
        self.values = lambda t: f(t)
        self.max_der = max_der


#############################################################################
######## Classe regroupant toutes les donnees fournies a l'etudiant #########
#############################################################################

class Student_data:

    def __init__(self):

        self.a, self.b = Borne_a, Borne_b
        def f_ab(x): return Base_function(amgis(Borne_a, Borne_b, x))
        F = Function_data(f_ab)
        self.Maximal_derivatives_values = F.max_der(
            Borne_a, Borne_b, number_points=Nmbre_derivees, total_number=Nmbre_pts)
        self.Taylor_points = [sigma(Borne_a, Borne_b, p)
                              for p in Taylor_points]
        self.Taylor_derivatives_values = {pt: F.approx_derivatives(
            pt, number_points=Nmbre_derivees) for pt in self.Taylor_points}
        self.f = f_ab


def compute_taylor_series(x, pt, derivs):
    sum = 0

    for nth_deg, nth_deriv in enumerate(derivs):
        sum += (nth_deriv / math.factorial(nth_deg)) * (x - pt)**nth_deg

    return sum

###########################################################################
#### Exemple d'utilisation de la classe, avec explication des donnees #####
###########################################################################


SD = Student_data()

print(f"bornes de l'intervalle considere: [{SD.a},{SD.b}]")
print()
print(f"Points ou on connait la derivee : {SD.Taylor_points}")
print(f"Valeurs des derivees 0 a {Nmbre_derivees-1}:")
for x in SD.Taylor_points:
    print(f"x = {x}, dérivées: {SD.Taylor_derivatives_values[x]}")
print()
print(f"Valeurs maximales des derivees 0 a {
      Nmbre_derivees-1} sur l'intervalle [{SD.a},{SD.b}]: {SD.Maximal_derivatives_values}")
print()
print("La fonction est dans la variable SD.f, si on doit l'évaluer en certains points")
print()

print(f"valeur de la fonction en x = {SD.Taylor_points} : {
      [SD.f(x) for x in SD.Taylor_points]}")
print(f"valeur de la fonction en a et b:  {SD.f(SD.a), SD.f(SD.b)}")

# Exemple de graphe de la fonction f.
# t = np.linspace(SD.a, SD.b, Nmbre_pts)
# y = compute_taylor_series(
#     t, SD.Taylor_points[1], SD.Taylor_derivatives_values[SD.Taylor_points[1]])

# plt.plot(t, SD.f(t), "k")
# plt.show()


###########################################################################
#### Exercices #####
###########################################################################


def errmax_taylor(plot_range, eval_pt):
    return (SD.Maximal_derivatives_values[-1] /
            math.factorial(len(SD.Maximal_derivatives_values) + 1)) * \
        (plot_range - eval_pt)**(len(SD.Maximal_derivatives_values) + 1)


def plot_taylor_poly():
    t = np.linspace(SD.a, SD.b, Nmbre_pts)
    y0 = compute_taylor_series(
        t, SD.Taylor_points[0], SD.Taylor_derivatives_values[SD.Taylor_points[0]])

    err_y0 = errmax_taylor(t, SD.Taylor_points[0])

    y1 = compute_taylor_series(
        t, SD.Taylor_points[1], SD.Taylor_derivatives_values[SD.Taylor_points[1]])

    err_y1 = errmax_taylor(t, SD.Taylor_points[1])

    y2 = compute_taylor_series(
        t, SD.Taylor_points[2], SD.Taylor_derivatives_values[SD.Taylor_points[2]])

    err_y2 = errmax_taylor(t, SD.Taylor_points[2])

    fig, axs = plt.subplots(2, 2, figsize=(18, 12))

    axs[0, 0].plot(t, SD.f(t), color='black')
    axs[0, 0].set_title('Fonction $f$ originale')
    axs[0, 0].grid()

    axs[0, 1].plot(t, SD.f(t), color='black', label='$f$')
    axs[0, 1].plot(t, y0, color='orange', label='$T_{f}$')
    axs[0, 1].plot(t, err_y0, '--', color='orange',
                   label='Erreur maximal de $T_{f}$')
    axs[0, 1].plot(SD.Taylor_points[0], SD.f(
        SD.Taylor_points[0]), "-o", color='red', label=f'a = {SD.Taylor_points[0]}')
    axs[0, 1].set_title(
        f'Développement de $T_f$ en $a = {SD.Taylor_points[0]}$')
    axs[0, 1].set_ylim([-1.2, 1.2])
    axs[0, 1].grid()
    axs[0, 1].legend()

    axs[1, 0].plot(t, SD.f(t), color='black', label='$f$')
    axs[1, 0].plot(t, y1, color='blue', label='$T_{f}$')
    axs[1, 0].plot(t, err_y1, '--', color='blue',
                   label='ErrMax de $T_{f}$')
    axs[1, 0].plot(SD.Taylor_points[1], SD.f(
        SD.Taylor_points[1]), "-o", color='red', label=f'a = {SD.Taylor_points[1]}')
    axs[1, 0].set_title(
        f'Développement de $T_f$ en $a = {SD.Taylor_points[1]}$')
    axs[1, 0].set_ylim([-1.2, 1.2])
    axs[1, 0].grid()
    axs[1, 0].legend()

    axs[1, 1].plot(t, SD.f(t), color='black', label='$f$')
    axs[1, 1].plot(t, y2, color='violet', label='$T_{f}$')
    axs[1, 1].plot(t, err_y2, '--', color='violet',
                   label='ErrMax de $T_{f}$')
    axs[1, 1].plot(SD.Taylor_points[2], SD.f(
        SD.Taylor_points[2]), "-o", color='red', label=f'a = {SD.Taylor_points[2]}')
    axs[1, 1].set_title(
        f'Développement de $T_f$ en $a = {SD.Taylor_points[2]}$')
    axs[1, 1].set_ylim([-1.2, 1.2])
    axs[1, 1].grid()
    axs[1, 1].legend()

    fig.tight_layout()
    plt.show()


def errmax_poly_interpolation(nb_points, interpolation_pts, plot_range):
    errs_range = []

    for pt in plot_range:
        prod = 1

        for inter_pt in interpolation_pts:
            prod *= pt - inter_pt

        err = (SD.Maximal_derivatives_values[nb_points] /
               math.factorial(nb_points + 1)) * np.abs(prod)

        errs_range.append(err)

    return errs_range


def chebyshev_pts(start, stop, nb_points):
    return (((start + stop) / 2) + ((stop - start) / 2) *
            np.cos(((2 * np.arange(nb_points) + 1) * np.pi) / (2 * (nb_points))))


def plot_lagrange_poly():
    nb_points = np.linspace(1, 12, 6, dtype=np.uint8)
    fig, axes = plt.subplots(2, 3, figsize=(20, 12))

    t = np.linspace(SD.a, SD.b, Nmbre_pts)

    for i, ax in enumerate(axes.flat):
        cheb_pts = chebyshev_pts(SD.a, SD.b, nb_points[i])

        interpolate_pts = np.linspace(SD.a, SD.b, nb_points[i])

        l_poly_uniform = lagrange(interpolate_pts, SD.f(interpolate_pts))
        l_poly_chebyshev_pts = lagrange(
            cheb_pts, SD.f(cheb_pts))

        ax.plot(t, SD.f(t), color='black', label='f')
        ax.plot(t, l_poly_uniform(t), color='red',
                label='$L_{f}$, intervalle équidistants')
        # ax.plot(t, np.abs(SD.f(t) - l_poly_uniform(t)), '--', color='red',
        #         label='$L_{f}$, intervalle équidistants, erreur')
        ax.plot(t, l_poly_chebyshev_pts(t), color='blue',
                label='$L_{f}$, points de Chebyshev')
        # ax.plot(t, np.abs(SD.f(t) - l_poly_chebyshev_pts(t)), '--', color='blue',
        #         label='$L_{f}$, points de Chebyshev, erreur')
        ax.plot(interpolate_pts, SD.f(interpolate_pts), 'o', color='red',
                label='Points équidistants')
        ax.plot(cheb_pts, SD.f(cheb_pts), 'o', color='blue',
                label='Points de Chebyshev')
        ax.set_title(f'n = {nb_points[i]}')
        ax.set_ylim([-1.2, 1.2])

        ax.legend()

    fig.suptitle(f'Polynômes d\'interpolation de Lagrange de $f$ avec 2 subdivisions différentes d\'intervalle: Équidistantes (rouge) / Points de Chebyshev (bleu)')

    fig.tight_layout()
    plt.show()


def plot_errmax_interpolation():
    nb_points = np.linspace(1, 12, 6, dtype=np.uint8)
    fig, axes = plt.subplots(2, 3, figsize=(20, 12))

    t = np.linspace(SD.a, SD.b, Nmbre_pts)

    for i, ax in enumerate(axes.flat):
        cheb_pts = chebyshev_pts(SD.a, SD.b, nb_points[i])

        interpolate_pts = np.linspace(SD.a, SD.b, nb_points[i])

        l_poly_uniform = lagrange(interpolate_pts, SD.f(interpolate_pts))
        l_poly_chebyshev_pts = lagrange(
            cheb_pts, SD.f(cheb_pts))

        uniform_err = errmax_poly_interpolation(
            nb_points[i], interpolate_pts, t)
        chebyshev_err = errmax_poly_interpolation(nb_points[i], cheb_pts, t)

        ax.plot(t, SD.f(t), color='black', label='f')
        ax.plot(t, l_poly_uniform(t), color='red',
                label='$L_{f}$, intervalle équidistants')
        ax.plot(t, uniform_err, '--', color='red',
                label='ErrMax de $L_{f}$, intervalle équidistants')
        ax.plot(t, l_poly_chebyshev_pts(t), color='blue',
                label='$L_{f}$, points de Chebyshev')
        ax.plot(t, chebyshev_err, '--', color='blue',
                label='ErrMax de $L_{f}$, points de Chebyshev')
        ax.plot(interpolate_pts, SD.f(interpolate_pts), 'o', color='red',
                label='Points équidistants')
        ax.plot(cheb_pts, SD.f(cheb_pts), 'o', color='blue',
                label='Points de Chebyshev')
        ax.set_title(f'n = {nb_points[i]}')

        ax.legend()

    fig.suptitle(f'Erreur maximale théorique lors de l\'interpolation de $f$ avec 2 subdivisions différentes d\'intervalle: Équidistantes (rouge) / Points de Chebyshev (bleu)')

    fig.tight_layout()
    plt.show()


plot_taylor_poly()
plot_lagrange_poly()
plot_errmax_interpolation()