为了账号安全,请及时绑定邮箱和手机立即绑定

Scipy.optimize.fsolve() 在我的复杂方程上

Scipy.optimize.fsolve() 在我的复杂方程上

翻翻过去那场雪 2021-07-26 17:47:28
我尝试用 scipy 解决一个复杂的函数。我读到 fsolve 只适用于实数。所以我对我的公式做了一些更改,现在看起来像这样:import numpy as npfrom scipy.optimize import fsolvedelta = 84.37228318652858 * np.pi / 180psi = 55.2217535040673 * np.pi / 180n_S = 2.6726 + 3.0375jphi_i = 70 * np.pi / 180d_L = 300  # thickness of layer in nmn_air = 1  # refractive index of airdef snell(phi, n1, n2):    phi_ref = np.arcsin((n1 / n2) * np.sin(phi))    return phi_refdef fresnel(n1, phi1, n2, phi2):    rs = (n1 * np.cos(phi1) - n2 * np.cos(phi2)) / (        n1 * np.cos(phi1) + n2 * np.cos(phi2))    rp = (n2 * np.cos(phi1) - n1 * np.cos(phi2)) / (        n2 * np.cos(phi1) + n1 * np.cos(phi2))    return rs, rpdef calc_rho(n_k):    n = n_k[0]    k = n_k[1]    n_L = n + 1j * k    phi_L = snell(phi_i, n_air, n_L)    phi_S = snell(phi_L, n_L, n_S)    rs_al, rp_al = fresnel(n_air, phi_i, n_L, phi_L)    rs_ls, rp_ls = fresnel(n_L, phi_L, n_S, phi_S)    beta = 2 * np.pi * d_L * n_L * np.cos(phi_L) / lambda_vac    rp_L = (rp_al + rp_ls * np.exp(-2 * 1j * beta)) / (        1 + rp_al * rp_ls * np.exp(-2 * 1j * beta))    rs_L = (rs_al + rs_ls * np.exp(-2 * 1j * beta)) / (        1 + rs_al * rs_ls * np.exp(-2 * 1j * beta))    rho_L = rp_L / rs_L    return abs(rho_L - rho_given)lambda_vac = 300n_S = 2.6726 + 3.0375jrho_given = np.tan(psi) * np.exp(    1j * delta)  # should be 1.4399295435287844+0.011780279522394433jinitialGuess = [1.5, 0.1]nsolve = fsolve(calc_rho, initialGuess)print(nsolve)这是我第一次使用求解器(和 scipy),但我认为这是使用它的正确方法。上面的代码给出了以下错误:TypeError: fsolve: there is a mismatch between the input and output shape of the 'func' argument 'calc_rho'.Shape should be (2,) but it is (1,).我不知道如何解决这个错误以及如何获得一个复杂的 n 来解决我的方程。
查看完整描述

1 回答

?
波斯汪

TA贡献1811条经验 获得超4个赞

正如评论中所讨论的,minimize似乎更适合此类问题。以下工作并成功终止:


import numpy as np

from scipy.optimize import minimize


delta = 84.37228318652858 * np.pi / 180

psi = 55.2217535040673 * np.pi / 180

n_S = 2.6726 + 3.0375j

phi_i = 70 * np.pi / 180

d_L = 300  # thickness of layer in nm

n_air = 1  # refractive index of air



def snell(phi, n1, n2):

    phi_ref = np.arcsin((n1 / n2) * np.sin(phi))

    return phi_ref



def fresnel(n1, phi1, n2, phi2):

    rs = (n1 * np.cos(phi1) - n2 * np.cos(phi2)) / (

        n1 * np.cos(phi1) + n2 * np.cos(phi2))

    rp = (n2 * np.cos(phi1) - n1 * np.cos(phi2)) / (

        n2 * np.cos(phi1) + n1 * np.cos(phi2))

    return rs, rp



def calc_rho(n_k):

    n = n_k[0]

    k = n_k[1]

    n_L = n + 1j * k

    phi_L = snell(phi_i, n_air, n_L)

    phi_S = snell(phi_L, n_L, n_S)

    rs_al, rp_al = fresnel(n_air, phi_i, n_L, phi_L)

    rs_ls, rp_ls = fresnel(n_L, phi_L, n_S, phi_S)

    beta = 2 * np.pi * d_L * n_L * np.cos(phi_L) / lambda_vac

    rp_L = (rp_al + rp_ls * np.exp(-2 * 1j * beta)) / (

        1 + rp_al * rp_ls * np.exp(-2 * 1j * beta))

    rs_L = (rs_al + rs_ls * np.exp(-2 * 1j * beta)) / (

        1 + rs_al * rs_ls * np.exp(-2 * 1j * beta))

    rho_L = rp_L / rs_L

    return abs(rho_L - rho_given)



lambda_vac = 300

n_S = 2.6726 + 3.0375j

rho_given = np.tan(psi) * np.exp(

    1j * delta)  # should be 1.4399295435287844+0.011780279522394433j

initialGuess = [1.5, 0.1]

nsolve = minimize(calc_rho, initialGuess, method='Nelder-Mead')

print(nsolve)

这将打印:


 final_simplex: (array([[1.50419259, 0.10077416],

       [1.50419591, 0.10076434],

       [1.50421181, 0.1007712 ]]), array([0.00015298, 0.00015507, 0.00022935]))

           fun: 0.00015297827023618208

       message: 'Optimization terminated successfully.'

          nfev: 48

           nit: 25

        status: 0

       success: True

             x: array([1.50419259, 0.10077416])

我们使用的method='Nelder-Mead'是一种无梯度优化方法,这似乎是此类问题所必需的。


查看完整回答
反对 回复 2021-08-03
  • 1 回答
  • 0 关注
  • 291 浏览
慕课专栏
更多

添加回答

举报

0/150
提交
取消
微信客服

购课补贴
联系客服咨询优惠详情

帮助反馈 APP下载

慕课网APP
您的移动学习伙伴

公众号

扫描二维码
关注慕课网微信公众号