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!

|improve this question
$\endgroup$
  • 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 May 15 '18 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 May 16 '18 at 1:04
  • $\begingroup$ Of course. Thank you. $\endgroup$ – Gallagher May 16 '18 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 May 21 '18 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 May 21 '18 at 23:28
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.

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

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

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.