Почему аппроксимация кривой иногда работает, а не работает в интерактивном моделировании Python? ⇐ Python

[Anonymous]

Я относительно новичок в графическом интерфейсе Python и работаю над простым проектом, связанным с моими занятиями по физике.
Я пытаюсь смоделировать экспериментальные данные резерфордовского рассеяния с помощью Монте Карло, а также использование ipywidgets для интерактивного моделирования, где вы можете использовать ползунки для изменения значений в уравнении.
Проблема заключается в том, что когда я пытаюсь выполнить аппроксимация кривой смоделированных данных, изменение параметры значений N_i/Кинетическая энергия/Толщина/Элемент фольги, подгонка кривой иногда работает, но когда вы меняете параметры на что-то другое, она выравнивается после 10 градусов (x_values). Затем, если вы измените параметры обратно на прежние значения, иногда они будут давать сбой или работать снова.
Я искал в Интернете, и там говорится, что проблема может возникнуть из-за отключения исходных параметров, но когда Я пытаюсь использовать другие значения или каким-либо образом обновлять код, но это все равно приводит к нарушению подгонки кривой. Это код, наиболее близкий к коду эксперимента по резерфордскому рассеянию.
Это код, с которым я работаю. При каждом нажатии кнопки «Выполнить» отображаются смоделированные данные уравнения Резерфорда, которые были согласованы для различных значений N_i. Затем вы нажимаете кнопку «Аналитическая кривая», которая дает теоретическую кривую уравнения. Затем кнопка «Подогнать кривую» должна дать кривую, которая точно соответствует смоделированным данным. Есть ползунки для значений N_i, кинетической энергии, толщины фольги и элемента, из которого она сделана.

Код: Выделить всё

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
import ipywidgets as widgets
from IPython.display import display, Latex
import time

# Constants
Z1 = 2  # Atomic number of alpha particle
e = 1.602e-19  # Elementary charge (C)
epsilon_0 = 8.854e-12  # Vacuum permittivity (F/m)
M_alpha = 6.644e-27  # Mass of alpha particle (kg)

# Simulation Function: Rutherford's scattering formula
def rutherford_scattering(theta, N_i, a, t, r, Z2, EK):
"""Compute N(theta) based on Rutherford's law."""
first = (N_i * a * t) / (16 * r**2)
second = ((Z1 * Z2 * e**2) / (4 * np.pi * epsilon_0 * EK))**2
return first * second * (1 / np.sin(theta / 2)**4)

# Run the simulation (Monte Carlo simulation)
def monte_carlo_simulation(analytical, theta, EK_slider_value, Ni_slider_value, noise_level=0.2):
"""Add random noise to the analytical data to simulate experimental results."""
# Apply more noise at larger angles (factor based on angle)
angle_factor = 1 / np.sin(theta / 2)**4  # Using angle factor for noise

# Scale noise with respect to kinetic energy (lower noise at higher KE)
# Adjust the scaling factor based on a more realistic energy value in Joules (not multiplied by 1e-13 here)
scaled_noise_level = 1 / np.sqrt(EK_slider_value)  # Inversely proportional to KE

# Noise scaling based on the number of particles (Ni)
noise_scaling = 1 / np.sqrt(Ni_slider_value)  # More noise for smaller Ni

# Combine noise scaling
total_noise_level = scaled_noise_level + noise_scaling + angle_factor

# Generate noise (normally distributed)
noise = np.random.normal(0, total_noise_level, size=theta.shape)  # Generate noise for each data point

# Add noise to the analytical values
simulated_data = analytical + noise

# Ensure that N(theta) does not go below zero
simulated_data = np.maximum(simulated_data, 0)  # Clamp values to be >= 0

return simulated_data

