#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Wed May 20 08:11:56 2026

@author: gb
"""

import numpy as np
import matplotlib.pyplot as plt

"""
Exercice 1
"""
# Q1. 
M = np.zeros((4, 4))
## Les commandes sont à taper dans la console (voir corrigé pdf).

# Q3.
def F(n): # en remplissant une matrice de 0.
    M = np.zeros((n, n), dtype = complex)
    for i in range(n):
        for j in range(n):
            M[i, j] = np.exp(2j * i * j * np.pi/n)
    return M

def F_bis(n): # par conversion d'une liste de listes en compréhension
    L = [[np.exp(2j * i * j * np.pi/n) for j in range(n)] for i in range(n)]
    return np.array(L)

def F_ter(n): # par conversion, puis redimensionnement d'une liste
    L = [np.exp(2j * i * j * np.pi/n) for i in range(n) for j in range(n)]
    return np.array(L).reshape(n, n)

## dans la console : F(3)

# Q4
import numpy.linalg as alg

M = F(5)
vp, VP = alg.eig(M)
D = np.diag(vp)
Delta = alg.inv(VP).dot(M).dot(VP) - D
d = abs(alg.det(VP))

# Q5
M = np.ones((5, 5))
vp, VP = alg.eig(M)
spectre = set(vp) 

## complément
A = np.zeros((2, 2))
A[0, 1] = 1
B = np.zeros((5, 5))
for i in range(4):
    B[i, i+1] = 1
vpa, VPA = alg.eig(A)
vpb, VPB = alg.eig(B)
Da = alg.inv(VPA).dot(A).dot(VPA) 
Db = alg.inv(VPB).dot(A).dot(VPB)
da, db = abs(alg.det(VPA)), abs(alg.det(VPB))



"""
Exercice 2
"""
# Q2.
import scipy.integrate as sci

def f(a, t):
    return np.exp(-t*t - a**2 / t**2)

def I(a):
    return sci.quad(lambda t: f(a, t), 0, np.inf)[0]

def aJ(a):
    return a * sci.quad(lambda t: f(a, t)/t**2, 0, np.inf)[0]

JJ = np.vectorize(aJ)

A = np.linspace(0.1, 4, 80)
B = [I(a) for a in A]
C = JJ(A)
D = abs(B - C) / B
plt.figure()
plt.title("$I(a)$ et $aJ(a)$")
plt.plot(A, B, label = "$I(a)$")
plt.plot(A, C, label = "$aJ(a)$")
plt.savefig('/Users/mbpadmin/Desktop/PSI/Centrale - oral 2/Prepa 2026/Mansouri-1.png')
plt.legend()
plt.show()

plt.figure()
plt.title("$I(a)$ vs. $aJ(a)$, erreur relative")
plt.plot(A, D)
plt.savefig('/Users/mbpadmin/Desktop/PSI/Centrale - oral 2/Prepa 2026/Mansouri-2.png')
plt.show()



"""
Exercice 3
"""
from numpy.polynomial import Polynomial

# Q1.
def can(n):
    return [Polynomial([0]*k + [1]) for k in range(n+1)]
    
# Q2.
def Tchebychev(n):
    T = [Polynomial([1]), Polynomial([0, 1])]
    for k in range(n-1):
        T.append(2*T[1] * T[-1] - T[-2])
    return T

# Q3.
# import scipy.integrate as sci

def produit_scalaire(P, Q):
    return sci.quad(lambda t: P(t) * Q(t) / np.sqrt(1 - t*t), -1, 1)[0]

# Q4. 
T, maxi = Tchebychev(20), 0
for i in range(20):
    for j in range(i+1, 21):
        maxi = max(maxi, abs(produit_scalaire(T[i], T[j])))
print(maxi) # à faire en principe dans la console

# Q5.
def Gram_Schmidt(n):
    CAN = can(n)
    BON = []
    for i in range(n+1):
        proj = CAN[i]
        for j in range(i):
            proj -= produit_scalaire(CAN[i], BON[j]) * BON[j]
        BON.append(proj / np.sqrt(produit_scalaire(proj, proj)))
    return BON

# Q6.
BON = Gram_Schmidt(20)
BON2 = [P / np.sqrt(produit_scalaire(P, P)) for P in Tchebychev(20)]

Diff2 = []
for k in range(len(BON)):
    Diff2.append(np.amax(abs(BON[k].coef - BON2[k].coef)))

# Q7.
## Vérification formelle
P, T, L = Polynomial([1, 0, -1]), Tchebychev(20), []
for n in range(21):
    Q = P * T[n].deriv(2) - T[1] * T[n].deriv() + n**2 * T[n]
    L.append(np.amax(Q.coef))

## Vérification numérique
def f(n, x, t):
    return np.array([x[1], (t * x[1] - n**2 * x[0])/(1 - t**2)])

V = []
for n in range(21):
    tps = np.arange(0, 0.99, 0.01)
    a, b = np.cos(n * np.pi / 2), n * np.sin(n * np.pi/2)
    X = sci.odeint(lambda x, t: f(n, x, t), np.array([a, b]), tps)
    V.append(X[:, 0])
  
W = [abs(T[k](tps) - V[k]) for k in range(21)]
Err = [np.amax(W[k]) for k in range(21)]



"""
Exercice 4
"""
# Q1.
def matrice_extraite(A, i, j): # en utilisant une concaténation de blocs
    B, C, D, E = A[:i, :j], A[:i, j+1:], A[i+1:, :j], A[i+1:, j+1:]
    F, G = np.concatenate((B, C), axis = 1), np.concatenate((D, E), axis = 1)
    return np.concatenate((F, G), axis = 0)

def matrice_extraite2(A, i, j): # en remplissant une matrice de la bonne taille
    n, p = np.shape(A)
    M = np.zeros((n-1, p-1))
    M[:i, :j], M[:i, j:] = A[:i, :j], A[:i, j+1:], 
    M[i:, :j], M[i:, j:] = A[i+1:, :j], A[i+1:, j+1:]
    return M

# Q2.
def determinant_rec(A):
    n = A.shape[0]
    if n == 1:
        return A[0, 0]
    D = 0
    for j in range(n):
        if A[0, j] != 0:
            D += (-1)**(j+1) * A[0, j] * determinant_rec(matrice_extraite(A, 0, j))
    return D

from random import random

def verification(): # cette fonction ne renvoie rien, elle affiche le résultat
    for k in range(5):
        B = np.array([8*random() + k/3 for k in range(25)]).reshape(5, 5)
        d = alg.det(B)
        print(d, d - determinant_rec(B))



"""
Exercice 5
"""
# Q2.
import scipy.integrate as sci
import matplotlib.pyplot as plt

A, omega = [k*0.5 for k in range(1, 6)], 1

def f(x, t): # on vectorialise l'ED, qui devient une ED d'ordre 1
    return np.array([x[1], - omega**2 * np.sin(x[0])])
T = np.linspace(0, 10, 100)

plt.figure()
for a in A:
    X = sci.odeint(f, np.array([0, a]), T)
    plt.plot(T, X[:, 0])
plt.title('Equation du pendule pour différentes vitesses angulaires initiales')    
plt.show()

# Q3c.
TF, r = np.arange(-5, 5.01, 0.01), 0.9 

def pend(u):
    return 1/np.sqrt(1 - r**2 * np.sin(omega*u)**2)

def F(t):
    return sci.quad(pend, 0, t)[0]

XF = [F(t) for t in TF]
plt.plot(TF, XF)
plt.show()


# Q4.
import scipy.optimize as sco

borne = sco.fsolve(lambda t: F(t) - 5, 5.)[0]
TF2 = np.arange(-borne, borne, 0.02)
XF2 = [F(t) for t in TF2]
YF2 = [2*np.arcsin(r*np.sin(omega*u)) for u in TF2]
plt.figure()
plt.plot(XF2, YF2)
plt.title('Le pendule non linéarisé')
plt.show()