Mathematica Asked by V. Arjona on November 2, 2020
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 -> {}]
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]]
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]
Correct answer by Bob Hanlon on November 2, 2020
Get help from others!
Recent Questions
Recent Answers
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP