2
$\begingroup$

I am trying to apply @bbgodfrey method for finding the right initial guess for the following problem but I get an error message and I don't understand where it comes from.

I would also like to know why FindRoot in Do loop doesn't work in this case.

Below is my minimal working code:

l1 = 0.81;
Z = 500;
x0 = 10;
v0 = 0.02;
\[Epsilon] = $MachineEpsilon;
yl = -12;
yu = 0;

l0 = 0.0714`20.;

ps = ParametricNDSolveValue[{y''[r] + 
      2 y'[r]/r == -4 \[Pi] l k Exp[-y[r]], y[\[Epsilon]] == y0, 
    y'[\[Epsilon]] == 0, WhenEvent[r == 1, y'[r] -> y'[r] + Z l]}, {y,
     y'}, {r, \[Epsilon], R}, {k, l}, 
   Method -> {"StiffnessSwitching"}, AccuracyGoal -> 5, 
   PrecisionGoal -> 4, WorkingPrecision -> 15];

Do[x = i x0;
  v = i^3 v0;
  yl = -12;
  yu = 0;
  R = Rationalize[v^(-1/3), 0];
  l = Rationalize[l1/(i x0), 0];
  fy := (Do[yc = (yl + yu)/2;  (* guess finder function *)

     test = First[ps[yc]]["Domain"][[1, 2]];
     If[test == 1, Throw[yc]];
     If[Last[ps[yc]][test] > 0, yu = yc, yl = yc], {i, 50}]; yc);
  yint = 
   Which[1 == First[ps[yl]]["Domain"][[1, 2]], yl, 
    1 == First[ps[yu]]["Domain"][[1, 2]], yu, True, Catch[fy]];
  nn = FindRoot[Last[ps[y0]][R], {y0, yint}, Evaluated -> False][[1, 
     2]];
  Tot = 4 \[Pi] nn NIntegrate[
     r^2 Exp[-First[ps[nn, l]][r]], {r, \[Epsilon], R}, 
     PrecisionGoal -> 4];
  Print[NumberForm[i*1., 5], "  ", NumberForm[Tot, 5]];, {i, 292/100, 
   31/10, 1/100}]
Improve this question
$\endgroup$
17
  • 1
    $\begingroup$ The first error message your code gives me says that ps was called with one parameter having a value of -12. That would be in the line ψint = ... ps[ yl ] ... Your definition of ps expects 2 parameters. However, the definition of ps contains y0 and R, which are undefined. Perhaps y0 and R should be a third and fourth parameters of ps ? If they are not parameters, they need numeric values. Also, to debug loops, start with With[ {i=3}, ... ] and add statements one at a time, making sure each one works for the value of i. Then try to replace the With with Table instead of Do. $\endgroup$
    – LouisB
    Commented Jan 25, 2020 at 3:25
  • $\begingroup$ @LouisB First, I need to work with Do. And in Do statements R is numerical value. For y0, I belived that the code should take yint and pass it to y0. Also, there is a typo ψint should be yint. $\endgroup$
    – aluuzz
    Commented Jan 25, 2020 at 3:37
  • 1
    $\begingroup$ Just to understand the syntax, try this immediately after your ps = statement, y0 = 1; R = 1.2; Plot[ ps[0.5, 0.027] [[1]] [r], {r, 0, R}, PlotRange -> {{0, R}, {0, 4}}]. This indicates we must set values of y0 and R, then evaluate ps with 2 parameters, followed by a subscript of 1 or 2 followed by [ r ]. Of course we could use First and Last instead of the [[ ... ]] notation. $\endgroup$
    – LouisB
    Commented Jan 25, 2020 at 5:00
  • $\begingroup$ @LouisB OK. I just want to let know that the code that I am using is taken from here [mathematica.stackexchange.com/questions/183148/…, and it works fine. What I am asking here is to implemented bbgodfrey method for finding the right initial guess to the above code. $\endgroup$
    – aluuzz
    Commented Jan 25, 2020 at 5:12
  • $\begingroup$ @aluuzz 1) what problem do you want to solve? 2) it looks so that you have syntax errors there. $\endgroup$
    – Alex Trounev
    Commented Jan 28, 2020 at 20:53

1 Answer 1

Reset to default
1
$\begingroup$

We can use my code without changing for Z<=15000

l1 = 0.81;
Z = 15000;
x0 = 10;
v0 = 0.02;
\[Epsilon] = $MachineEpsilon;

l0 = 0.0714`20.;

