6
$\begingroup$

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!

$\endgroup$
9
  • 2
    $\begingroup$ 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. $\endgroup$
    – bbgodfrey
    Commented May 15, 2018 at 19:39
  • $\begingroup$ 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! $\endgroup$
    – Chris K
    Commented May 16, 2018 at 1:04
  • $\begingroup$ Of course. Thank you. $\endgroup$
    – Gallagher
    Commented May 16, 2018 at 5:21
  • 1
    $\begingroup$ 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. $\endgroup$
    – bbgodfrey
    Commented May 21, 2018 at 23:14
  • 1
    $\begingroup$ Reduce[(eq1 /. {k -> 3, θ -> 100}) == 0 && 0 < R < 1, {R}] // N also works and may use the same algorithm in this case. $\endgroup$
    – bbgodfrey
    Commented May 21, 2018 at 23:28

2 Answers 2

7
$\begingroup$
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

enter image description here

For: n=15

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

enter image description here

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

enter image description here

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

enter image description here


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

enter image description here

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

$\endgroup$
8
  • $\begingroup$ Thanks Mariusz. Please how to plot $R$ vs $k$ for any $θ$ in your case ? $\endgroup$
    – Gallagher
    Commented May 15, 2018 at 21:00
  • $\begingroup$ Mariusz, thank you so much. $\endgroup$
    – Gallagher
    Commented May 15, 2018 at 21:45
  • $\begingroup$ @Gallagher. You're welcome! :) $\endgroup$ Commented May 15, 2018 at 21:58
  • 1
    $\begingroup$ @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. $\endgroup$ Commented May 15, 2018 at 22:40
  • 1
    $\begingroup$ @Gallagher.See I edit my answer. $\endgroup$ Commented May 15, 2018 at 23:02
6
$\begingroup$

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

enter image description here

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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