Non linear system of equations with discretization on k-space

Question

I want to numerically solve the following system of differential equations at the steady state: \begin{equation} \begin{aligned} \frac{\partial \rho_{11\mathbf{k}}}{\partial t} =& +\frac{i}{\hbar} V_{0}^{2121}\sum_\mathbf{k'}(\rho_{21\mathbf{k'}}\rho_{21\mathbf{k}}^* - \rho_{21\mathbf{k'}}^*\rho_{21\mathbf{k}}) + \\ & -\Gamma (\rho_{1\mathbf{k}}-\rho^0_{11\mathbf{k}})+ \\ & -\eta\Gamma_\mathbf{k}^{out}\rho_{11\mathbf{k}} + (1-\eta)\Gamma_\mathbf{k}^{in}(1-\rho_{11\mathbf{k}})\\ \end{aligned} \end{equation} \begin{equation} \begin{aligned} \frac{\partial \rho_{22\mathbf{k}}}{\partial t} =& -\frac{i}{\hbar} V_{0}^{2121}\sum_\mathbf{k'}(\rho_{21\mathbf{k'}}\rho_{21\mathbf{k}}^* - \rho_{21\mathbf{k'}}^*\rho_{21\mathbf{k}}) + \\ & -\Gamma (\rho_{22\mathbf{k}}-\rho_{22\mathbf{k}}^0)+\\ &+\eta\Gamma_\mathbf{k}^{in}(1-\rho_{22\mathbf{k}}) - (1-\eta)\Gamma_\mathbf{k}^{out}\rho_{22\mathbf{k}}\\ \end{aligned} \end{equation} \begin{equation} \begin{aligned} \frac{\partial \rho_{21\mathbf{k}}}{\partial t} = & - i\omega_{12} \rho_{21\mathbf{k}} + \\ &+ \frac{i}{\hbar} V_{0}^{2121}\sum_\mathbf{k'} \rho_{21\mathbf{k'}}(\rho_{22\mathbf{k}}-\rho_{11\mathbf{k}}) +\\ & -\Gamma_\rho \,\rho_{21\mathbf{k}} \end{aligned} \end{equation}

$V_0^{2121}, \Gamma, \eta, \Gamma_\mathbf{k}^{in}, \Gamma_\mathbf{k}^{out}, \Gamma_\rho$ are parameters.

The variables $\rho_{11\mathbf{k}}, \rho_{22\mathbf{k}}$ and $\rho_{21\mathbf{k}}$ are complex, and defined on a discretized space $\mathbf{k}$.

I implemented a Newton-Raphson method in the complex plane, but I am not sure of the results. Could someone suggest me a package (Python, Matlab, Mathematica?) that I could use to confirm my results? I know that python can solve nonlinear equations, but I am not sure of how to deal with the $\mathbf{k}$ discretization.

Thank you!

Answer 1

The discretization in $\mathbf{k}$-space is just a frequency discretization; if you took the inverse Fourier transform of your equations over $\mathbf{k}$-space, you would be back in "real" space (from the looks of it, your coordinates would be spatial coordinates).

Nonlinear equation solvers operate on systems of equations; how you obtain those equations is mostly immaterial, except for properties of the equations (for instance, conditioning of the Jacobian matrix, avoiding roundoff error, etc.). You could use:

• fsolve in MATLAB
• NSolve in Mathematica
• various functions from scipy.optimize in Python
• petsc4py has more advanced nonlinear solves (Newton-Krylov methods with many different linear and nonlinear preconditioners); it is a Python interface to PETSc

I strongly recommend using nonlinear solvers from a reputable library instead of writing them yourself, to save yourself time.