0
$\begingroup$

I'm a beginner at Mathematica, and I have a problem finding the roots of a function for a specific variable (using Jens' findAllRoots function).

Due to the nature of my problem, I need to compute these zeros inside a table on the fly. I decided to define a function to calculate them, which depends on a variable. Nevertheless, after some iterations, the function keeps information from the previous value given to the variable. What is going wrong? Thanks for your help.

Here is a sample code:

ClearAll[V, g, tem, lim1, fun]

V = 100;
g = 3;
lim1 = 20;

fun[x_,tem_] = D[ 
                 TT[ x , g ] * fte[ x , +V/2, tem]*( 1 - fte[ x , +V/2, tem]), 
                x ]; (* TT is a transmission function and fte is the Fermi distribution  *)

zeros[tem_] := findAllRoots[
                            fun[x , tem ], { x , -V/2 - 20*tem, V/2 + 20*tem}, 
                           "ShowPlot" -> False]

The derivative (function fun) should be computed only once, keeping x and tem as parameters. The function zeros must be computed every time it is called. My first problem is that zeros[tem] does not work, it does not provide any number. To solve it, I have tried using:

zeros1 = Function[tem,
                findAllRoots[fun[ x , tem ], { x , -V/2 - 20*tem, V/2 + 20*tem},
                             "ShowPlot" -> False]
                 ]

zeros2 =   (findAllRoots[fun[ x , #], {x, -V/2 - 20*#, V/2 + 20*#}, "ShowPlot" -> False]) &

But after some iterations the functions maintain the previous given value

zeros1[45] = {-10.8508, 4.43734*10^-31, 50.562}

zeros1[50] = {-10.8508, 4.43734*10^-31, 50.562}[50]

What is the best form to define this zeros function?

(* The definition of the other functions *)

TT[x_, g_] := e^2/(e^2 + g^2) 

fte[x_, m_, t_] := 1/(Exp[(e - m)/t] + 1)  

SyntaxInformation[findAllRoots] = {"LocalVariables" -> {"Plot", {2, 2}}, 
"ArgumentsPattern" -> {_, _, OptionsPattern[]}};
SetAttributes[findAllRoots, HoldAll];

Options[findAllRoots] = Join[{"ShowPlot" -> False, PlotRange -> All},
FilterRules[Options[Plot], Except[PlotRange]]];

findAllRoots[fn_, {l_, lmin_, lmax_}, opts : OptionsPattern[]] :=
Module[{pl, p, x, localFunction, brackets}, localFunction = ReleaseHold[Hold[fn] /. HoldPattern[l] :> x]; 

If[lmin != lmax, pl = Plot[localFunction, {x, lmin, lmax}, 
Evaluate@FilterRules[Join[{opts}, Options[findAllRoots]], Options[Plot]]];

p = Cases[pl, Line[{x__}] :> x, Infinity];

If[OptionValue["ShowPlot"], Print[Show[pl, PlotLabel -> "Finding roots for this function", ImageSize -> 200, BaseStyle -> {FontSize -> 8}]]], p = {}];

brackets = Map[First, Select[(*This Split trick pretends that two points on the curve are "equal" if the function values have _opposite _ sign.Pairs of such sign-changes form the brackets for the subsequent FindRoot*)

Split[p, Sign[Last[#2]] == -Sign[Last[#1]] &], Length[#1] == 2 &], {2}];

x /. Apply[FindRoot[localFunction == 0, {x, ##1}] &, brackets, {1}] /. x -> {}]
$\endgroup$
3
  • $\begingroup$ Please provide the definitions for all functions (e.g., TT, fte, findAllRoots) $\endgroup$ – Bob Hanlon Sep 22 '20 at 21:06
  • 2
    $\begingroup$ I assume that the definitions for TT and fte should have x instead of e on the right hand side? Also, why can't you just use Solve? $\endgroup$ – Carl Woll Sep 22 '20 at 21:45
  • $\begingroup$ Yes, it was a mistake while I was writing the post. If I am right, I tried using solve as first option, but the answer didn't provide all the possible solutions. $\endgroup$ – V. Arjona Sep 23 '20 at 15:36
2
$\begingroup$
Clear["Global`*"]

TT[x_, g_] := x^2/(x^2 + g^2)

fte[x_, m_, t_] := 1/(Exp[(x - m)/t] + 1)

V = 100;
g = 3;
lim1 = 20;

fun[x_, tem_] = 
 D[TT[x, g]*fte[x, +V/2, tem]*(1 - fte[x, +V/2, tem]), x] // FullSimplify
(*TT is a transmission function and fte is the Fermi distribution*)

(* -((x Sech[(-50 + x)/(
   2 tem)]^2 (-18 tem + x (9 + x^2) Tanh[(-50 + x)/(2 tem)]))/(
 4 tem (9 + x^2)^2)) *)

Using NSolve to find the roots

zeros1[tem_?NumericQ] := 
 NSolve[{fun[x, tem] == 0, -V/2 - 20*tem <= x <= V/2 + 20*tem}, x]

zeros1[45]

(* {{x -> -10.8508}, {x -> 0.}, {x -> 50.562}} *)

x == 0 is always a root

fun[0, tem]

(* 0 *)

zeros1[50]

(* {{x -> -11.5445}, {x -> 0.}, {x -> 50.6887}} *)

To plot zeros1

data = Transpose[
   Table[Thread[{tem, x /. zeros1[tem]}], {tem , 5/2, 50, 5/2}]];

MinMax[#[[All, 2]]] & /@ data

(* {{-11.5445, -2.73266}, {0., 0.}, {50.0018, 50.6887}} *)

ListLinePlot[data,
 Frame -> True,
 FrameLabel -> (Style[#, 14, Bold] & /@ {tem, x}),
 PlotLegends -> Automatic,
 PlotLabel -> Style["fun[x, tem] \[Equal] 0", 14, Bold]]

enter image description here

Comparing with ContourPlot

ContourPlot[fun[x, tem] == 0, {tem, 5/2, 50}, {x, -15, 54},
 PlotRange -> {{-1, 51}, {-15, 54}},
 WorkingPrecision -> 50,
 PlotPoints -> 50,
 MaxRecursion -> 5,
 FrameLabel -> (Style[#, 14, Bold] & /@ {tem, x}),
 PlotLabel -> Style["fun[x, tem] \[Equal] 0", 14, Bold],
 AspectRatio -> 1/GoldenRatio]

enter image description here

$\endgroup$
1
  • $\begingroup$ Really appreciate it! This works nicely. $\endgroup$ – V. Arjona Sep 25 '20 at 18:43

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.