Skip to content
Extraits de code Groupes Projets
Valider 1fcb7a82 rédigé par Jérôme de Favereau de Jeneret's avatar Jérôme de Favereau de Jeneret
Parcourir les fichiers

Update NM3_ex_solution.ipynb

parent e6a79d94
Aucune branche associée trouvée
Aucune étiquette associée trouvée
Aucune requête de fusion associée trouvée
Masquer les modifications d'espaces
En ligne Côte à côte
Affichage de
avec 1 ajout et 1 suppression
NM3_ex_solution.ipynb
+ 1
− 1
Voir le fichier @ 1fcb7a82
%% Cell type:markdown id:narrow-newspaper tags:
## Exercice : méthode de Romberg
Estimation de
$$
\int_0^{\frac{3\pi}{4}} \sin x\ dx
$$
Par la méthode de Romberg. Cette intégrale vaut $1 + 1/\sqrt{2} \approx 1.7071067811865475$
%% Cell type:code id:accomplished-coordinator tags:
``` python
import numpy as np
import math
def f(x):
return np.sin(x)
# méthode des trapèzes avec 2 points:
x_0 = 0
x_n = 3 * math.pi / 4
R11 = (f(x_0) + f(x_n))/2 * (x_n - x_0)
print(f"{R11=}")
```
%% Output
R11=0.8330405509046936
%% Cell type:markdown id:annoying-access tags:
Calcul de R21 par la méthode itérative des trapèzes:
%% Cell type:code id:fewer-partition tags:
``` python
R21 = R11/2 + f(x_n/2) * (x_n - x_0)/2
print(f"{R21=}")
```
%% Output
R21=1.5049402075046334
%% Cell type:markdown id:returning-mission tags:
calcul de R22 par Richardson avec $p=2$ puisque l'erreur est dominée par $h^2$:
calcul de R22 par Richardson avec $p=2$ puisque l'erreur sur R11 et R21 est dominée par $h^2$:
%% Cell type:code id:liable-forestry tags:
``` python
R22 = (4*R21 - R11)/3
print(f"{R22=}")
```
%% Output
R22=1.7289067597046133
%% Cell type:markdown id:streaming-tucson tags:
Calcul de R31 par la méthode itérative des trapèzes:
%% Cell type:code id:cathedral-investigation tags:
``` python
h = (x_n - x_0)/2 # espace entre les nouveaux points
x = x_0 + h/2 # position x du premier point
R31 = R21/2 + (f(x) + f(x+h)) * h/2
print(f"{R31=}")
```
%% Output
R31=1.6574582026781939
%% Cell type:markdown id:sacred-pension tags:
calcul de R32 par Richardson avec $p=2$ puisque l'erreur sur R21 et R31 est dominée par $h^2$:
%% Cell type:code id:compliant-zimbabwe tags:
``` python
R32 = (4*R31 - R21)/3
print(f"{R32=}")
```
%% Output
R32=1.708297534402714
%% Cell type:markdown id:elect-desktop tags:
calcul de R33 par Richardson avec $p=4$ puisque l'erreur sur R22 et R32 est dominée par $h^4$:
%% Cell type:code id:floppy-alignment tags:
``` python
R33 = (16*R32 - R22)/15
print(f"{R33=}")
```
%% Output
R33=1.706923586049254
%% Cell type:markdown id:surgical-syndrome tags:
Même chose en utilisant les fonctions définies au cours:
%% Cell type:code id:green-league tags:
``` python
import matplotlib.pyplot as plt
import numpy as np
from math import sin, pi, sqrt
def trapezes_rec(f, x_0, x_n, i_old, k):
if k == 1:
return (f(x_0) + f(x_n))*(x_n - x_0)/2
else:
n = 2**(k-2)
h = (x_n - x_0)/n
x = x_0 + h/2
sum = 0
for i in range(n):
sum += f(x)
x += h
return i_old/2 + h*sum/2
def richardson(i_1, i_2, k):
fact = 2**k
return (fact*i_2 - i_1)/(fact - 1)
# ligne 1
R11 = trapezes_rec(np.sin, 0, 3*pi/4, 0, 1)
print(f"{R11=}")
# ligne 2
R21 = trapezes_rec(np.sin, 0, 3*pi/4, R11, 2)
print(f"{R21=}")
R22 = richardson(R11, R21, 2)
print(f"{R22=}")
# ligne 3
R31 = trapezes_rec(np.sin, 0, 3*pi/4, R21, 3)
print(f"{R31=}")
R32 = richardson(R21, R31, 2)
print(f"{R32=}")
R33 = richardson(R22, R32, 4)
print(f"{R33=}")
# ligne 4
R41 = trapezes_rec(np.sin, 0, 3*pi/4, R31, 4)
print(f"{R41=}")
R42 = richardson(R31, R41, 2)
print(f"{R42=}")
R43 = richardson(R32, R42, 4)
print(f"{R43=}")
R44 = richardson(R33, R43, 6)
print(f"{R44=}")
print("\nEn utilisant la formule analytique:")
print('INT={}'.format(1+1/math.sqrt(2)))
```
%% Output
R11=0.8330405509046936
R21=1.5049402075046334
R22=1.7289067597046133
R31=1.6574582026781939
R32=1.708297534402714
R33=1.706923586049254
R41=1.6947487165588795
R42=1.7071788878524412
R43=1.7071043114157562
R44=1.7071071800723674
En utilisant la formule analytique:
INT=1.7071067811865475
%% Cell type:code id:leading-montana tags:
``` python
```
%% Cell type:code id:printable-neutral tags:
``` python
```
......
0% ou .
You are about to add 0 people to the discussion. Proceed with caution.
Terminez d'abord l'édition de ce message.
Veuillez vous inscrire ou vous pour commenter