# How to get accurate solutions from FindRoot?

Gallagher 05/15/2018. 2 answers, 255 views

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!

### 2 Answers

Mariusz Iwaniuk 05/15/2018.
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.

Henrik Schumacher 05/16/2018.

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}}}