# Widgets for input parameters
Ni_slider = widgets.IntSlider(value=13000, min=6000, max=20000, step=500, description="Ni: No.  alpha particles", style={'description_width': 'initial'}, layout={'width': '400px'})
t_input = widgets.FloatText(value=400, min=200, max=1400, step=100, description="t: thickness of foil (nm)", style={'description_width': 'initial'})
EK_slider = widgets.FloatSlider(value=5, min=1, max=10, step=1, description="EK: Kinetic energy (J)", style={'description_width': 'initial'}, layout={'width': '400px'})
Z2_dropdown = widgets.Dropdown(
options={
'Gold (79)': 79,
'Silver (47)': 47,
'Copper (29)': 29,
'Aluminum (13)': 13,
'Iron (26)': 26,
'Titanium (22)': 22,
'Carbon (6)': 6,
'Oxygen (8)': 8,
'Lead (82)': 82,
'Mercury (80)': 80
},
value=79,
description="Z2: Element of foil",
style={'description_width': 'initial'}
)

# Buttons for running simulation, curve fitting, and displaying results
simulate_button = widgets.Button(description="Run Monte Carlo Simulation")
add_analytical_button = widgets.Button(description="Add Analytical Curve")
fit_button = widgets.Button(description="Fit Curve")

# Output plots and error messages
output_plot = widgets.Output()
output_error = widgets.Output()
output_analytical = widgets.Output()

# Layout for widgets
ui = widgets.VBox([
widgets.HBox([Ni_slider]),
widgets.HBox([EK_slider]),
widgets.VBox([t_input]),
widgets.VBox([Z2_dropdown]),
widgets.HBox([simulate_button, add_analytical_button, fit_button]),
output_plot,  # Main plot for simulation
output_analytical,  # Dedicated plot for analytical curve
output_error
])

# Global variables for storing simulated data and analytical data
theta_values = np.linspace(np.radians(1), np.radians(179), 500)
simulated_data = None
analytical_data = None

# Debounce function to prevent rapid successive calls
def debounce(f, wait=1):
"""Debounce a function to delay execution until user stops adjusting sliders."""
last_call_time = [0]  # Use mutable object to retain state between calls

def wrapped_function(change):
now = time.time()
if now - last_call_time[0] > wait:
last_call_time[0] = now
f(change)

return wrapped_function

@debounce # Updated functions with debouncing applied
# Run the simulation (Monte Carlo simulation)
def run_simulation(change):
global simulated_data
with output_plot:
output_plot.clear_output()  # Clear previous content in the output area

# Read values from widgets
Ni = Ni_slider.value
t_nm = t_input.value  # Thickness in nanometers
t = t_nm * 1e-9  # Convert thickness to meters
EK = EK_slider.value * 1e-13  # Ensure the kinetic energy is properly converted to joules
Z2 = Z2_dropdown.value
r = 0.1  # Fixed distance from foil (m)
a = 1e28  # Approx.  atoms/m^3 for dense materials

# Compute the theoretical Rutherford scattering distribution (analytical curve)
analytical = rutherford_scattering(theta_values, Ni, a, t, r, Z2, EK)
global analytical_data
analytical_data = analytical  # Store analytical data for future use

# Simulate the data by adding noise to the theoretical model
simulated_data = monte_carlo_simulation(analytical_data, theta_values, EK_slider.value, Ni_slider.value, noise_level=0.1)

# Scatter plot: Simulated data
plt.figure(figsize=(8, 6))
plt.plot(np.degrees(theta_values), simulated_data, '.', color='green', label="Simulated Data")
plt.xlabel("Scattering Angle (degrees)")
plt.ylabel(r"$N(\theta)$")
plt.yscale('log')
plt.title("Simulated Data")
plt.grid()
plt.tight_layout()
plt.legend()
plt.show()

@debounce
# Add analytical Rutherford scattering curve to the simulated data
def add_analytical_curve(change):
global simulated_data, analytical_data
if analytical_data is None or simulated_data is None:
with output_error:
output_error.clear_output()
print("Run the simulation first to generate simulated data and analytical curve.")
return  # Ensure both simulated and analytical data exist

# Use the main output area to display both curves together
with output_plot:
output_plot.clear_output()  # Clear previous content in the main plot area

# Create the plot
plt.figure(figsize=(8, 6))
plt.plot(np.degrees(theta_values), simulated_data, '.', color='green', label="Simulated Data")
plt.plot(np.degrees(theta_values), analytical_data, '-', color='blue', label="Analytical Curve")
plt.xlabel("Scattering Angle (degrees)")
plt.ylabel(r"$N(\theta)$")
plt.yscale('log')
plt.title("Simulated Data with Analytical Curve")
plt.grid()
plt.tight_layout()
plt.legend()
plt.show()

