"""
Thème 4 : calcul intégral
"""

## Importations
import numpy as np
import matplotlib.pyplot as plt
###Q1


# Faites un dessin à la main

# rem : pas de conclusion sur d'éventuelles considérations
# de sous-estimations de l'aire par les rectangles, sauf,
# bien entendu, si la fonction est monotone sur [a,b]


## Q2

def Riemann(f,a,b,n):
    h = (b-a)/n   # largeur du rectangle/ pas de la méthode des rectangles
    S  = 0        # initialisation du calcul de somme
    for k in range(n):
        S+=h*f(a+k*h)
    return S


# Rem : on peut boucler directement sur la liste des xk :

def Riemann(f,a,b,n):
    h = (b-a)/n
    Lx=[a + k*h for k in range(n)] # liste en compréhension
    S = 0
    for xk in Lx:
        S+=h*f(xk)
    return S

## Q3a)
# Vérif du programme pour "intégrale de 0 à x de x dx" :
def f1(x):
    return x

a,b,n=0,1,100
print(Riemann(f1,a,b,n)) # vaut environ 0.5


# test en augmentant le nombre de rectangles
for n in [10,20,30,40,100,150,1000]:
    An = Riemann(f1,a,b,n)
    print(f"Avec {n} rectangles,l'aire estimée est : {An}")



## Q3b)
# Si on essaie de calculer l'intégrale de a à b, avec
# a>b, on obtient -intégrale de b à a, puisque h <0. De plus,
# la somme se calcule  sur les valeurs f(b), f(b+h), etc.
# Ainsi le calcul de cette somme est celui du calcul
# l'intégrale de b à a  calculée avec la méthode des
# rectangles  **à droite**.
# Autrement dit, on n'a pas numériquement l'égalité :
# Riemann(a,b)=-Rimeann(b,a).

## Q4
def g(x):
    return 4/(1+x**2)

# rem : g s'intègre en 4Arctan(x) donc
#       son intégrale sur [0,1] vaut pi.

for p in range(1,1+6):
    n = 10**p
    a = Riemann(g,0,1,n)
    print(p,a)

# En gros, pour avoir p décimales exactes, il faut
# 10**(p+1) rectangles : trop coûteux en termes
# de nombre de calculs : 1 million de rectangles pour
#  obtenir 5 décimales !


## Q6a

# np.arange(a,b,h) :  comme range(a,b,h), mais
# a,b,h n'ont pas besoin d'être entiers !

def RG(f,a,b,h):
    S = 0
    for x in np.arange(a,b,h):
        S+= h *f(x)
    return S

# Remarque : suivant la formulation de l'énoncé (ex : calculer
# l'intégrale avec n=100 rectangles/calculer avec la méthode des
#rectangles avec un pas de h=10**-3), Riemann/RG sera plus
# adapté aux données de l'énoncé

## Q7
def RD(f,a,b,h):
    S = 0
    for x in np.arange(a+h,b,h):
        S+= h *f(x)
    return S

# Si a>b, comme h est ici >0, la liste np.arange est vide
# donc le calcul de la somme par RG renvoie 0 : il faudrait
# dans ce cas un pas h négatif.

# Adpatation de RG pour traiter le cas b<a :

def RG(f,a,b,h):
    pas =h
    if a>b:
        pas = -h
    S = 0
    for x in np.arange(a,b,pas):
        S+= pas*f(x)
    return S

# Rem : on doit faire la même adaptation avec RD

## Q5) Comparaison RD et  RG

for p in range(1,1+6):
    h = 10**(-p)
    gauche = RG(g,0,1,h) # sur-estime pi car g est décroissante
    droite = RD(g,0,1,h) # sous-estime
    print(f"pas = {h}  {droite} < pi < {gauche}")



## Q8. Fonction F

def f(t):
    return 1/np.sqrt(1+t**2+t**4)


def F(x):
    # Le nombre F(x) s'obtient par calcul d'une intégrale
    h=10**(-3) # choix du pas de la méthode complètement arbitraire
    return RG(f,x,2*x,h) # si x~1e-6, le segment  [x,2x] est de longueur 1e-6,
                         # donc le choix de h est très mauvais

# Tracé sur [-10,10] de la courbe de F

X = np.linspace(-10,10)
Y =[F(x) for x in X]     # liste en compréhension
plt.grid('on')
plt.plot(X,Y)
plt.show()

## Q9
# Le problème ici est que l'intégrande est une fonction de
# deux variables. Or, notre Riemann est programmé pour
# prendre en entrée une fonction d'une seule variable

# Solution -> Approche naïve : reprogrammer Riemann
#                              pour l'adapter à ce
#                              cas particulier

# SOLUTION RECOMMANDEE POUR LA PLUPART DES CANDIDATS
#

def f(t,x):
    return np.sin(t*x)/(x+t)

def F(x):
    n = 100
    h = 1/n
    S = 0
    for k in range(n):
        S+=h*f(h,x) # la fonction f est de deux variables en effet
    return S

# Approche plus conceptuelle (indique un bon niveau de maîtrise théorique)
def G(x):
    """
    G renvoie une fonction : autrement dit : G
    est une fonction, et G(x) est aussi une fonction
    """
    def g(t):
        return np.sin(t*x)/(x+t)
    return g          # g est une fonction. G(x) est bien une fonction

def F(x):
    Return Riemann (G(x),0,1,100)