# How to get accurate solutions from FindRoot?

I am a new user in Wolfram language.

Please, I am trying to calculate accurate solutions $R$ of an algébric equation which contains the Dawson function. The equation is:

$$1+\frac{\theta }{k^2}+\frac{R^2}{k^2}\left[\frac{\sqrt{\frac{2}{\theta }} \left(\frac{1}{2} \theta ^2 \left(\frac{R^2}{2}+1\right)^2+\theta \left(\frac{R^2}{2}+1\right)+1\right) \left(-2 F\left(\sqrt{\frac{\theta }{2}} R\right)\right)}{R}+1+\theta \left(\frac{1}{4 \theta }+\frac{R^2}{4}+1\right)\right]=0$$

where,

$F$ : is the Dawson function

$1<k<3$ : real

$R<1$ : real

$\theta \ggg 1$ : real >>1

To do it, I think we have to use Findroot. The code that I used is the following:

eq[R_?NumericQ, k_?NumericQ, \[Theta]_?NumericQ]:=1+\[Theta]/k^2+R^2/k^2 (1+\[Theta](1+1/(4\[Theta])+R^2/4)+(Sqrt[2/\[Theta]] 1/R (1+\[Theta](1+(R^2)/2)+\[Theta]^2/2 (1+(R^2)/2)^2)(-2 DawsonF[Sqrt[\[Theta]/2]R])))

solr[k_,\[Theta]_]:=Re[R/.FindRoot[eq[R,k,\[Theta]],{R,1.1/k},AccuracyGoal->16,PrecisionGoal->16,WorkingPrecision->25]]

The solutions $R$ that I'm looking for must satisfy: for all values $1<k<3$ we have always $R<1$.

But for $k<2.1$ all solutions satisfy $R>>1$, so the condition $R<1$ is not satisfied.

Block[{\[Theta]=100},Table[{k,solr[k,\[Theta]]},{k,1,3,0.1}]]
(* {{1.,-402.0221127442610480453402},{1.1,162.7992531392247315125129},{1.2,560.0440880401854997217024},{1.3,149.3137124001590352102585}, {1.4,927.9240354439233066791802},{1.5,1204.227451826562972959649},{1.6,  352.3919638872815653450775},{1.7,-578.9018078239970096174140},{1.8,2.284779688985257597042901},{1.9,1.366209370422262795909812},{2.,1.059718464650025788056352},{2.1,0.8928149869232555556054552},{2.2,0.7840644848720738071277466},{2.3,0.7061209842390898952755792},{2.4,0.6468555823332280593457180},{2.5,0.5999539682600659528508929},{2.6,0.5617652337478441609444316},{2.7,0.5300143374213611634335440},{2.8,0.5031957080178527133949467},{2.9,0.4802563992569976549295289},{3.,0.4604211571053869996277583}} *)

I think this is due to the precision that I can't control.

Any help please on this problem!

• Try plotting the left side of your equation for some sample parameters to see where its zeros actually are. Also, although it may not matter, do not include \[ThinSpace] in your code. It is a formatting quantity. May 15, 2018 at 19:39
• Welcome to Mathematica.SE! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Take the tour and check the faqs! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! May 16, 2018 at 1:04
• Of course. Thank you. May 16, 2018 at 5:21
• FindRoot[eq[R, k, θ], {R, 1/2, 0, 1}] will not give roots outside 0 < R < 1, but it will quit if its search reaches either end-point, returning the value of that end-point. So, I do not view this as much of an improvement. Instead, use (for instance) R /. FindInstance[(eq1 /. {k -> 3, θ-> 100}) == 0 && 0 < R < 1, {R}, 5] // N, as suggested by @Mariusz. There is no answer to your comment of which of the two solutions is better. Both are valid and can be found by FindRoot too, although with more effort. May 21, 2018 at 23:14
• Reduce[(eq1 /. {k -> 3, θ -> 100}) == 0 && 0 < R < 1, {R}] // N also works and may use the same algorithm in this case. May 21, 2018 at 23:28

eq = 1 + θ/k^2 + (R (5 R + 4 R θ + R^3 θ - Sqrt[2] Sqrt[1/θ] (8 +
4 (2 + R^2) θ + (2 + R^2)^2 θ^2) DawsonF[(R Sqrt[θ])/Sqrt[2]]))/(4 k^2);

n = 5;(*Increase the value to find more solution*)

FindInstance[eq == 0 && 1 < k < 3 && 0 < R < 1 && 1 < θ < 1000,
{R, k, θ}, n, RandomSeeding -> Automatic] // N // MatrixForm

For: n=15

n = 15;
FindInstance[
eq == 0 && 1 < k < 3 && 0 < R < 1 && 1 < θ < 1000, {R,
k, θ}, n, RandomSeeding -> Automatic,
WorkingPrecision -> 20] // MatrixForm

