import numpy as np
import matplotlib as mpl
mpl.rcParams['figure.figsize'] = (10, 7)
import matplotlib.pyplot as plt
from scipy import integrate
import numpy.random as rd
import seaborn as sns
sns.set(context="notebook", style="whitegrid", palette="hls", font="sans-serif", font_scale=1.1)
import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate
import numpy.random as rd
def I(n):
def f(t):
return 1 / ((1+t)**n * np.sqrt(1-t))
i, err = integrate.quad(f, 0, 1)
return i
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
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()
<Figure size 432x288 with 0 Axes>
[<matplotlib.lines.Line2D at 0x7fc29ae167f0>]
Text(0.5, 1.0, 'Valeurs de $I_n$')
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()
<Figure size 432x288 with 0 Axes>
[<matplotlib.lines.Line2D at 0x7fc29adbbd30>]
Text(0.5, 1.0, 'Valeurs de $\\ln(I_n)$ en fonction de $\\ln(n)$')
valeurs_n = np.arange(1, 500)
valeurs_In = np.array([I(n) for n in valeurs_n])
valeurs_Jn = np.array([J(n) for n in valeurs_n])
alpha = 0.9
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()
<Figure size 432x288 with 0 Axes>
[<matplotlib.lines.Line2D at 0x7fc29a86a358>]
[<matplotlib.lines.Line2D at 0x7fc29a8f7320>]
<matplotlib.legend.Legend at 0x7fc29ac1f400>
Text(0.5, 1.0, 'Valeurs de $n^{\\alpha} I_n$ et $n^{\\alpha} J_n$')
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()
<Figure size 720x504 with 0 Axes>
[<matplotlib.lines.Line2D at 0x7f1d46487358>]
<matplotlib.legend.Legend at 0x7f1d4656df98>
Text(0.5, 1.0, 'Valeurs de $n^{\\alpha} (I_n - J_n)$')
X = np.linspace(0, 100, 1000)
plt.plot(X, np.log(1 + X), 'ro-', label=r'$\log(1+x)$', markevery=50)
plt.plot(X, X / (1 + X), 'b+-', label=r'$\frac{x}{1+x}$', markevery=50)
plt.legend()
plt.title("Comparaison entre deux fonctions")
plt.show()
[<matplotlib.lines.Line2D at 0x7f1d462721d0>]
[<matplotlib.lines.Line2D at 0x7f1d462fcfd0>]
<matplotlib.legend.Legend at 0x7f1d46272710>
Text(0.5, 1.0, 'Comparaison entre deux fonctions')
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).
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 :
M = max_de_f = max(Ys)
print("Sur [0, 1], avec 2000 points, le maximum de f(x) est M =", M)
Sur [0, 1], avec 2000 points, le maximum de f(x) est M = 0.4812601808328304
i_maximisant_f = np.argmax(Ys)
x_maximisant_f = Xs[i_maximisant_f]
print("Sur ces 2000 points, le maximum est atteint en x =", x_maximisant_f)
Sur ces 2000 points, le maximum est atteint en x = 0.4062031015507754
On affiche la fonction, comme demandé, avec un titre :
plt.figure()
plt.plot(Xs, Ys)
plt.vlines(x_maximisant_f, min(Ys), max(Ys), color="b", linestyles="dotted")
plt.hlines(max(Ys), min(Xs), max(Xs), color="b", linestyles="dotted")
plt.title("Fonction $f(x)$ sur $[0,1]$")
plt.show()
<Figure size 720x504 with 0 Axes>
[<matplotlib.lines.Line2D at 0x7f1d45b94f60>]
<matplotlib.collections.LineCollection at 0x7f1d46061fd0>
<matplotlib.collections.LineCollection at 0x7f1d45b744e0>
Text(0.5, 1.0, 'Fonction $f(x)$ sur $[0,1]$')
Pour calculer l'intégrale, on utilise scipy.integrate.quad
:
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 :
for n in range(10):
print("In(x,", n, ") =", In(Xs, n))
In(x, 0 ) = 1.0 In(x, 1 ) = 0.16666666666666663 In(x, 2 ) = 0.049987680821294386 In(x, 3 ) = 0.017874498250810132 In(x, 4 ) = 0.0069477925896280456 In(x, 5 ) = 0.0028398023087289567 In(x, 6 ) = 0.0012011683591199485 In(x, 7 ) = 0.00052074827894252 In(x, 8 ) = 0.00022991112607879838 In(x, 9 ) = 0.00010290185158968158
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), markevery=20)
plt.legend()
plt.show()
<Figure size 720x504 with 0 Axes>
[<matplotlib.lines.Line2D at 0x7f1d458e57b8>]
[<matplotlib.lines.Line2D at 0x7f1d458e5b70>]
[<matplotlib.lines.Line2D at 0x7f1d458e5eb8>]
[<matplotlib.lines.Line2D at 0x7f1d458ee278>]
[<matplotlib.lines.Line2D at 0x7f1d458ee5f8>]
<matplotlib.legend.Legend at 0x7f1d45944860>
$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$.
def un(n):
return In(Xs, n + 1) / In(Xs, n)
for n in range(10):
print("un =", un(n), "pour n =", n)
un = 0.16666666666666663 pour n = 0 un = 0.2999260849277664 pour n = 1 un = 0.35757806637822104 pour n = 2 un = 0.38869860804697826 pour n = 3 un = 0.4087344681198935 pour n = 4 un = 0.4229760485185215 pour n = 5 un = 0.43353479550864377 pour n = 6 un = 0.44150146121592054 pour n = 7 un = 0.4475723004132197 pour n = 8 un = 0.4522404636909894 pour n = 9
Ici, un
ne peut pas être utilisé comme une fonction "numpy" qui travaille sur un tableau, on stocke donc les valeurs "plus manuellement" :
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(True)
plt.show()
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 ?
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.
case_max = 12
univers = list(range(case_max))
def prochaine_case(case):
return (case + rd.randint(1, 6+1)) % case_max
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 :
[Yn(1) for _ in range(10)]
[4, 2, 5, 2, 4, 6, 3, 6, 6, 1]
[Yn(100) for _ in range(10)]
[9, 7, 9, 2, 4, 5, 3, 2, 5, 5]
Pour l'histogramme, on triche un peu en utilisant numpy.bincount
. Mais on peut le faire à la main très simplement !
observations = [Yn(100) for _ in range(10)]
print(observations)
np.bincount(observations, minlength=case_max)
[4, 9, 3, 10, 0, 7, 1, 0, 8, 10]
array([2, 1, 0, 1, 1, 0, 0, 1, 1, 1, 2, 0])
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)
histogramme(50)
array([0.0902, 0.085 , 0.0832, 0.0838, 0.083 , 0.0822, 0.082 , 0.0894, 0.0808, 0.0728, 0.0792, 0.0884])
On va afficher des histogrammes :
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()
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.
case_max = 12
P = np.zeros((case_max, case_max))
for k in range(case_max):
for i in range(k - 6, k):
P[k, i] = 1/6
P
array([[0. , 0. , 0. , 0. , 0. , 0. , 0.16666667, 0.16666667, 0.16666667, 0.16666667, 0.16666667, 0.16666667], [0.16666667, 0. , 0. , 0. , 0. , 0. , 0. , 0.16666667, 0.16666667, 0.16666667, 0.16666667, 0.16666667], [0.16666667, 0.16666667, 0. , 0. , 0. , 0. , 0. , 0. , 0.16666667, 0.16666667, 0.16666667, 0.16666667], [0.16666667, 0.16666667, 0.16666667, 0. , 0. , 0. , 0. , 0. , 0. , 0.16666667, 0.16666667, 0.16666667], [0.16666667, 0.16666667, 0.16666667, 0.16666667, 0. , 0. , 0. , 0. , 0. , 0. , 0.16666667, 0.16666667], [0.16666667, 0.16666667, 0.16666667, 0.16666667, 0.16666667, 0. , 0. , 0. , 0. , 0. , 0. , 0.16666667], [0.16666667, 0.16666667, 0.16666667, 0.16666667, 0.16666667, 0.16666667, 0. , 0. , 0. , 0. , 0. , 0. ], [0. , 0.16666667, 0.16666667, 0.16666667, 0.16666667, 0.16666667, 0.16666667, 0. , 0. , 0. , 0. , 0. ], [0. , 0. , 0.16666667, 0.16666667, 0.16666667, 0.16666667, 0.16666667, 0.16666667, 0. , 0. , 0. , 0. ], [0. , 0. , 0. , 0.16666667, 0.16666667, 0.16666667, 0.16666667, 0.16666667, 0.16666667, 0. , 0. , 0. ], [0. , 0. , 0. , 0. , 0.16666667, 0.16666667, 0.16666667, 0.16666667, 0.16666667, 0.16666667, 0. , 0. ], [0. , 0. , 0. , 0. , 0. , 0.16666667, 0.16666667, 0.16666667, 0.16666667, 0.16666667, 0.16666667, 0. ]])
import numpy.linalg as LA
spectre, vecteur_propres = LA.eig(P)
On a besoin d'éliminer les erreurs d'arrondis, mais on voit que $1$ est valeur propre, associée au vecteur $[1,\dots,1,\dots,1]$ avec un $1$ seulement à la 8ème composante.
np.round(spectre,10)
array([ 1. +0.j , -0.16666667+0.62200847j, -0.16666667-0.62200847j, -0.16666667+0.16666667j, -0.16666667-0.16666667j, -0.16666667+0.0446582j , -0.16666667-0.0446582j , 0. +0.j , 0. +0.j , -0. +0.j , 0. +0.j , 0. -0.j ])
np.round(vecteur_propres[0])
array([ 0.+0.j, 0.+0.j, 0.-0.j, -0.+0.j, -0.-0.j, -0.+0.j, -0.-0.j, 0.+0.j, 0.+0.j, 0.+0.j, -0.+0.j, -0.-0.j])
$P$ n'est pas diagonalisable, à prouver au tableau si l'examinateur le demande.
def f(x):
return 1 / (2 - np.exp(x))
Donc en $x=0$, on utilise que $g^{(k)}(x) = - \exp(x)$, qui donne que $g^{(k)}(0) = 1$ si $k=0$ ou $-1$ sinon, pour trouver que $\sum_{k=0}^n {n \choose k} f^{(k)}(0) = f^{(n)}(0)$. En écrivant ${n \choose k} = \frac{k! (n-k)!}{n!}$ et avec la formule définissant $a(n)$, cela donne directement la somme recherchée : $$ a(n) = \sum_{k=1}^n \frac{a(n-k)}{k!}.$$
math.factorial
pour calculer $k!$.Il nous faut aussi $a(0) = f(0) = 1$ :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
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])
Pour n = 0 on a a(n) = 1.0 Pour n = 1 on a a(n) = 1.0 Pour n = 2 on a a(n) = 1.5 Pour n = 3 on a a(n) = 2.1666666666666665 Pour n = 4 on a a(n) = 3.1249999999999996 Pour n = 5 on a a(n) = 4.508333333333334 Pour n = 6 on a a(n) = 6.504166666666667 Pour n = 7 on a a(n) = 9.383531746031748 Pour n = 8 on a a(n) = 13.537574404761907 Pour n = 9 on a a(n) = 19.530591380070547 Pour n = 10 on a a(n) = 28.176687334656087
plt.figure()
plt.plot(valeurs_n, valeurs_a, 'ro-', label=r'$a(n)$', markersize=12)
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()
<Figure size 720x504 with 0 Axes>
[<matplotlib.lines.Line2D at 0x7f1d45585588>]
[<matplotlib.lines.Line2D at 0x7f1d460c6a90>]
[<matplotlib.lines.Line2D at 0x7f1d4556e5f8>]
Text(0.5, 1.0, '$a(n)$ et deux autres suites')
<matplotlib.legend.Legend at 0x7f1d4556eda0>
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 :
x = 0.5
for n in range(0, 6 + 1):
print("Pour n =", n, "S_n(x) =", Sn(x, n))
Pour n = 0 S_n(x) = 1.0 Pour n = 1 S_n(x) = 1.5 Pour n = 2 S_n(x) = 1.875 Pour n = 3 S_n(x) = 2.1458333333333335 Pour n = 4 S_n(x) = 2.3411458333333335 Pour n = 5 S_n(x) = 2.4820312500000004 Pour n = 6 S_n(x) = 2.583658854166667
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]$ !
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, 'o:', label='$S_' + str(n) + '(x)$', markevery=50)
plt.plot(valeurs_x, valeurs_f, '+-', label='$f(x)$', markevery=50)
plt.title("$f(x)$ et $S_n(x)$ pour $n = 0$ à $n = 6$")
plt.legend()
plt.show()
<Figure size 720x504 with 0 Axes>
[<matplotlib.lines.Line2D at 0x7f1d455ca748>]
[<matplotlib.lines.Line2D at 0x7f1d4536def0>]
[<matplotlib.lines.Line2D at 0x7f1d455cae48>]
[<matplotlib.lines.Line2D at 0x7f1d455caeb8>]
[<matplotlib.lines.Line2D at 0x7f1d455b8550>]
[<matplotlib.lines.Line2D at 0x7f1d455b8828>]
[<matplotlib.lines.Line2D at 0x7f1d453bf898>]
[<matplotlib.lines.Line2D at 0x7f1d4536d860>]
Text(0.5, 1.0, '$f(x)$ et $S_n(x)$ pour $n = 0$ à $n = 6$')
<matplotlib.legend.Legend at 0x7f1d455b8fd0>
def u(n):
return np.arctan(n+1) - np.arctan(n)
valeurs_n = np.arange(50)
valeurs_u = u(valeurs_n)
plt.figure()
plt.plot(valeurs_n, valeurs_u, "o-")
plt.xlabel("Entier $n$")
plt.ylabel("Valeur $u_n$")
plt.title("Premières valeurs de $u_n$")
<Figure size 720x504 with 0 Axes>
[<matplotlib.lines.Line2D at 0x7f1d457a2a90>]
Text(0.5, 0, 'Entier $n$')
Text(0, 0.5, 'Valeur $u_n$')
Text(0.5, 1.0, 'Premières valeurs de $u_n$')
On peut vérifier le prognostic quand à la somme de la série $\sum u_n$ :
pi/2
1.5707963267948966
sum(valeurs_u)
1.550798992821746
somme_serie = pi/2
somme_partielle = sum(valeurs_u)
erreur_relative = abs(somme_partielle - somme_serie) / somme_serie
erreur_relative
0.01273069820194558
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 :
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.xlabel("Entier $n$")
plt.title(r"Valeurs de $u_n / \frac{1}{n(n+1)}$")
<Figure size 720x504 with 0 Axes>
[<matplotlib.lines.Line2D at 0x7f1d453292b0>]
Text(0.5, 0, 'Entier $n$')
Text(0.5, 1.0, 'Valeurs de $u_n / \\frac{1}{n(n+1)}$')
from math import ceil, sqrt, pi
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
def e010101(n):
return 1 if n % 2 == 0 else 0
delta = 1e-5
Se010101 = Se(e010101, delta)
print("Pour delta =", delta, "on a Se010101(delta) ~=", round(Se010101, 5))
Pour delta = 1e-05 on a Se010101(delta) ~= 1.06408
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.
from random import random
def pile(proba):
""" True si pile, False si face (false, face, facile à retenir)."""
return random() < proba
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
import numpy as np
lances = [ En(2, 0.5) for _ in range(100) ]
np.bincount(lances)
array([30, 70])
def pn(n, p, nbSimulations=100000):
return np.mean([ En(n, p) for _ in range(nbSimulations) ])
pn(2, 0.5)
0.75054
pn(4, 0.5)
0.56445
pn(4, 0.1)
0.97283
pn(4, 0.9)
0.10047
pn(6, 0.2)
0.86535
pn(20, 0.2)
0.80335
pn(100, 0.2)
0.80053
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).
from math import floor, log, pi
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
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))
Pour x = -0.75 f(x) = -0.64276 Pour x = -0.5 f(x) = -0.44841 Pour x = 0.25 f(x) = 0.26765 Pour x = 0 f(x) = 0 Pour x = 0.25 f(x) = 0.26765 Pour x = 0.5 f(x) = 0.58224 Pour x = 0.75 f(x) = 0.97847
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 :
from scipy import integrate
def g(x):
def h(t):
return log(1 - t) / t
integrale, erreur = integrate.quad(h, 0, x)
return integrale
import numpy as np
import matplotlib.pyplot as plt
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)$", markevery=50)
plt.plot(domaine, valeurs_g, "+-", label="$g(x)$", markevery=50)
plt.legend()
plt.grid()
plt.title("Représentation de $f(x)$ et $g(x)$")
plt.show()
<Figure size 720x504 with 0 Axes>
[<matplotlib.lines.Line2D at 0x7f1d442255f8>]
[<matplotlib.lines.Line2D at 0x7f1d442c9e80>]
<matplotlib.legend.Legend at 0x7f1d44225780>
Text(0.5, 1.0, 'Représentation de $f(x)$ et $g(x)$')
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 !
Ce document est distribué sous licence libre (MIT), comme les autres notebooks que j'ai écrit depuis 2015.