0
$\begingroup$

My code is

(*constants*)
c = 1; 
G = 1; 
R0 = 1.473; 
k = 4.012/10^4; 
γ = 2; 
α = 0.2; 
ϵ0 = (((1/k)*(R0/α)^γ)^(1 - γ))^(-((1 - γ)^2)^(-1)); 
ϵ0s = 1.603*10^38; 
ScaleMeVfmtoMs = 8.92/10^7; 
ScaleErgcmtoMs = 5.56/10^40; 
β = ((4*Pi*ϵ0)/(k*ϵ0^(γ - 1))^(1/γ))*ScaleMeVfmtoMs; 
(*Solving*)
TOV = Derivative[1][p][r] + (((α*M[r])/r^2)*p[r]^(1/γ)*(1 + p[r]/(p[r]/(k*ϵ0^(γ - 1)))^(1/γ))*(1 + (4*Pi*r^3*p[r])/(M[r]*c^2)))/
     (1 - (2*G*M[r])/(c^2*r)); 
rTOV[p0_] := NDSolve[{TOV == 0, Derivative[1][M][r] == β*r^2*p[r]^(1/γ), p[0.001] == p0, 
     M[0.001] == (4/3)*Pi*0.001^3*((p[0.001]/(k*ϵ0^(γ - 1)))^(1/γ)/c^2)}, {p, M}, {r, 0.001, 30}, 
    Method -> {"FixedStep", Method -> "ExplicitRungeKutta"}]; 
pTOV[r_, p0_] := p[r] /. rTOV[p0]
MTOV[r_, p0_] := M[r] /. rTOV[p0]
NR[p0_] := FindRoot[pTOV[r, p0] == 0, {r, 10}][[1,2]]
NM[p0_] := MTOV[NR[p0], p0][[1]]

{NR[10^(-4)], NM[10^(-4)]}

.

The errors message are

FindRoot::nlnum: The function value {InterpolatingFunction[{{0.001,30.}},{5,15,1,{2013207},{4},0,0,0,0,Automatic,{},{},False},{{<<1>>}},{Developer`PackedArrayForm,{<<1>>},{0.0001 +0. I,-1.4911*10^-8+0. I,0.0001 +0. I,-1.46601*10^-8+0. I,0.0001 +0. I,-1.44227*10^-8+0. I,0.0001 +0. I,<<37>>,0.0001 +0. I,-1.163*10^-8+0. I,0.0001 +0. I,-1.15577*10^-8+0. I,0.0001 +0. I,-1.14898*10^-8+0. I,<<4026364>>}},{Automatic}][<<1>>]} is not a list of numbers with dimensions {1} at {r} = {18.569 -6.99405*10^-7 I}.
FindRoot::nlnum: The function value {InterpolatingFunction[{{0.001,30.}},{5,15,1,{2013207},{4},0,0,0,0,Automatic,{},{},False},{{<<1>>}},{Developer`PackedArrayForm,{<<1>>},{0.0001 +0. I,-1.4911*10^-8+0. I,0.0001 +0. I,-1.46601*10^-8+0. I,0.0001 +0. I,-1.44227*10^-8+0. I,0.0001 +0. I,<<37>>,0.0001 +0. I,-1.163*10^-8+0. I,0.0001 +0. I,-1.15577*10^-8+0. I,0.0001 +0. I,-1.14898*10^-8+0. I,<<4026364>>}},{Automatic}][<<1>>]} is not a list of numbers with dimensions {1} at {r} = {18.569 -6.99405*10^-7 I}.

This seems like some steps in FindRoot return a small imaginary value for next step and my function contains only real numbers. Is there any way to prevent FindRoot going to give this small number?

Also, the results are

{18.569 + 0. I, 0.903448 + 1.86409*10^-14 I}
$\endgroup$
  • $\begingroup$ Your ODE has p[r]^(1/γ) where γ=2, so it's no surprise that complex numbers arise when you are trying to find when p[r] is zero. Probably using something like Abs[p[r]]^(1/γ) might help. $\endgroup$ – Carl Woll Jun 12 at 16:43
  • $\begingroup$ Hi Carl. I've edited as your suggestion. The code runs fine without any errors. But I'm not sure if doing this will affect the result? $\endgroup$ – Panithi Nakkhruea Jun 12 at 16:56
  • $\begingroup$ @PanithiNakkhruea You can use NR[p0_] := First[NSolve[pTOV[r, p0] == 0, r]][[1, 2]] // Quiet $\endgroup$ – Alex Trounev Jun 13 at 4:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.