ContourPlot[Evaluate@Table[(eq /. θ -> j) == 0, {j, 1, 1000, 50}], {k, 1,
3}, {R, 0, 1}, FrameLabel -> Automatic]

ContourPlot[Evaluate@Table[(eq /. θ -> j) == 0, {j, 1, 100}], {k, 1,
3}, {R, 0, 1}, FrameLabel -> Automatic]

OP request: Solution by FindRoot (borrowed code form user: Henrik):

Block[{R, k, \[Theta]},
eq = {R, k, \[Theta]} \[Function]
Evaluate@
Simplify[
1 + \[Theta]/k^2 +
R^2/k^2 (1 + \[Theta] (1 + 1/(4 \[Theta]) + R^2/4) + (Sqrt[
2/\[Theta]] 1/
R (1 + \[Theta] (1 + (R^2)/2) + \[Theta]^2/
2 (1 + (R^2)/2)^2) (-2 DawsonF[
Sqrt[\[Theta]/2] R]))), \[Theta] > 0];]
solr = {k, \[Theta]} \[Function]
Block[{R}, R /. FindRoot[eq[R, k, \[Theta]], {R, 1/10, 3/4}, Method -> "Secant"]] // Quiet;
data = Block[{\[Theta] = 10}, Table[{k, solr[k, \[Theta]]}, {k, 1, 3, 1/10}]] // MatrixForm

As you can see values for R is not in domain: 0<R<1.

• Thanks Mariusz. Please how to plot $R$ vs $k$ for any $θ$ in your case ? May 15, 2018 at 21:00
• Mariusz, thank you so much. May 15, 2018 at 21:45
• @Gallagher. You're welcome! :) May 15, 2018 at 21:58
• @Gallagher. FindRootin yours case is unsuitable for your calculations( inadequate function) to solve. For example we can't add assumption like this 0<R<1. May 15, 2018 at 22:40
• @Gallagher.See I edit my answer. May 15, 2018 at 23:02

The division by R in the function eq is superfluous and makes problems in the root finding. Fortunately, some simplification can get you rid of that:

Block[{R, k, θ},
eq = {R, k, θ} \[Function] Evaluate@Simplify[
1 + θ/k^2 + R^2/k^2 (1 + θ (1 + 1/(4 θ) + R^2/4) + (Sqrt[ 2/θ] 1/ R (1 + θ (1 + (R^2)/2) + θ^2/ 2 (1 + (R^2)/2)^2) (-2 DawsonF[ Sqrt[θ/2] R]))),
θ > 0
];
]

Second issue: For some parameter values of k, there may be more than one solution or no solutions at all. NSolve can find (hopefully) all solution:

solr = {k, θ} \[Function] Block[{R}, R /. NSolve[eq[R, k, θ] == 0 && 0 <= R <= 1, R, Reals, WorkingPrecision -> 25]];
data = Block[{θ = 100}, Table[{k, solr[k, θ]}, {k, 1, 3, 1/10}]]
ListPlot[Join @@ Thread /@ data, AxesLabel -> {"k", "R"}]

{{1, {0.1294234646803674672924767}}, {11/10, {0.1296843073043101973816057}}, {6/5, {0.1299709530763769620463354}}, {13/10, {0.1302836799112623098323127}}, {7/5, {0.1306227971663856528622202}}, {3/2, {0.1309886477076770723646643}}, {8/5, {0.1313816102517342127410522}}, {17/10, {0.1318021020175089797099938}}, {9/5, {0.1322505817265513300747067}}, {19/10, {0.1327275529978164798164777}}, {2,{0.1332335681913876276756329}}, {21/ 10, {0.1337692327655012292943036, 0.8928149869232539863495161}}, {11/5, {0.1343352102233942222258231, 0.7840644848720745156904206}}, {23/10, {0.1349322277412426656153432, 0.7061209842390900229330888}}, {12/5, {0.1355610825864964635669894, 0.6468555823332288287855214}}, {5/2, {0.1362226494580973891652682, 0.5999539682600664300172137}}, {13/5, {0.1369178889075194101436787, 0.5617652337478442737838912}}, {27/10, {0.1376478570337623102004883, 0.5300143374213614722559567}}, {14/5, {0.1384137166883060426333669, 0.5031957080178530201387233}}, {29/10, {0.1392167504801858849697920, 0.4802563992569981158981134}}, {3, {0.1400583759402635089617742, 0.4604211571053870525255899}}}

• Thanks Henrik. Please how to plot $R$ vs $k$ for any $θ$? May 15, 2018 at 20:59
• Henrik, thank you so much. May 15, 2018 at 21:45
• You're welcome! May 15, 2018 at 21:47