0
$\begingroup$

I'm trying to solve a pair of uncoupled differential equations using NDSolve, where the initial conditions are determined by FindRoot. Here's the relevant section of the code, including error message. As you can see, the result evaluates perfectly well pointwise, but fails to plot (I suspect due to some issue with hold and the order in which things are evaluated).

V0 = 10^-10;
\[Mu] = 1;
\[Omega] = 2;
amp = 0.1;
freq = 0.1;

V[\[Phi]_] := 
  V0*(Tanh[\[Phi]/\[Omega]] + amp*Sin[Tanh[\[Phi]/\[Omega]]/freq])^2;

aini = 1;
\[Phi]ini = 5;
d\[Phi]ini = -(V'[\[Phi]ini]/V[\[Phi]ini]);

sol = NDSolve[{dummy\[Phi]''[n] + 3*dummy\[Phi]'[n] - 
       1/2*(dummy\[Phi]'[n])^3 + (3 - 1/2*(dummy\[Phi]'[n])^2)*
        V'[dummy\[Phi][n]]/V[dummy\[Phi][n]] == 0, 
     dummy\[Phi][0] == \[Phi]ini, dummy\[Phi]'[0] == d\[Phi]ini}, 
    dummy\[Phi], {n, 0, 56}, MaxStepSize -> Infinity][[1, 1]];

\[Phi][n_] := dummy\[Phi][n] /. sol;
d\[Phi][n_] := Derivative[1][dummy\[Phi]][n] /. sol;
dd\[Phi][n_] := (dummy\[Phi]^\[Prime]\[Prime])[n] /. sol;

\[Epsilon]1[n_] := 1/2*d\[Phi][n]^2;
\[Epsilon]2[n_] := -(1/\[Epsilon]1[n])*d\[Phi][n]*dd\[Phi][n];
\[Epsilon]3[n_] := 
  1/(\[Epsilon]1[n]*\[Epsilon]2[n])*(dd\[Phi][n]^2 + 
      d\[Phi][n]*ddd\[Phi][n]) - \[Epsilon]2[n];

H[n_] := Sqrt[V[\[Phi][n]]/(3 - \[Epsilon]1[n])];
ascale[n_] := aini*Exp[n];

Nini[k_?NumericQ] := 
  ni /. FindRoot[
    k == 10^2*H[ni]*ascale[ni] // Rationalize[#, 0] & // 
     Evaluate, {ni, 0, 56}, WorkingPrecision -> 20];
Nfin[k_?NumericQ] := 
  nf /. FindRoot[
    k == 10^-2*H[nf]*ascale[nf] // Rationalize[#, 0] & // 
     Evaluate, {nf, 0, 56}, WorkingPrecision -> 20];

ddzoverz[n_] := 
  2 - \[Epsilon]1[n] + 3/2 \[Epsilon]2[n] + 1/4 \[Epsilon]2[n]^2 - 
   1/2 \[Epsilon]1[n]*\[Epsilon]2[n] + 
   1/2 \[Epsilon]2[n]*\[Epsilon]3[n];

eqs[k_, n_] := {(vR^\[Prime]\[Prime])[n] + (1 - \[Epsilon]1[n])*
      Derivative[1][vR][n] + (k^2/(ascale[n]^2*H[n]^2) - ddzoverz[n])*
      vR[n] == 
    0, (vI^\[Prime]\[Prime])[n] + (1 - \[Epsilon]1[n])*
      Derivative[1][vI][n] + (k^2/(ascale[n]^2*H[n]^2) - ddzoverz[n])*
      vI[n] == 0, vR[Nini[k]] == 1/Sqrt[2 k], 
   Derivative[1][vR][Nini[k]] == 0, vI[Nini[k]] == 0, 
   Derivative[1][vI][Nini[k]] == -50 Sqrt[2/k]};

vRsol[k_] := 
  NDSolve[eqs[k, n], {vR, vI}, {n, Nini[k], Nfin[k]}, 
    MaxStepSize -> Infinity][[1, 1]];
vIsol[k_] := 
  NDSolve[eqs[k, n], {vR, vI}, {n, Nini[k], Nfin[k]}, 
    MaxStepSize -> Infinity][[1, 2]];

vRk[k_?NumericQ, n_?NumericQ] := vR[n] /. vRsol[k];
vIk[k_?NumericQ, n_?NumericQ] := vI[n] /. vIsol[k];

vRk[10^11, 33]
-7.88251*10^-7
vRk[10^11, 34]
1.84531*10^-6
vRk[10^11, 35]
2.05082*10^-6

Plot[vRk[10^11, n], {n, 33, 37}]
NDSolve::dsvar: 33.00008171428572` cannot be used as a variable.

ReplaceAll::reps: {7307.25 vR[33.0001]+0.999996 (vR^\[Prime])[33.0001]+(vR^\[Prime]\[Prime])[33.0001]==0} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing.

ReplaceAll::reps: {7307.25 vR[33.0001]+0.999996 (vR^\[Prime])[33.0001]+(vR^\[Prime]\[Prime])[33.0001]==0.} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing.

NDSolve::dsvar: 33.08171436734695` cannot be used as a variable.

ReplaceAll::reps: {6206.23 vR[33.0817]+0.999996 (vR^\[Prime])[33.0817]+(vR^\[Prime]\[Prime])[33.0817]==0} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing.

General::stop: Further output of ReplaceAll::reps will be suppressed during this calculation.

NDSolve::dsvar: 33.163347020408175` cannot be used as a variable.

General::stop: Further output of NDSolve::dsvar will be suppressed during this calculation.

I'd generally define "eqs" with set rather than making it a function and using setdelay, but when I tried it this way I got the error message:

"FindRoot: The function value {-0.000544959+k} is not a list of numbers with dimensions {1} at {ni} = {0}"

I suspect my issue is the fact that everything is a function of the parameter k, so I tried a variant on this using ParametricNDSolveValue, but because k appears as part of the initial conditions, I got an error that a starting point for the variable vR couldn't be found.

Any help or workaround would be greatly appreciated!!

$\endgroup$
7
  • $\begingroup$ "ascale"and "H" and "epsiloni" are not defined. $\endgroup$ May 6, 2022 at 8:41
  • $\begingroup$ Sorry--I've edited the question to include all definitions now. $\endgroup$
    – anireland
    May 6, 2022 at 14:13
  • $\begingroup$ Derivatives are not written; "vR^[Prime][Prime][n]", but "vR''[n]". Or the long form: "Derivative[2][vR][n]" $\endgroup$ May 6, 2022 at 15:22
  • $\begingroup$ I'm not sure what's up with the formatting here, but in the actual notebook they are written as "vR''[n]" $\endgroup$
    – anireland
    May 6, 2022 at 19:48
  • $\begingroup$ \[Prime] has character code 8242, however the character used in "vR''[n]" is a single quote with character code 39 $\endgroup$ May 6, 2022 at 20:07

1 Answer 1

1
$\begingroup$

Just in case anyone runs into a similar issue, I figured out a fix! The issue was that NDSolve was setting n to a numeric value, and then trying to evaluate. So I used Block so that it would be evaluated with local n. Here is the last segment of my code after having implemented this:

vRsol[k_] := 
  Block[{n}, 
    NDSolve[eqs[k, n], {vR, vI}, {n, Nini[k], Nfin[k]}, 
     MaxStepSize -> Infinity]][[1, 1]];
vIsol[k_] := 
  Block[{n}, 
    NDSolve[eqs[k, n], {vR, vI}, {n, Nini[k], Nfin[k]}, 
     MaxStepSize -> Infinity]][[1, 2]];

vRk[k_?NumericQ, n_?NumericQ] := vR[n] /. vRsol[k];
vIk[k_?NumericQ, n_?NumericQ] := vI[n] /. vIsol[k];

Table[vRk[10^11, n], {n, 33, 37, 1}]
{-7.88037*10^-7, 1.84941*10^-6, 2.0586*10^-6, 
 6.50416*10^-7, -2.31497*10^-6}
$\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.