🎯 But pédagogique¶
Étudier l’influence de la discrétisation spatiale sur la précision du calcul de la variance de blocs à partir d’un modèle de covariance.
⚙️ Impact de la discrétisation¶
Dans la pratique numérique, cette intégrale double est approchée par une discrétisation spatiale du bloc en un nombre fini de points. La précision du calcul de la variance de bloc dépend donc de la résolution choisie :
- Une faible résolution (peu de points) entraîne une approximation grossière et une estimation moins précise de la variance.
- Une résolution élevée (beaucoup de points) améliore la précision mais augmente le coût de calcul.
Ce scénario permet d’étudier quantitativement cet impact, en observant comment la variance estimée varie avec la densité de points de discrétisation.
Source
import numpy as np
import matplotlib.pyplot as plt
import ipywidgets as widgets
from IPython.display import display, clear_output
from numpy.polynomial.legendre import leggauss
# --- Modèles de covariance ---
def spherical_covariance(h, range_, sill):
cov = np.zeros_like(h)
mask = h <= range_
hr = h[mask] / range_
cov[mask] = sill * (1 - 1.5 * hr + 0.5 * hr**3)
return cov
def exponential_covariance(h, range_, sill):
return sill * np.exp(-3 * h / range_)
def gaussian_covariance(h, range_, sill):
return sill * np.exp(-3 * (h / range_)**2)
def get_covariance_model(name):
return {'sphérique': spherical_covariance,
'exponentiel': exponential_covariance,
'gaussien': gaussian_covariance}[name]
# --- Quadrature de Gauss avec coordonnées des points ---
def block_variance_gauss_quad(geometry, lx, ly, lz, sill,
range_x, range_y, range_z,
model_name, n_points=5):
cov_func = get_covariance_model(model_name)
pts_1D, w_1D = leggauss(n_points)
pts_1D = 0.5 * (pts_1D + 1)
w_1D = 0.5 * w_1D
if geometry == 'ligne':
points = pts_1D * lx
weights = w_1D
diff = (points[:, None] - points[None, :]) / range_x
h = np.abs(diff)
cov_mat = cov_func(h, range_=1.0, sill=sill)
variance = np.sum(weights[:, None] * weights[None, :] * cov_mat)
return variance, points, np.zeros_like(points), None
elif geometry == 'surface':
X, Y = np.meshgrid(pts_1D * lx, pts_1D * ly)
points = np.column_stack([X.ravel(), Y.ravel()])
weights = np.outer(w_1D, w_1D).ravel()
diff = (points[:, None, :] - points[None, :, :]) / np.array([range_x, range_y])
h = np.linalg.norm(diff, axis=2)
cov_mat = cov_func(h, range_=1.0, sill=sill)
variance = np.sum(weights[:, None] * weights[None, :] * cov_mat)
return variance, points[:, 0], points[:, 1], None
elif geometry == 'cube':
X, Y, Z = np.meshgrid(pts_1D * lx, pts_1D * ly, pts_1D * lz)
points = np.column_stack([X.ravel(), Y.ravel(), Z.ravel()])
weights = np.outer(np.outer(w_1D, w_1D), w_1D).ravel()
diff = (points[:, None, :] - points[None, :, :]) / np.array([range_x, range_y, range_z])
h = np.linalg.norm(diff, axis=2)
cov_mat = cov_func(h, range_=1.0, sill=sill)
variance = np.sum(weights[:, None] * weights[None, :] * cov_mat)
return variance, points[:, 0], points[:, 1], points[:, 2]
return sill, None, None, None
# --- Widgets UI ---
style_desc = {'description_width': 'initial'}
# Section 1 : Anisotropie
range_x_input = widgets.FloatSlider(value=30, min=1, max=100, step=1, description='Portée X ($a_x$)', style=style_desc)
range_y_input = widgets.FloatSlider(value=30, min=1, max=100, step=1, description='Portée Y ($a_y$)', style=style_desc)
range_z_input = widgets.FloatSlider(value=30, min=1, max=100, step=1, description='Portée Z ($a_z$)', style=style_desc)
# Section 2 : Longueurs
lx_input = widgets.FloatSlider(value=10, min=1, max=100, step=1, description='Longueur X ($l_x$)', style=style_desc)
ly_input = widgets.FloatSlider(value=10, min=1, max=100, step=1, description='Longueur Y ($l_y$)', style=style_desc)
lz_input = widgets.FloatSlider(value=10, min=1, max=100, step=1, description='Longueur Z ($l_z$)', style=style_desc)
# Section 3 : Paramètres restants
geometry_input = widgets.Dropdown(options=['ligne', 'surface', 'cube'], value='surface', description='Géométrie', style=style_desc)
n_points_input = widgets.IntSlider(value=5, min=1, max=20, step=1, description='Points Gauss par dim.', style=style_desc)
sill_input = widgets.FloatSlider(value=1.0, min=0.1, max=5.0, step=0.1, description='Plateau (Sill)', style=style_desc)
model_input = widgets.ToggleButtons(options=['sphérique', 'exponentiel', 'gaussien'], value='sphérique', style=style_desc)
def update_n_points_range(change):
geom = change['new']
n_points_input.max = {'ligne': 50, 'surface': 20, 'cube': 10}[geom]
if n_points_input.value > n_points_input.max:
n_points_input.value = n_points_input.max
geometry_input.observe(update_n_points_range, names='value')
update_n_points_range({'new': geometry_input.value})
# --- Bouton d'action ---
calc_button = widgets.Button(description="Exécuter le Calcul", button_style='success', icon='cogs')
output_area = widgets.Output()
def on_calc_button_clicked(b):
with output_area:
clear_output(wait=True)
base_params = {
'geometry': geometry_input.value,
'lx': lx_input.value,
'ly': ly_input.value,
'lz': lz_input.value,
'sill': sill_input.value,
'model_name': model_input.value,
'range_x': range_x_input.value,
'range_y': range_y_input.value,
'range_z': range_z_input.value,
}
# Pour la courbe de convergence
if base_params['geometry'] == 'ligne':
sample_range = np.arange(1, 51)
elif base_params['geometry'] == 'surface':
sample_range = np.arange(1, 21)
else:
sample_range = np.arange(1, 11)
variances = []
for n in sample_range:
var, *_ = block_variance_gauss_quad(**base_params, n_points=n)
variances.append(var)
var_current, X, Y, Z = block_variance_gauss_quad(**base_params, n_points=n_points_input.value)
print(f"Variance estimée du bloc : {var_current:.4f}")
# --- Affichage des figures ---
fig = plt.figure(figsize=(8, 6))
# 1. Convergence
ax1 = fig.add_subplot(1, 2, 1)
ax1.plot(sample_range, variances, '-o', color='royalblue')
ax1.plot(n_points_input.value, var_current, 'o', color='red', markersize=10, label='Config actuelle')
ax1.set_title('Convergence de la variance')
ax1.set_xlabel('Points Gauss par dimension')
ax1.set_ylabel('Variance')
ax1.grid(True)
ax1.legend()
# 2. Points quadrature
ax2 = fig.add_subplot(1, 2, 2, projection='3d' if base_params['geometry'] == 'cube' else None)
if base_params['geometry'] == 'ligne':
ax2.plot(X, np.zeros_like(X), 'o', color='black')
ax2.set_ylim(-0.1 * base_params['lx'], 0.1 * base_params['lx'])
ax2.set_xlabel("X")
elif base_params['geometry'] == 'surface':
ax2.plot(X, Y, 'o', color='black')
ax2.set_aspect('equal', adjustable='box')
ax2.set_xlabel("X")
ax2.set_ylabel("Y")
else:
ax2.scatter(X, Y, Z, c='black', alpha=0.6)
ax2.set_xlabel("X")
ax2.set_ylabel("Y")
ax2.set_zlabel("Z")
ax2.set_box_aspect([1, 1, 1])
ax2.set_title("Points de quadrature")
plt.tight_layout()
plt.show()
calc_button.on_click(on_calc_button_clicked)
# --- UI finale ---
box_anisotropie = widgets.VBox([
widgets.HTML("<b>Anisotropie</b>"),
range_x_input, range_y_input, range_z_input
])
box_longueurs = widgets.VBox([
widgets.HTML("<b>Longueurs du bloc</b>"),
lx_input, ly_input, lz_input
])
box_reste = widgets.VBox([
widgets.HTML("<b>Paramètres supplémentaires</b>"),
geometry_input,
n_points_input,
sill_input,
widgets.HTML("<b>Modèle</b>"),
model_input
])
ui = widgets.VBox([box_anisotropie, box_longueurs, box_reste, calc_button, output_area])
display(ui)
Loading...