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!


1 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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.