In [1]:
import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate
import numpy.random as rd
In [2]:
def I(n):
def f(t):
return 1 / ((1+t)**n * np.sqrt(1-t))
i, err = integrate.quad(f, 0, 1)
return i
In [3]:
def J(n):
def f(t):
return 1 / ((1+t)**n * np.sqrt(1-t))
i, err = integrate.quad(f, 0, 0.5)
return i
In [12]:
valeurs_n = np.arange(1, 50)
valeurs_In = np.array([I(n) for n in valeurs_n])
plt.figure()
plt.plot(valeurs_n, valeurs_In, 'ro')
plt.title("Valeurs de $I_n$")
plt.show()
Out[12]:
Out[12]:
Out[12]:
In [13]:
plt.figure()
plt.plot(np.log(valeurs_n), np.log(valeurs_In), 'go')
plt.title(r"Valeurs de $\ln(I_n)$ en fonction de $\ln(n)$")
plt.show()
Out[13]:
Out[13]:
Out[13]:
In [14]:
valeurs_Jn = np.array([J(n) for n in valeurs_n])
alpha = 1
plt.figure()
plt.plot(valeurs_n, valeurs_n**alpha * valeurs_In, 'r+', label=r'$n^{\alpha} I_n$')
plt.plot(valeurs_n, valeurs_n**alpha * valeurs_Jn, 'b+', label=r'$n^{\alpha} J_n$')
plt.legend()
plt.title(r"Valeurs de $n^{\alpha} I_n$ et $n^{\alpha} J_n$")
plt.show()
Out[14]:
Out[14]:
Out[14]:
Out[14]:
Out[14]:
In [16]:
plt.figure()
plt.plot(valeurs_n, valeurs_n**alpha * (valeurs_In - valeurs_Jn), 'g+', label=r'$n^{\alpha} (I_n - J_n)$')
plt.legend()
plt.title(r"Valeurs de $n^{\alpha} (I_n - J_n)$")
plt.show()
Out[16]:
Out[16]:
Out[16]:
Out[16]:
In [18]:
X = np.linspace(0, 100, 10000)
plt.plot(X, np.log(1 + X), 'r-', label=r'$\log(1+x)$')
plt.plot(X, X / (1 + X), 'b-', label=r'$\frac{x}{1+x}$')
plt.legend()
plt.title("Comparaison entre deux fonctions")
plt.show()
Out[18]:
Out[18]:
Out[18]:
Out[18]:
On commence par définir la fonction, en utilisant numpy.cos et pas math.cos (les fonctions de numpy peuvent travailler sur des tableaux, c'est plus pratique).
In [ ]:
def f(x):
return x * (1 - x) * (1 + np.cos(5 * np.pi * x))
Xs = np.linspace(0, 1, 2000)
Ys = f(Xs)
Pas besoin de lire le maximum sur un graphique :
In [3]:
M = max_de_f = max(Ys)
print("Sur [0, 1], avec 2000 points, M =", M)
On affiche la fonction, comme demandé, avec un titre :
In [4]:
plt.figure()
plt.plot(Xs, Ys)
plt.title("Fonction $f(x)$ sur $[0,1]$")
plt.show()
Out[4]:
Out[4]:
Out[4]:
Pour calculer l'intégrale, on utilise scipy.integrate.quad :
In [11]:
def In(x, n):
def fn(x):
return f(x) ** n
return integrate.quad(fn, 0, 1)[0]
def Sn(x):
return np.sum([In(Xs, n) * x**n for n in range(0, n+1)], axis=0)
On vérifie avant de se lancer dans l'affichage :
In [12]:
for n in range(10):
print("In(x,", n, ") =", In(Xs, n))
In [13]:
a = 1/M + 0.1
X2s = np.linspace(-a, a, 2000)
plt.figure()
for n in [10, 20, 30, 40, 50]:
plt.plot(X2s, Sn(X2s), label="n =" + str(n))
plt.legend()
plt.show()
Out[13]:
Out[13]:
Out[13]:
Out[13]:
Out[13]:
Out[13]:
Out[13]:
$S_n(x)$ semble diverger pour $x\to2^-$ quand $n\to\infty$. Le rayon de convergence de la série $\sum In x^n$ semble être $2$.
In [14]:
def un(n):
return In(Xs, n + 1) / In(Xs, n)
In [15]:
for n in range(10):
print("un =", un(n), "pour n =", n)
Ici, un ne peut pas être utilisé comme une fonction "numpy" qui travaille sur un tableau, on stocke donc les valeurs "plus manuellement" :
In [26]:
def affiche_termes_un(N):
valeurs_un = [0] * N
for n in range(N):
valeurs_un[n] = un(n)
plt.figure()
plt.plot(valeurs_un, 'o-')
plt.title("Suite $u_n$")
plt.grid()
plt.show()
In [27]:
affiche_termes_un(30)
La suite $u_n$ semble être croissante (on peut le prouver), toujours plus petite que $1$ (se prouve facilement aussi, $I_{n+1} < I_n$), et semble converger. Peut-être vers $1/2$, il faut aller regarder plus loin ?
In [28]:
affiche_termes_un(100)
Pour conclure, on peut prouver que la suite est monotone et bornée, donc elle converge. Il est plus dur de calculer sa limite, et cela sort de l'exercice.
In [29]:
case_max = 12
univers = list(range(case_max))
In [30]:
def prochaine_case(case):
return (case + rd.randint(1, 6+1)) % case_max
In [48]:
def Yn(duree, depart=0):
case = depart
for coup in range(duree):
case = prochaine_case(case)
return case
Avant de s'en servir pour simuler plein de trajectoirs, on peut vérifier :
In [32]:
[Yn(1) for _ in range(10)]
Out[32]:
In [39]:
[Yn(100) for _ in range(10)]
Out[39]:
Pour l'histogramme, on triche un peu en utilisant numpy.bincount. Mais on peut le faire à la main très simplement !
In [40]:
np.bincount(_, minlength=case_max)
Out[40]:
In [44]:
def histogramme(duree, repetitions=5000):
cases = [Yn(duree) for _ in range(repetitions)]
frequences = np.bincount(cases, minlength=case_max)
# aussi a la main si besoin
frequences = [0] * case_max
for case in cases:
frequences[case] += 1
return frequences / np.sum(frequences)
In [45]:
histogramme(50)
Out[45]:
In [55]:
def voir_histogramme(valeurs_n):
for n in valeurs_n:
plt.figure()
plt.bar(np.arange(case_max), histogramme(n))
plt.title("Histogramme de cases visitées en " + str(n) + " coups")
plt.show()
In [58]:
voir_histogramme([1, 2, 3, 50, 100, 200])
On s'approche d'une distribution uniforme !
On a tout simplement l'expression suivante : $$\forall n \geq 0, \mathbb{P}(Y_{n+1} = k) = \frac{1}{6} \sum_{\delta = 1}^{6} \mathbb{P}(Y_n = k - \delta \mod 12).$$ Avec $k - 1 \mod 12 = 11$ si $k = 0$ par exemple.
On a donc la matrice suivante pour exprimer $U_n = (\mathbb{P}(Y_n = k))_{0\leq k \leq 11}$ en fonction de $U_{n-1}$ :
$$ P = \frac{1}{6} \begin{bmatrix} 0 & 0 & 0 & 0 & 0 & 0 & 1 & 1 & 1 & 1 & 1 & 1 \\ 1 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 1 & 1 & 1 & 1\\ 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 1 & 1 & 1\\ 1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 1 & 1\\ 1 & 1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 1\\ 1 & 1 & 1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 1\\ 1 & 1 & 1 & 1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 1 & 1 & 1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 1 & 1 & 1 & 1 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 1 & 1 & 1 & 1 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 1 & 1 & 1 & 1 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & 1 & 1 & 1 & 1 & 1 & 0 \\ \end{bmatrix}$$On va la définir rapidement en Python, et calculer ses valeurs propres notamment.
In [70]:
P = np.zeros((case_max, case_max))
In [71]:
for k in range(case_max):
for i in range(k - 6, k):
P[k, i] = 1
In [72]:
P
Out[72]:
In [73]:
import numpy.linalg as LA
In [75]:
spectre, vecteur_propres = LA.eig(P)
On a besoin d'éliminer les erreurs d'arrondis, mais on voit que $6$ est valeur propre, associée au vecteur $[0,\dots,1,\dots,0]$ avec un $1$ seulement à la 8ème composante.
In [84]:
np.round(spectre)
Out[84]:
In [83]:
np.round(vecteur_propres[0])
Out[83]:
$P$ n'est pas diagonalisable, à prouver au tableau si l'examinateur le demande.
In [33]:
def f(x):
return 1 / (2 - np.exp(x))
math.factorial pour calculer $k!$. Il nous faut aussi $a(0) = f(0) = 1$ :
In [22]:
from math import factorial
def a_0an(nMax):
valeurs_a = np.zeros(nMax+1)
valeurs_a[0] = 1.0
for n in range(1, nMax+1):
valeurs_a[n] = sum(valeurs_a[n-k] / factorial(k) for k in range(1, n+1))
return valeurs_a
In [26]:
nMax = 10
valeurs_n = np.arange(0, nMax + 1)
valeurs_a = a_0an(nMax)
for n in valeurs_n:
print("Pour n =", n, "on a a(n) =", valeurs_a[n])
In [29]:
plt.figure()
plt.plot(valeurs_n, valeurs_a, 'ro', label=r'$a(n)$')
plt.plot(valeurs_n, 1 / np.log(2)**valeurs_n, 'g+', label=r'$1/\log(2)^n$')
plt.plot(valeurs_n, 1 / (2 * np.log(2)**valeurs_n), 'bd', label=r'$1/(2\log(2)^n)$')
plt.title("$a(n)$ et deux autres suites")
plt.legend()
plt.show()
Out[29]:
Out[29]:
Out[29]:
Out[29]:
Out[29]:
Out[29]:
In [30]:
def Sn(x, n):
valeurs_a = a_0an(n)
return sum(valeurs_a[k] * x**k for k in range(0, n + 1))
On peut vérifie que notre fonction marche :
In [32]:
x = 0.5
for n in range(0, 6 + 1):
print("Pour n =", n, "S_n(x) =", Sn(x, n))
In [68]:
valeurs_x = np.linspace(0, 0.5, 1000)
valeurs_f = f(valeurs_x)
Je pense que l'énoncé comporte une typo sur l'intervale ! Vu le rayon de convergence, on ne voit rien si on affiche sur $[0,10]$ !
In [67]:
plt.figure()
for n in range(0, 6 + 1):
valeurs_Sn = []
for x in valeurs_x:
valeurs_Sn.append(Sn(x, n))
plt.plot(valeurs_x, valeurs_Sn, ':', label='$S_' + str(n) + '(x)$')
plt.plot(valeurs_x, valeurs_f, '-', label='$f(x)$')
plt.title("$f(x)$ et $S_n(x)$ pour $n = 0$ à $n = 6$")
plt.legend()
plt.show()
Out[67]:
Out[67]:
Out[67]:
Out[67]:
Out[67]:
Out[67]:
Out[67]:
Out[67]:
Out[67]:
Out[67]:
Out[67]:
In [2]:
def u(n):
return np.arctan(n+1) - np.arctan(n)
In [42]:
valeurs_n = np.arange(50)
valeurs_u = u(valeurs_n)
plt.figure()
plt.plot(valeurs_n, valeurs_u, "o-")
plt.title("Premières valeurs de $u_n$")
Out[42]:
Out[42]:
Out[42]:
On peut vérifier le prognostic quand à la somme de la série $\sum u_n$ :
In [44]:
pi/2
Out[44]:
In [43]:
sum(valeurs_u)
Out[43]:
In [45]:
somme_serie = pi/2
somme_partielle = sum(valeurs_u)
erreur_relative = abs(somme_partielle - somme_serie) / somme_serie
erreur_relative
Out[45]:
Avec seulement $50$ termes, on a moins de $1.5%$ d'erreur relative, c'est déjà pas mal !
$(u_n)_n$ semble être décroisante, et tendre vers $0$. On peut prouver ça mathématiquement.
On sait aussi que $\forall x\neq0, \arctan(x) + \arctan(1/x) = \frac{\pi}{2}$, et que $\arctan(x) \sim x$, donc on obtient que $u_n \sim \frac{1}{n} - \frac{1}{n+1} = \frac{1}{n(n+1)}$. On peut le vérifier :
In [4]:
valeurs_n = np.arange(10, 1000)
valeurs_u = u(valeurs_n)
valeurs_equivalents = 1 / (valeurs_n * (valeurs_n + 1))
plt.figure()
plt.plot(valeurs_n, valeurs_u / valeurs_equivalents, "-")
plt.title(r"Valeurs de $u_n / \frac{1}{n(n+1)}$")
Out[4]:
Out[4]:
Out[4]:
In [36]:
from math import ceil, sqrt, pi
In [37]:
def Se(e, delta=1e-5, borne_sur_n_0=10000):
borne_sur_n_1 = int(ceil(1 + sqrt(delta)/2.0))
borne_sur_n = max(borne_sur_n_0, borne_sur_n_1)
somme_partielle = 0
for n in range(0, borne_sur_n + 1):
somme_partielle += e(n) * u(n)
return somme_partielle
In [38]:
def e010101(n):
return 1 if n % 2 == 0 else 0
In [39]:
delta = 1e-5
Se010101 = Se(e010101, delta)
print("Pour delta =", delta, "on a Se010101(delta) ~=", round(Se010101, 5))
In [50]:
def inverse_Se(x, n):
assert 0 < x < pi/2.0, "Erreur : x doit être entre 0 et pi/2 strictement."
print("Je vous laisse chercher.")
raise NotImplementedError
Ca suffit pour la partie Python.
In [51]:
from random import random
def pile(proba):
""" True si pile, False si face (false, face, facile à retenir)."""
return random() < proba
In [52]:
def En(n, p):
lance = pile(p)
for i in range(n - 1):
nouveau_lance = pile(p)
if lance and nouveau_lance:
return False
nouveau_lance = lance
return True
In [55]:
import numpy as np
In [56]:
lances = [ En(2, 0.5) for _ in range(100) ]
np.bincount(lances)
Out[56]:
In [59]:
def pn(n, p, nbSimulations=100000):
return np.mean([ En(n, p) for _ in range(nbSimulations) ])
In [60]:
pn(2, 0.5)
Out[60]:
In [61]:
pn(4, 0.5)
Out[61]:
In [62]:
pn(4, 0.1)
Out[62]:
In [63]:
pn(4, 0.9)
Out[63]:
In [66]:
pn(6, 0.2)
Out[66]:
In [67]:
pn(20, 0.2)
Out[67]:
In [69]:
pn(100, 0.2)
Out[69]:
Le domaine de définition de $f(x) = \sum_{n \geq 1} \frac{x^n}{n^2}$ est $[-1, 1]$ car $\sum \frac{x^n}{n^k}$ converge si $\sum x^n$ converge (par $k$ dérivations successives), qui converge ssi $|x| < 1$. Et en $-1$ et $1$, on utilise $\sum \frac{1}{n^2} = \frac{\pi^2}{6}$.
Pour calculer $f(x)$ à $10^{-5}$ près, il faut calculer sa somme partielle $S_n(x) := \sum_{i=1}^n \frac{x^i}{i^2}$ en bornant son reste $S_n(x) := \sum_{i \geq n+1} \frac{x^i}{i^2}$ par (au moins) $10^{-5}$. Une inégalité montre rapidement que $R_n(x) \leq |x|^{n+1}\sum_{i\geq n+1} \frac{1}{i^2} $, et donc $R_n(x) \leq \delta$ dès que $|x|^{n+1} \leq \frac{\pi^2}{6} \delta$, puisque $\sum_{i\geq n+1} \frac{1}{i^2} \leq \sum_{i=0}^{+\infty} \frac{1}{i^2} = \frac{\pi^2}{6}$. En inversant pour trouver $n$, cela donne que le reste est contrôlé par $\delta$ dès que $n \leq \log_{|x|}\left( \frac{6}{\pi^2} \delta \right) - 1$ (si $x\neq 0$, et par $n \geq 0$ sinon).
In [72]:
from math import floor, log, pi
In [73]:
delta = 1e-5
def f(x):
if x == 0: return 0
borne_sur_n = int(floor(log((6/pi**2 * delta), abs(x)) - 1))
somme_partielle = 0
for n in range(1, borne_sur_n + 1):
somme_partielle += x**n / n**2
return somme_partielle
In [76]:
for x in [-0.75, -0.5, 0.25, 0, 0.25, 0.5, 0.75]:
print("Pour x =", x, "\tf(x) =", round(f(x), 5))
L'intégrale $g(t) = \int_0^x \frac{\ln(1 - t)}{t} \mathrm{d}t$ est bien défine sur $D = [-1, 1]$ puisque son intégrande existe, est continue et bien intégrable sur tout interval de la forme $]a, 0[$ ou $]0, b[$ pour $-1 < a < 0$ ou $0 < b < 1$. Le seul point qui peut déranger l'intégrabilité est en $0$, mais $\ln(1-t) \sim t$ quand $t\to0$ donc l'intégrande est $\sim 1$ en $0^-$ et $0^+$ et donc est bien intégrable. De plus, comme "intégrale de la borne supérieure" d'une fonction continue, $g$ est dérivable sur l'intérieur de son domaine, i.e., sur $]-1, 1[$.
Pour la calculer numériquement, on utilise évidemment le module scipy.integrate et sa fonction integrale, erreur = quad(f, a, b), qui donne une approximation de la valeur d'une intégrale en dimension 1 et une borne sur son erreur :
In [77]:
from scipy import integrate
In [78]:
def g(x):
def h(t):
return log(1 - t) / t
integrale, erreur = integrate.quad(h, 0, x)
return integrale
In [79]:
import numpy as np
import matplotlib.pyplot as plt
In [80]:
domaine = np.linspace(-0.99, 0.99, 1000)
valeurs_f = [f(x) for x in domaine]
valeurs_g = [g(x) for x in domaine]
plt.figure()
plt.plot(domaine, valeurs_f, label="$f(x)$")
plt.plot(domaine, valeurs_g, label="$g(x)$")
plt.legend()
plt.grid()
plt.title("Représentation de $f(x)$ et $g(x)$")
plt.show()
Out[80]:
Out[80]:
Out[80]:
Out[80]:
Out[80]:
La suite des questions est à faire au brouillon et sans Python :
On trouve que $f'(x) = \sum_{n\geq 1} \frac{n x^{n-1}}{n^2} = \frac{1}{x} \sum_{n\geq 1} \frac{x^n}{n}$ si $x\neq0$. Or on sait que $\log(1 + x) = \sum_{n\geq 1} \frac{x^n}{n}$ et donc cela montre bien que $g(x) = \int_0^x - f'(t) \mathrm{d}t = f(0) - f(x) = f(x)$ comme observé.
On trouve que $g(1) = - f(1) = - \frac{\pi^2}{6}$.
Par ailleurs, un changement de variable $u=1-x$ donne $g(1-x) = \int_x^1 \frac{\ln(u)}{1-u} \mathrm{d} u$, et une intégration par partie avec $a(u) = \ln(u)$ et $b'(u) = \frac{1}{1-u}$ donne $g(1-x) = [\ln(u)\ln(1-u)]_x^1 + \int_x^1 \frac{\ln(1-u)}{u} \mathrm{d}u$ et donc on reconnaît que $$g(1-x) = \ln(x)\ln(1-x) + g(1) - g(x).$$
Je vous laisse la fin comme exercice !
Se préparer aux oraux de "maths avec Python" (maths 2) du concours Centrale Supélec peut être utile.
Après les écrits et la fin de l'année, pour ceux qui seront admissibles à Centrale-Supélec, ils vous restera les oraux (le concours Centrale-Supélec a un oral d'informatique, et un peu d'algorithmique et de Python peuvent en théorie être demandés à chaque oral de maths et de SI).
Je vous invite à lire cette page avec attention, et à jeter un œil aux documents mis à disposition :
Pour réviser : voir ce tutoriel Matplotlib (en anglais), ce tutoriel Numpy (en anglais). Ainsi que tous les TP, TD et DS en Python que j'ai donné et corrigé au Lycée Lakanal (Sceaux, 92) en 2015-2016 !
Ces 5 sujets sont corrigés, et nous les avons tous traité en classe durant les deux TP de révisions pour les oraux (10 et 11 juin).
Ce document est distribué sous licence libre (MIT), comme les autres notebooks que j'ai écrit depuis 2015.