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$
4
  • 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, 2017 at 16:32
  • $\begingroup$ Sorry: Subscript[m, p] = 1.22*10^(19) * 1.52*10^(24) $\endgroup$
    – Kamog
    Oct 11, 2017 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, 2017 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, 2017 at 17:28

1 Answer 1

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$
7
  • $\begingroup$ My FehlbergCoefficients is defined... this plot is for $\psi$ or $b$ ? $\endgroup$
    – Kamog
    Oct 12, 2017 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, 2017 at 1:29
  • $\begingroup$ Sorry, I forgot to put it in my question, but I did it now.... thank you $\endgroup$
    – Kamog
    Oct 12, 2017 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, 2017 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, 2017 at 1:52

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.