# Curve fitting function
def fit_curve_function(theta, A, B):
"""Model function for fitting a curve to the data."""
return A / np.sin(theta / 2)**4 + B

@debounce
# Fit the simulated data
def fit_simulation(change):
global simulated_data, analytical_data
if analytical_data is None or simulated_data is None:
with output_error:
output_error.clear_output()
print("Run the simulation first to generate simulated data and analytical curve.")
return

try:
# Better initial parameter estimates
A_init = np.max(simulated_data) * (np.sin(np.pi/4))**4
B_init = np.min(simulated_data)

# Try multiple initial guesses if the first fit fails
initial_guesses = [
[A_init, B_init],
[A_init * 10, B_init],
[A_init * 0.1, B_init],
[A_init, 0]
]

best_fit = None
best_residual = np.inf

for guess in initial_guesses:
try:
# Wider bounds
bounds = ([0, 0], [np.inf, np.max(simulated_data)])

# Add method='trf' and increase max iterations
popt_try, pcov_try = curve_fit(
fit_curve_function,
theta_values,
simulated_data,
p0=guess,
bounds=bounds,
method='trf',
maxfev=20000
)

# Calculate fit quality
fitted = fit_curve_function(theta_values, *popt_try)
residual = np.sum((simulated_data - fitted)**2)

if residual <  best_residual:
best_residual = residual
best_fit = (popt_try, pcov_try)

except:
continue

if best_fit is None:
raise RuntimeError("Failed to find a good fit with any initial parameters")

popt, pcov = best_fit

# Rest of your plotting code remains the same
A, B = popt
fitted_curve = fit_curve_function(theta_values, *popt)

with output_plot:
output_plot.clear_output()

plt.figure(figsize=(8, 6))
plt.plot(np.degrees(theta_values), simulated_data, '.', color='green', label="Simulated Data")
plt.plot(np.degrees(theta_values), analytical_data, '-', color='blue', label="Analytical Curve")
plt.plot(np.degrees(theta_values), fitted_curve, '--', color='red', label="Fitted Curve")
plt.xlabel("Scattering Angle (degrees)")
plt.ylabel(r"$N(\theta)$")
plt.yscale('log')
plt.title("Curve Fitting to Simulated Data")
plt.legend()
plt.grid()

plt.tight_layout()
plt.show()

# Display fit parameters and covariance matrix
with output_error:
output_error.clear_output()
print(f"Fitted Parameters: A = {A:.3e}, B = {B:.3e}")
print(f"Covariance Matrix: \n{pcov}")

except Exception as e:
with output_error:
print(f"Error during curve fitting: {e}")

# Connect buttons
simulate_button.on_click(run_simulation)
add_analytical_button.on_click(add_analytical_curve)
fit_button.on_click(fit_simulation)

# Display the widgets
display(ui)

Этот код приводит к тому, что первые две кнопки работают совершенно нормально, а ползунки показывают графики, относительно близкие к тому, как это было бы в реальной жизни. Однако иногда при нажатии кнопки «Подогнать кривую» кривая соответствует ровным линиям, превышающим 10 градусов. Это будет работать для некоторых значений, но затем не будет работать при повторном использовании этих значений.
Что-то не так с форматом или макетом, что вызывает эту проблему? Или это логарифмическая проблема, поскольку мои значения y являются логарифмическими для больших значений? Я попробовал все, что мог, и не знаю, куда идти дальше.
Спасибо.
редактирование: я также заметил проблема с уровнем шума для KE. Обычно при высоком KE должно быть меньше шума, поскольку меньше разброс, а при низком KE необходимо больше шума для большего разброса. Я заметил, что он делает прямо противоположное, но я сделал обратное, и это все равно не меняет ситуацию, поэтому я предполагаю, что он накапливает шум откуда-то еще

Подробнее здесь: https://stackoverflow.com/questions/792 ... n-simulati