You need to enable JavaScript to run this app.
优惠活动
大模型
产品
解决方案
定价
更多

使用lmfit库拟合双缝实验数据时遇报错求助

双缝实验数据拟合问题(基于lmfit库)

我有双缝实验的实验数据,希望用lmfit库拟合数值曲线,因为相比Python其他拟合工具,我更偏好它的拟合结果输出形式。以下是我的代码和遇到的问题:

初始代码与Key Error: N报错

from PIL import Image
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy.signal import argrelextrema
from lmfit.models import Model

#Determine the observed intensity distribution

img_data = pd.read_csv('Image31 intensity profile.csv')

pixels = np.array(img_data['Distance_(_)'])
counts = np.array(img_data['Gray_Value'])

norm_intensity = counts / np.linalg.norm(counts)


# Double slit interferometry coherent beam

def keV2m(keV):
    """Converts Photon Energy [keV] to wavelength [m].
    .. note:: Calculation in Vacuum.
    """
    wl = 1./(keV*1000)*4.1356*(10**(-7))*2.9998
    return wl

ener = 10#kev
I_0 = 1
a = 2.0e-6 # slit width in  m
d = 10e-6 # slit seperation in m
L = 6.0 # slit to detec dist in m
detec_pixel_size = 3.1e-6 #pixel size

x = np.arange(-500,500)*detec_pixel_size + 1e-9

def coh_inter_patt(x):
    lam = keV2m(ener)

    sinq_part = ((np.sin((np.pi*x*a)/(lam*L))/((np.pi*x*a)/(lam*L)))**2)

    sinq_part[sinq_part == np.nan] = 1

    inter_pat = I_0*(np.cos((np.pi*x*d)/(lam*L))**2)*sinq_part

    return inter_pat

def coh_inter_patt_sft(x,detec_pixel_size,sft=10):
    lam = keV2m(ener)

    sinq_part = ((np.sin((np.pi*(x+(sft*detec_pixel_size))*a)/(lam*L))/((np.pi*(x+(sft*detec_pixel_size))*a)/(lam*L)))**2)

    sinq_part[sinq_part == np.nan] = 1

    inter_pat = I_0*(np.cos((np.pi*(x+(sft*detec_pixel_size))*d)/(lam*L))**2)*sinq_part

    return inter_pat


inter_pat = coh_inter_patt(x)
inter_pat_sft = coh_inter_patt_sft(x,detec_pixel_size,sft=10)

max_vals = np.sort(argrelextrema(inter_pat, np.greater)[0])

min_vals = np.sort(argrelextrema(inter_pat, np.less)[0])

visibility = (inter_pat[max_vals] - inter_pat[min_vals[1:]])/(inter_pat[max_vals] + inter_pat[min_vals[1:]])


# 'Blurred' patterns
extension = 40 # 10*detec_pixel_size long extended source
blur_patt  = np.zeros_like(coh_inter_patt(x))

for i in range(extension):

    blur_patt += coh_inter_patt_sft(x,detec_pixel_size,sft=i)
    blur_patt += coh_inter_patt_sft(x,detec_pixel_size,sft=-i)

blur_patt = blur_patt/np.linalg.norm(blur_patt)

max_vals = np.sort(argrelextrema(inter_pat, np.greater)[0])
min_vals = np.sort(argrelextrema(inter_pat, np.less)[0])

visibility = (inter_pat[max_vals] - inter_pat[min_vals[1:]])/(inter_pat[max_vals] + inter_pat[min_vals[1:]])

blur_max_vals = np.sort(argrelextrema(blur_patt, np.greater)[0])
blur_min_vals = np.sort(argrelextrema(blur_patt, np.less)[0])
blur_visibility = (blur_patt[max_vals] - blur_patt[min_vals[1:]])/(blur_patt[max_vals] + blur_patt[min_vals[1:]])

#Look at plot output
fig = plt.figure(figsize=(10,6))
ax1 = fig.add_subplot(111)
ax1.plot(pixels, norm_intensity, color='red', label='Observed')
ax1.plot(blur_patt,label='Analytical, vis %f'%blur_visibility.mean(), color='k')
ax1.legend()


#iterative curve fitting

def Blurred_fringes(N):
  extension = N # 10*detec_pixel_size long extended source
  blur_patt  = np.zeros_like(coh_inter_patt(x))

  for i in range(extension):

    blur_patt += coh_inter_patt_sft(x,detec_pixel_size,sft=i)
    blur_patt += coh_inter_patt_sft(x,detec_pixel_size,sft=-i)

  #inter_pat = inter_pat/np.linalg.norm(inter_pat)
  blur_patt = blur_patt/np.linalg.norm(blur_patt)
  return blur_patt

model = Model(Blurred_fringes)

params1 = model.make_params(N=35)
results1 = model.fit(counts, params1, x=pixels)

运行上述代码时出现Key Error: N报错。我清楚lmfit的fit()函数需要make_params()生成的字典对象,且模型函数的第一个参数会被默认视为自变量。

修改后的代码与NaN值报错

为解决上述问题,我修改了曲线拟合部分的代码:

X = np.arange(0,1000,1)

def Blurred_fringes(x, N):
  extension = N # 10*detec_pixel_size long extended source
  blur_patt  = np.zeros_like(coh_inter_patt(x))

  for i in range(int(extension)):

    blur_patt += coh_inter_patt_sft(x,detec_pixel_size,sft=i)
    blur_patt += coh_inter_patt_sft(x,detec_pixel_size,sft=-i)

  #inter_pat = inter_pat/np.linalg.norm(inter_pat)
  blur_patt = blur_patt/np.linalg.norm(blur_patt)
  return blur_patt


test_y = Blurred_fringes(x, 20)

plt.plot(X,test_y)
plt.show()


model = Model(Blurred_fringes)

params1 = model.make_params(N=35)
results1 = model.fit(counts, params1, x=pixels)

修改后运行出现**ValueError: The model function generated NaN values and the fit aborted! Please check your model function and/or set boundaries on parameters where applicable.**报错,同时伴随以下运行警告:

:30: RuntimeWarning: invalid value encountered in divide
sinq_part = ((np.sin((np.pi*x*a)/(lam*L))/((np.pi*x*a)/(lam*L)))**2)
:41: RuntimeWarning: invalid value encountered in divide
sinq_part = ((np.sin((np.pi*(x+(sft*detec_pixel_size))*a)/(lam*L))/((np.pi*(x+(sft*detec_pixel_size))*a)/(lam*L)))**2)

我知道sinc函数容易出现除零错误,尝试用sinq_part[sinq_part == np.nan] = 1处理,但问题仍未解决。

最终目标

先让函数完成拟合,再添加迭代步骤,通过查看一系列N值的置信水平找到最优拟合结果。N是唯一的拟合参数,其他提前定义的变量均为固定实验参数(可能因误差略有变动)。恳请提供解决思路。

内容的提问来源于stack exchange,提问作者Baba Booey

相关产品推荐
方舟 Agent Plan

超全模态模型 × Harness 升级,最新支持 Deepseek-V4.1-Flash、GLM-5.3 系列、Doubao-Seedream-5.0-pro、Kimi-K3 (部分), 限时 9.9 元起

最近更新时间:2026.06.29 05:34:57