4
$\begingroup$

First I will generate such kind of monotonic sawtooth like function called func

n = 50;
block = Symmetrize[RandomReal[1., {n, n}], Symmetric[{1, 2}]];
diagF = SparseArray[Band[{1, 1}] -> Normal /@ #] &;
mat = diagF[{block, block, block, block}];
dim = Length@mat;
eigs = Transpose@Sort@Transpose@Eigensystem@mat;
eigsdensity = {eigs[[1]], Abs[eigs[[2]]]^2};
fermiCom = 
  Compile[{{e, _Real}, {u, _Real}, {t, _Real}}, 
   Module[{tmp}, tmp = Exp[(u - e)/t]; tmp = 1./(1. + tmp); 
    tmp = 1. - tmp], RuntimeOptions -> "Speed"];
Clear[fermi];
fermi[e_?NumericQ, u_?NumericQ, t_?NumericQ] := fermiCom[e, u, t];
num = Length@Select[eigs[[1]], # < 0 &] + 1;
Clear[func];
func[u_] := 
 Total[Sum[
   eigsdensity[[2, i]] fermi[eigs[[1, i]], u, 0.0001], {i, 1, dim}]]

The plot of func is like this

Plot[func[u], {u, -3, 3}]

enter image description here

This is what I call monomotic sawtooth like function

Now I want to know at which u value, func[u] is closest to num ( num = Length@Select[eigs[[1]], # < 0 &] + 1)

Normally, FindRoot fastest for this kind of task. However it is not working for this problem

FindRoot[func[u] - num == 0, {u, 0}]

gives error

FindRoot::jsing: Encountered a singular Jacobian at the point {u} = {0.}. Try perturbing the initial point(s). >>

This is because num is a value guaranteed to be between two plateaus (the fermi function is very close to step function, while not exactly the same), so it can't find such a root.

Then, we can transform the problem into a minimization problem. That is to minimize Abs[func[u] - num]

HoweverFindMinimum is also not working.

FindMinimum[Abs[func[u] - num], u]

FindMinimum::fmgz: Encountered a gradient that is effectively zero. The result returned may not be a minimum; it may be a maximum or a saddle point. >>

The only working, I know at the moment is NMinimize

NMinimize[Abs[func[u] - num], u] // AbsoluteTiming
(*{0.934564, {1., {u -> 0.0116856}}}*)

But as you can see NMinimize is really expensive. So I wonder, since the func such a good monotonic property, is there a better way than direct NMinimize?

$\endgroup$
5
$\begingroup$

From the FindRoot documentation:

FindRoot[lhs == rhs, {x, x0, x1}] searches for a solution using $x_0$ and $x_1$ as the first two values of $x$, avoiding the use of derivatives. ... If you specify two starting values, FindRoot uses a variant of the secant method.

Indeed, that does the job pretty much instantly:

FindRoot[func[u] - num == 0, {u, -3, 3}]
(* {u -> 0.00249415} *)
$\endgroup$
  • $\begingroup$ I was googling secant before I post this. Somehow I didn't find how to do it in mathematica.... I should have read the doc more carefully. Thank you very much for quick response : ) $\endgroup$ – matheorem Apr 21 '16 at 7:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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