1
$\begingroup$

I'm trying to solve the system

Subscript[m, p] = 1.22*10^(19) * 1.52*10^(24);
Ag = Derivative[1][b][t]-1.5608310063810777811971842694*10^-43 b[t] Sqrt[(1 + 2.9079917009173806522983649735*10^-87 \[Psi][t]^2) (1.18253*10^149 f^4 (1 - (2.90799*10^-87 \[Psi][t]^2)/h^2)^2 + 1/2 Derivative[1][\[Psi]][t]^2)]

Bg = -((1.37552*10^63 f^4 \[Psi][t] (1 - (2.90799*10^-87 \[Psi][t]^2)/h^2))/h^2) - (3 \[Psi][t] Derivative[1][b][t]^2)/(4 b[t]^2 (1 + 2.9079917009173806522983649735*10^-87 \[Psi][t]^2)^2) + (3 Derivative[1][b][t] Derivative[1][\[Psi]][t])/ b[t] + Derivative[2][\[Psi]][t]

Using

sol = ParametricNDSolve[{Bg == 0, Ag == 0, \[Psi][0] == 2.1*(Subscript[m, p]),b[0] == 1, \[Psi]'[0] == 0}, {\[Psi], b}, {t, 0.0, 2.2*10^(-31)}, {f,h}, MaxSteps -> 100000, Method -> {"ExplicitRungeKutta", "DifferenceOrder" -> 4,"Coefficients" -> FehlbergCoefficients,"StiffnessTest" -> False}];
Manipulate[ParametricPlot[{\[Psi][f, h][t], b[f, h][t]} /. sol, {t, 0,6*10^(-32)}, AspectRatio -> 1], {f, 0.1, 10}, {h, 0.1, 10}]

Where the FehlbergCoefficients is:

Fehlbergamat = {{1/4}, {3/32, 9/32}, {1932/2197, -7200/2197, 
7296/2197}, {439/216, -8, 3680/513, -845/4104}, {-8/27, 
2, -3544/2565, 1859/4104, -11/40}};
Fehlbergbvec = {25/216, 0, 1408/2565, 2197/4104, -1/5, 0};
Fehlbergcvec = {1/4, 3/8, 12/13, 1, 1/2};
Fehlbergevec = {-1/360, 0, 128/4275, 2197/75240, -1/50, -2/55};
FehlbergCoefficients[4, p_] := 
N[{Fehlbergamat, Fehlbergbvec, Fehlbergcvec, Fehlbergevec}, p];

But it fails with $f=5$ and $h=1$, (I got "\$abort"). Using $f = 5$ and $h = 1$ before integration with jut "NDsolve" the result is what I'm looking for.

What I need is a plot $\psi \times t$ and $b \times t$ with variable coefficients h and f.

$\endgroup$
  • 1
    $\begingroup$ The quantity Subscript[m, p] is undefined. In general, it is wise not to use superscript or subscript variables in computations. $\endgroup$ – bbgodfrey Oct 11 '17 at 16:32
  • $\begingroup$ Sorry: Subscript[m, p] = 1.22*10^(19) * 1.52*10^(24) $\endgroup$ – Kamog Oct 11 '17 at 16:53
  • $\begingroup$ With v 11.2, I cannot reproduce the error you describe. InterpolatingFunction does complain that some numerical values of t lie outside its range, but that is not serious and can be fixed. $\endgroup$ – bbgodfrey Oct 11 '17 at 17:17
  • $\begingroup$ @bbgodfrey, it is wierd, setting $f = 5$ and $h=1$ before integration the numerical values of $t$ lies inside its range. My version is 10. $\endgroup$ – Kamog Oct 11 '17 at 17:28
2
$\begingroup$

The computation as described in the question fails, because FehlbergCoefficients is undefined. It also seems unnecessary. To produce the desired results, replace the second block of code by

sol = ParametricNDSolve[{Bg == 0, Ag == 0, ψ[0] == 2.1*(Subscript[m, p]), 
    b[0] == 1, ψ'[0] == 0}, {ψ, b}, {t, 0.0, 2.2*10^(-31)}, {f, h}];
Manipulate[ParametricPlot[{ψ[f, h][t], b[f, h][t]} /. sol, {t, 0, 
    6*10^(-32)}, AspectRatio -> 1], {{f, 5}, 0.1, 10}, {{h, 1}, 0.1, 10}]

enter image description here

If it is necessary for some reason to employ FehlbergCoefficients, use the code provided under "A Method of Fehlberg" here.

Addendum

Defining FehlbergCoefficients, as just added to the question, allows

sol = ParametricNDSolve[{Bg == 0, Ag == 0, ψ[0] == 2.1*(Subscript[m, p]), 
    b[0] == 1, ψ'[0] == 0}, {ψ, b}, {t, 0.0, 2.2*10^(-31)}, {f, h}, 
    Method -> {"ExplicitRungeKutta", "DifferenceOrder" -> 4, 
               "Coefficients" -> FehlbergCoefficients, "StiffnessTest" -> False}];

to be used instead. It gives the same result.

$\endgroup$
  • $\begingroup$ My FehlbergCoefficients is defined... this plot is for $\psi$ or $b$ ? $\endgroup$ – Kamog Oct 12 '17 at 1:23
  • $\begingroup$ @Kamog Where in your question is FehlbergCoefficients defined? The plot is of b as a function of ψ, as you specified by using ParametricPlot. If, instead, you want b and ψ as functions of t, use Plot instead. $\endgroup$ – bbgodfrey Oct 12 '17 at 1:29
  • $\begingroup$ Sorry, I forgot to put it in my question, but I did it now.... thank you $\endgroup$ – Kamog Oct 12 '17 at 1:31
  • $\begingroup$ Using: Manipulate[ Plot[{[Psi][f, h][t]} /. sol, {t, 0, 6*10^(-32)}, AspectRatio -> 1], {f, 0.1, 10}, {h, 0.1, 10}] I got the same error... =\ $\endgroup$ – Kamog Oct 12 '17 at 1:41
  • $\begingroup$ I cannot reproduced the error. I suggest that you copy the definitions from your question plus the code in my answer and run them in a new notebook. $\endgroup$ – bbgodfrey Oct 12 '17 at 1:52

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.