假設您有一個離散的r數組，并且該數組的每個值z(r) 。 您想擬合曲線以估計非球面透鏡的參數。 我將使用lmfit提到這里展現做到這一點使用python的一種方式。

導入用于此的模塊：

import numpy as np

import matplotlib.pyplot as plt

from lmfit import Model, Parameters

定義非球面鏡的功能：

def asphere_complete(x, r0, k, a2, a4, a6, a8, a10, a12):

r_squared = x ** 2.

z_even_r = r_squared * (a2 + (r_squared * (a4 + r_squared * (a6 + r_squared * (a8 + r_squared * (a10 + (r_squared * a12)))))))

square_root_term = 1 - (1 + k) * ((x / r0) ** 2)

zg = (x ** 2) / (r0 * (1 + np.sqrt(square_root_term)))

return z_even_r + zg

由于您未提供任何數據，因此我將使用以下內容創建一些示例數據，包括人為噪聲：

def generate_dummy_data(x, asphere_parameters, noise_sigma, seed=12345):

np.random.seed(seed)

return asphere_complete(x, **asphere_parameters) + noise_sigma * np.random.randn(x.shape[0])

以下函數進行擬合并繪制結果曲線：

def fit_asphere(r, z, fit_parameters):

# create two subplots to plot the original data and the fit in one plot and the residual in another

fig, axarr = plt.subplots(1, 2, figsize=(10, 5))

fit_plot = axarr[0]

residuum_plot = axarr[1]

# configure first plot:

fit_plot.set_xlabel("r")

fit_plot.set_ylabel("z")

fit_plot.grid()

# configure second plot:

residuum_plot.set_xlabel("r")

residuum_plot.set_ylabel("$\Delta$z")

residuum_plot.grid()

# plot original data

fit_plot.plot(r, z, label="Input")

# create an lmfit model and the parameters

function_model = Model(asphere_complete)

# The fitting procedure may throw ValueErrors, if the radicand gets negative

try:

result = function_model.fit(z, fit_parameters, x=r)

# To plot the resulting curve remove the parameters which were just used for the constraints

opt_parameters = dict(result.values)

opt_parameters.pop('r_max', None)

opt_parameters.pop('radicand', None)

# calculate z-values of fitted curve:

z_fitted = asphere_complete(r, **opt_parameters)

# calculate residual values

z_residual = z - z_fitted

# plot fit and residual:

fit_plot.plot(r, z_fitted, label="Fit")

residuum_plot.plot(r, z_residual, label="Residual")

# legends:

fit_plot.legend(loc="best")

residuum_plot.legend(loc="best")

print(result.fit_report())

except ValueError as val_error:

print("Fit Failed: ")

print(val_error)

要設置示例數據的參數，我使用lmfit的Parameters對象：

if __name__ == "__main__":

parameters_dummy = Parameters()

parameters_dummy.add('r0', value=-34.4)

parameters_dummy.add('k', value=-0.98)

parameters_dummy.add('a2', value=0)

parameters_dummy.add('a4', value=-9.67e-9)

parameters_dummy.add('a6', value=1.59e-10)

parameters_dummy.add('a8', value=-5.0e-12)

parameters_dummy.add('a10', value=0)

parameters_dummy.add('a12', value=-1.0e-19)

創建示例數據：

r = np.linspace(0, 35, 1000)

z = generate_dummy_data(r, parameters_dummy, 0.00001)

使用lmfit而不是scipy的curve_fit是，平方根的radicand可能變為負數。 我們需要確保：

為此，我們需要提到定義的約束在這里 。 讓我們開始定義我們要在擬合中使用的參數。 基本半徑的添加非常簡單：

parameters = Parameters()

parameters.add('r0', value=-30, vary=True)

為了服從不等式，請添加變量radicand ，該變量不得小于零。 而不是讓k參與擬合常態，而是使r0取決于r0 ， r_max和radicand 。 我們需要使用r_max因為對于最大r ，不等式最成問題。 解決k的不等式導致

用作下面的expr 。 我使用布爾標志來打開/關閉約束：

keep_radicand_safe = True

if keep_radicand_safe:

r_max = np.max(r)

parameters.add('r_max', r_max, vary=False)

parameters.add('radicand', value=0.98, vary=True, min=0)

parameters.add('k', expr='(r0/r_max)**2*(1-radicand)-1')

else:

parameters.add('k', value=-0.98, vary=True)

其余參數可以直接添加：

parameters.add('a2', value=0, vary=False)

parameters.add('a4', value=0, vary=True)

parameters.add('a6', value=0, vary=True)

parameters.add('a8', value=0, vary=True)

parameters.add('a10', value=0, vary=False)

parameters.add('a12', value=0, vary=True)

現在我們準備開始并獲得我們的結果：

fit_asphere(r, z, parameters)

plt.show()

在控制臺上，您應該看到輸出：

[[Variables]]

r0: -34.3999435 +/- 6.1027e-05 (0.00%) (init = -30)

r_max: 35 (fixed)

radicand: 0.71508611 +/- 0.09385813 (13.13%) (init = 0.98)

k: -0.72477176 +/- 0.09066656 (12.51%) == '(r0/r_max)**2*(1-radicand)-1'

a2: 0 (fixed)

a4: 7.7436e-07 +/- 2.7872e-07 (35.99%) (init = 0)

a6: 2.5547e-10 +/- 6.3330e-11 (24.79%) (init = 0)

a8: -4.9832e-12 +/- 1.7115e-14 (0.34%) (init = 0)

a10: 0 (fixed)

a12: -9.8670e-20 +/- 2.0716e-21 (2.10%) (init = 0)

使用我上面使用的數據，如果將keep_radicand_safe設置為False ，則應該看到擬合失敗。

總結

以上是生活随笔為你收集整理的matlab做横截面回归,matlab - 将横截面表面轮廓拟合到通用的已知公式以获得系数并对表面进行数学建模 - 堆栈内存溢出...的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。