ps = ParametricNDSolveValue[{y''[r] + 
      2 y'[r]/r == -4 \[Pi] l k Exp[-y[r]], y[\[Epsilon]] == y0, 
    y'[\[Epsilon]] == 0, WhenEvent[r == 1, y'[r] -> y'[r] + Z l]}, {y,
     y'}, {r, \[Epsilon], R}, {k, l}, 
   Method -> {"StiffnessSwitching"}, AccuracyGoal -> 5, 
   PrecisionGoal -> 4, WorkingPrecision -> 15];

Do[x = i x0;
  v = i^3 v0;
  R = Rationalize[v^(-1/3), 0];
  l = Rationalize[l1/(i x0), 0];
  nn = FindRoot[Last[ps[y0, l]][R], {y0, -1}, Evaluated -> False][[1, 
    2]];
  Tot = 4 \[Pi] nn NIntegrate[
     r^2 Exp[-First[ps[nn, l]][r]], {r, \[Epsilon], R}, 
     PrecisionGoal -> 4];
  Print[NumberForm[i*1., 5], "  ", NumberForm[Tot, 5]];, {i, 292/100, 
   31/10, 1/100}] // Quiet

Out: 

2.92  15007.

2.93  15007.

2.94  15006.

2.95  15006.

2.96  15006.

2.97  15006.

2.98  15006.

2.99  15006.

3.  15006.

3.01  15006.

3.02  15006.

3.03  15006.

3.04  15006.

3.05  15006.

3.06  15006.

3.07  15006.

3.08  15005.

3.09  15005.

3.1  15005.

I checked code @bbgodfrey (question as it is). The result was much worse than with my code, but the code works for some k up to Z = 20000:

p[Z0_, g0_, k0_, R0_] := 
 Block[{Z = Z0, g = Rationalize[g0, 0], k2 = Rationalize[k0, 0], 
   yl = -8, yu = 0, ps, fy, y00, sol}, \[Epsilon] = 
   10^-4; R = Rationalize[R0, 0]; 
  ps = ParametricNDSolveValue[{y''[r] + 
       2 y'[r]/r == -4 Pi k2 Exp[-y[r]], y[\[Epsilon]] == y0, 
     y'[\[Epsilon]] == 0, 
     WhenEvent[r == 1, y'[r] -> y'[r] + Z g]}, {y, 
     y'}, {r, \[Epsilon], R}, {y0}, Method -> "StiffnessSwitching", 
    WorkingPrecision -> 20];
  fy := (Do[yc = (yl + yu)/2;
     tst = First[ps[yc]]["Domain"][[1, 2]];
     If[tst == R, Throw[yc]];
     If[Last[ps[yc]][tst] > 0, yu = yc, yl = yc], {i, 50}]; yc);
  y00 = Which[R == First[ps[yl]]["Domain"][[1, 2]], yl, 
    R == First[ps[yu]]["Domain"][[1, 2]], yu, True, Catch[fy]];
  sol = FindRoot[Last[ps[y0]][R], {y0, y00}, Evaluated -> False][[1, 
    2]]; L = 
   4 \[Pi] k2/g0 NIntegrate[
     r^2 Exp[-First[ps[sol]][r]], {r, \[Epsilon], R}]]

Using of p[] 

l1 = 0.81;
Z0 = 20000;
x0 = 10;
v0 = 0.02;
l0 = 0.0714`20.; k = 1/29;
 Do[x = i x0;
 v = i^3 v0;
 R = Rationalize[v^(-1/3), 0];
 l = Rationalize[l1/(i x0), 0];
 Print[NumberForm[i*1., 5], "  ", 
  NumberForm[p[Z0, l, k l, R] // Quiet, 5]];, {i, 292/100, 31/10, 
  1/100}]

Out

2.92  20010.

2.93  20010.

2.94  20010.

2.95  20009.

2.96  20009.

2.97  20009.

2.98  20009.

2.99  20008.

3.  20008.

3.01  8.3397*10^(21)

3.02  1.309*10^(22)

3.03  2.0575*10^(22)

3.04  3.2386*10^(22)

3.05  5.1048*10^(22)

3.06  8.0576*10^(22)

3.07  1.2736*10^(23)

3.08  2.0159*10^(23)

3.09  3.1952*10^(23)

3.1  5.0714*10^(23)
Improve this answer
$\endgroup$

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.