Scipy.optimize.fsolve () в моем сложном уравнении - PullRequest
Новая
       45

Scipy.optimize.fsolve () в моем сложном уравнении

0 голосов
/ 30 октября 2018

Я пытаюсь решить сложную функцию с помощью scipy. Я читал, что fsolve работает только с действительными числами. Поэтому я внес некоторые изменения в свою формулу, и теперь она выглядит следующим образом: 

import numpy as np
from scipy.optimize import fsolve

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 = 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 Ответ

0 голосов
/ 31 октября 2018

Как обсуждалось в комментариях, 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', так как это метод оптимизации без градиента, который, по-видимому, необходим для решения этой проблемы.

...