Search code examples
pythonscipysolverscientific-computing

Scipy.optimize.fsolve() on my complex equation


I try to solve a complex function with scipy. I read that fsolve only works with real numbers. So I made some changes to my formula and now it looks like this:

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)

It's my first time using a solver (and scipy), but I think this is the right way to use it.

Above code gives the following error:

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,).

I dont know how to fix this error and how to get a complex n that solves my equation.


Solution

  • As discussed in the comments, minimize seems to be more appropriate for this kind of problem. The following works and terminates successfully:

    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)
    

    This will print:

     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])
    

    We use method='Nelder-Mead' as this is a gradient-free optimization method which seems to be required for this kind of problem.