I have the following equation (the Kurz-Giovanola-Trivedi model [1]) $$ v^2 \frac{\pi^2 \Gamma}{P^2 D^2} + v \frac{mC_0(1-k)\xi}{D[1-(1-k)Iv(P)]} + G = 0, $$ where $Iv(P)=P \cdot \exp(P) \cdot E(P)$, $P=\frac{Rv}{2D}$, $E(P)=\int_{P}^{\infty} \frac{\exp(-x)}{x} dx$ is the exponential integral function, $\xi=1-\frac{2k}{[1+(\frac{2\pi}{P})^2]^{\frac{1}{2}}-1+2k}$, $v$ and $R$ are unknowns, and everything else is constant. To solve this equation, I use MATLAB (the core function is fzero) and suppose that $v$ is an array within a range $(10^{-2};10^2)$. The fzero function attempts to find a zero of one equation with one variable and might be called with either a one-element starting point or a two-element vector (starting interval).

My question is: how to effectively find a starting point?

I have started by solving this equation relative to $P$, and it was relatively easy to determine starting points $P_0$ for different intervals of the $v$ array by trial-and-error method. After finding $P$, it was easy to find $R$.

However, now I need to solve the equation relative to $R$ directly, without the $P$-step. The trial-and-error method is not the best option here, probably due to the type of the $R(v)$-dependence (please see below).

enter image description here


1 Answer 1


I recommend using a global optimization solver.

You should be able to solve this, or at least try, using the BMIBNB global optimization solver https://yalmip.github.io/solver/bmibnb/ under YALMIP https://yalmip.github.io/tutorial/basics/ under MATLAB.

In particular, note that YALMIP supports the exponential integral, expint, https://yalmip.github.io/command/expint/

For this problem:

Declare the variables v and R as sdpvar.

Enter the equation, and lower and upper bounds on the variables, all inside [] , as the constraints in the 1st argument of the optimize command

Set the objective to [] as the 2nd argument in the optimize command.

Set 'bmibnb' as the solver using sdpsettings as the 3rd argument in the optimize command.

You do not need a starting value. You will need a local nonlinear solver installed under MATLAB which will be called by BMIBNB, as well as a MILP (Mixed Integer Linear Programming) solver to be called by BMIBNB.

  • 1
    $\begingroup$ I upvoted with the note that, this might be very expensive for their problem. I would suggest, in addition to Mark's suggestion, to somehow linearize system, do one step and use it as an initial condition. For example, consider the following; $$v_0^2 \frac{\pi^2 \Gamma}{P^2 D^2} + v_0 \frac{mC_0(1-k)\xi}{D[1-(1-k)Iv(P)]} + G = 0$$ to obtain a guess for $R$ and $$v^2 \frac{\pi^2 \Gamma}{P_0^2 D^2} + v \frac{mC_0(1-k)\xi}{D[1-(1-k)Iv(P_0)]} + G = 0$$ to obtain a guess for $v$. $v_0$ and $P_0$ (hence $R_0$) can be chosen arbitrarily, $1$ probably would work. Or you can do multiple iterations... $\endgroup$ Aug 5, 2020 at 21:56

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.