0
$\begingroup$

What could possibly went wrong in my code? Basically, I am solving the differential equation $\textbf{ode}$ using $\textbf{NDSolve}$. But mathematica says, NDSolve::mxst: Maximum number of 10000 steps reached at the point y == -0.431184926619225550730163530627, and InterpolatingFunction::dmval: Input value {1/1000} lies outside the range of data in the interpolating function. Extrapolation will be used.

What might cause these errors and how to address these? Also, there is an infinite expression $\frac{1}{0}$ encountered in my code. Is there a better way to execute my code?

 rl[\[Rho]_, b0_, q_] := \[Rho] (1 - (Sqrt[\[Pi]]Gamma[1/(q - 1)])/((1 - q) Gamma[1/2 ((q + 1)/(q - 1))])b0/\[Rho])
 r1 = FullSimplify[rl[\[Rho], a, -2]];
 f0 = (1 - a^3/r1^3) /. \[Rho] -> a y0;

 ode = -l (1 + l) R[\[Rho]] + 2 (\[Rho] + (a Sqrt[\[Pi]] Gamma[2/3])/Gamma[1/6]) Derivative[1][][\[Rho]] + (\[Rho] + (a Sqrt[\[Pi]] Gamma[2/3])/Gamma[1/6])^2 (R^\[Prime]\[Prime])[\[Rho]]
 odey = ode /. R -> (R[#/a] &) /. \[Rho] -> a y // Simplify

 R[y_, l_] := (1/y)^(1 + l) - ((2 + 3 l + l^2) \[Pi] (1/y)^(3 + l) Gamma[2/3]^2)/(2 (3 + 2 l) Gamma[1/6]^2) + ((24 + 50 l + 35 l^2 + 10 l^3 + l^4) \[Pi]^2 (1/y)^(5 + l) Gamma[2/3]^4)/(8 (15 + 16 l + 4 l^2) Gamma[1/6]^4); 
 dR[y_, l_] := -(1 + l) (1/y)^(2 + l) + ((6 + 11 l + 6 l^2 + l^3) \[Pi] (1/y)^(4 + l)Gamma[2/3]^2)/(2 (3 + 2 l) Gamma[1/6]^2) - ((120 + 274 l + 225 l^2 + 85 l^3 + 15 l^4 + l^5) \[Pi]^2 (1/y)^(6 + l) Gamma[2/3]^4)/(8 (15 + 16 l + 4 l^2) Gamma[1/6]^4);

 rules = {AccuracyGoal -> Infinity, PrecisionGoal -> 20, WorkingPrecision -> 30, MaxSteps -> 10000};
 rat = 10^-30;
 yP = 10^3;
 yM = -10^3;

 a = 1;
 Q = 1;
 f[r_] := 1 - a^3/r^3;
 \[Psi][r_] := 0;
 aRP = -1/(r^2 Sqrt[f[r]]);
 bRP = -1/(2 r^2) (1 + r \[Psi]'[r]);
 dRP = -1/(16 r^2) ((1 + r \[Psi]'[r]) - (1 + r \[Psi]'[r] + 3 r^2 \[Psi]'[r]^2 + r^3 \[Psi]'[r]^3 - 6 r^2 \[Psi]''[r] - 2 r^3 D[\[Psi][r], {r, 3}]) f[r] + (1 + 4 r \[Psi]'[r] + 3 r^2 \[Psi]''[r]) r f'[r] + (1 + r \[Psi]'[r]) r^2 f''[r]);

 y0 = 10^-3;

 For[el = 0, el <= 5, el++, {
 R0p = Rationalize[R[yP, el], rat];
 dR0p = Rationalize[dR[yP, el], rat];
 R0m = Rationalize[R[yM, el], rat];
 dR0m = Rationalize[dR[yM, el], rat];

 Rsolp = R /. First@NDSolve[{(odey /. l -> el) == 0, R[yP] == R0p, R'[yP] == dR0p}, {R}, {y, yP, y0}, rules];
 Rsolm = R /. First@NDSolve[{(odey /. l -> el) == 0, R[yM] == R0m, R'[yM] == dR0m}, {R}, {y, yM, y0}, rules];

rp = Rationalize[ Rsolp[y0], rat];
rm = Rationalize[ Rsolm[y0], rat];
drp = Rationalize[ Rsolp'[y0], rat];
drm = Rationalize[ Rsolm'[y0], rat];
s = Rationalize[(-Q a Sqrt[4 Pi (2 el + 1)])/(r1 /. \[Rho] -> a y0)^2,rat];

c1f = (rm s)/(drp rm - drm rp);(*jump condition*)
baremode[el] = Sqrt[(2 el + 1)/(4 \[Pi])] (c1f Rsolp'[y0]);}]
$\endgroup$
  • $\begingroup$ DSolve instead of NDSolve works for me after cleaning up your ode. But I don't see that you are saving any results in your For loop. $\endgroup$ – Bill Watts Feb 1 at 1:21
  • $\begingroup$ The for loop should serve as interpolating function I think? The ode that enters in the for loop is $\textbf{odey}$, which is the nondimensionalized form of the ode. If I run the for loop, it should not return any output. But in my case, there are many errors. What did you do to fix it up? I should use NDSolve since I want to get the numerical solutions. $\endgroup$ – user583893 Feb 1 at 1:42
  • $\begingroup$ I assume you really mean R'[p] and R''[p] in your ode. I did not understand Derivative[1][][p] and MMA didn't like it either. You can plug numbers in results returned by DSolve. $\endgroup$ – Bill Watts Feb 1 at 1:53
  • $\begingroup$ can i see your code? $\endgroup$ – user583893 Feb 1 at 1:59
  • $\begingroup$ ah yes, it's the first and second derivative of R. Sorry for that $\endgroup$ – user583893 Feb 1 at 2:03
1
$\begingroup$
Clear["Global`*"]

rl[ρ_, b0_, 
  q_] := ρ (1 - (Sqrt[π] Gamma[1/(q - 1)])/((1 - q) Gamma[
         1/2 ((q + 1)/(q - 1))]) b0/ρ)
r1 = FullSimplify[rl[ρ, a, -2]];
f0 = (1 - a^3/r1^3) /. ρ -> a y0;

ode = -l (1 + l) R[ρ] + 
   2 (ρ + (a Sqrt[π] Gamma[2/3])/
      Gamma[1/6]) R'[ρ] + (ρ + (a Sqrt[π] Gamma[2/3])/
      Gamma[1/6])^2 R''[ρ];

odey = ode /. R -> (R[#/a] &) /. ρ -> a y // Simplify;

R[y_, l_] := (1/y)^(1 + 
      l) - ((2 + 3 l + l^2) π (1/y)^(3 + l) Gamma[2/3]^2)/(2 (3 + 
        2 l) Gamma[1/6]^2) + ((24 + 50 l + 35 l^2 + 10 l^3 + 
        l^4) π^2 (1/y)^(5 + l) Gamma[2/3]^4)/(8 (15 + 16 l + 
        4 l^2) Gamma[1/6]^4);
dR[y_, l_] := -(1 + l) (1/y)^(2 + 
       l) + ((6 + 11 l + 6 l^2 + l^3) π (1/y)^(4 + l) Gamma[
        2/3]^2)/(2 (3 + 2 l) Gamma[1/6]^2) - ((120 + 274 l + 
        225 l^2 + 85 l^3 + 15 l^4 + l^5) π^2 (1/y)^(6 + l) Gamma[
        2/3]^4)/(8 (15 + 16 l + 4 l^2) Gamma[1/6]^4);

rules = {AccuracyGoal -> Infinity, PrecisionGoal -> 20, 
   WorkingPrecision -> 30, MaxSteps -> 10000};
rat = 10^-30;
yP = 10^3;
yM = -10^3;

a = 1;
Q = 1;
f[r_] := 1 - a^3/r^3;
ψ[r_] := 0;
aRP = -1/(r^2 Sqrt[f[r]]);
bRP = -1/(2 r^2) (1 + r ψ'[r]);
dRP = -1/(16 r^2) ((1 + 
       r ψ'[r]) - (1 + r ψ'[r] + 3 r^2 ψ'[r]^2 + 
        r^3 ψ'[r]^3 - 6 r^2 ψ''[r] - 
        2 r^3 D[ψ[r], {r, 3}]) f[
       r] + (1 + 4 r ψ'[r] + 3 r^2 ψ''[r]) r f'[
       r] + (1 + r ψ'[r]) r^2 f''[r]);

y0 = 10^-3;

el = 0

dR0p = Rationalize[dR[yP, el], rat];
R0m = Rationalize[R[yM, el], rat];
dR0m = Rationalize[dR[yM, el], rat];

Rsolp[y_] = 
  R[y] /. DSolve[{(odey /. l -> el) == 0, R[yP] == R0p, 
       R'[yP] == dR0p}, R[y], y][[1]] // Simplify;

Rsolp[y] // N // Simplify

(*(1.00086 - 4.31061*10^-7 y)/(1. y + 0.431185)*)

Rsolm[y_] = 
  R[y] /. DSolve[{(odey /. l -> el) == 0, R[yM] == R0m, 
       R'[yM] == dR0m}, R[y], y][[1]] // Simplify;

Rsolm[y] // N // Simplify

(*(0.999137 - 4.31309*10^-7 y)/(1. y + 0.431185)*)

The you can set a new el and do it again.

$\endgroup$

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.