How to solve self-consistent equations numerically?

Question

I am working in condensed matter, where I'm required to solve an integral self consistently which is of the form, $$ \Delta = \int f(x,y,\Delta,\mu,h)dxdy $$ Basically I need to find value of $\mu,\,\Delta$ such that the above equation is satisfied for the specific value of $h$. Also the solution ($\mu^*, \Delta^*$) should be such that $$ \int g(x,y,\Delta^*,\mu^*,h) dx dy=1 $$ And the $h$ is a free parameter. How to solve such a problem numerically?

What I have tried until now is to take a range of $\mu$ and $\Delta$ and check if both the equations are satisfied. This is usually taking very very long time numerically. If anyhow I reach my solution how to know if that is the true minima of the solution space?

Answer 1

It seems to me that you're looking for constrained optimization. These are often reasonably similar to more common optimization methods, such as Nelder-Mead or Newton's method, except they usually add a penalty function that increases the 'cost' of a point (in $\{\eta,\,\Delta\}$ space) when the constraint is not satisfied (in your case, when the integral over $g$ is not 1).

Python's scipy package has several options in this regard. If you are using another language, the list at the bottom of this Wikipedia entry may be of use.

Some optimization functions in the wild will return a global minima of the parameter, while others will return only a local one; so your key task in using one of the existing packages is to determine if it will return the local or global. That said, I believe that you could use particle swarm methods to increase your chances of finding a true global minima using any of the local minima algorithms.

Answer 2

We can treat this as a problem of solving 2 equations of 2 variables. If we define $F_1 \left(\mu,\Delta\right)=\Delta-\int f\left(x,y,\Delta,\mu,h\right)dxdy$ and $F_2 \left(\mu,\Delta\right)=1-\int g\left(x,y,\Delta,\mu\right)dxdy$ then the problem becomes finding a $\left(\mu^{*},\Delta^{*}\right)$ which is a simultaneous root of both $F_1$ and $F_2$. For some vector function $\vec{F}:\mathbb{R}^{n}\to\mathbb{R}^{n}$, we can find a root (i.e. solve $\vec{F}\left(\vec{x}\right)=\vec{0}$) with the multivariable Newton-Raphson method. Starting with some guess $\vec{x}_{old}$, calculate the Jacobian matrix $J(\vec{x}_{old})$ of $\vec{F}$. If an analytic expression for the Jacobian is known, we use that (evaluate it at $\vec{x}_{old}$). Otherwise, we compute the derivatives in the Jacobian as finite differences, $\frac{\partial F_{i}}{\partial x_{j}}\left(\vec{x}_{old}\right)=\frac{F_{i}\left(\vec{x}_{old}+h\hat{e}_{j}\right)-F_{i}\left(\vec{x}_{old}\right)}{h}$ (here $\hat{e}_j$ is the $j$ unit vector). Now we look for a $\vec{\delta}$ such that $\vec{x}_{new}:=\vec{x}_{old}+\vec{\delta}$ will be closer to a root. So we stipulate $\vec{F}(\vec{x}_{new})=\vec{0}$. Applying a linear approximation around $\vec{x}_{old}$, we have $\vec{0}=\vec{F}(\vec{x}_{new})\approx \vec{F}\left(\vec{x}_{\text{old}}\right)+J\left(\vec{x}_{\text{old}}\right)\vec{\delta}$, and hence to find $\vec{\delta}$ we must solve the matrix equation $J\left(\vec{x}_{\text{old}}\right)\vec{\delta}=-\vec{F}\left(\vec{x}_{\text{old}}\right)$. This is done with Gaussian elimination. Now $\vec{x}_{new}=\vec{x}_{old}+\vec{\delta}$ is our new guess, and we repeat this process until convergence.

In general there can be multiple such roots, and this method will converge to just one of them. This algorithm is detailed well in the book Numerical Recipes (see chapters 9.6, 9.7). The book discusses some silghtly fancier algorithms too, which are likelier to converge when the initial guess is a bit farther from the root. All of the necessary code can also be found there.