Mathematica NDSolvederlen: - Mathematica error - TOV equations

[tau1777]

Hi All,

I'm trying to solve a set of three coupled ODES. I have given initial values for all functions, and I believe I have entered in the code correctly but I still get this error telling me that the derivative operator is not the same length as the number of arguments. What does that mean? And how I fix it? Any and all help is greatly appreciated.

I've included all of my code below. Along with some comments in the lines I think need it. Previously mathematica was giving me an error about the stepsize at r =1, so I just added the option for the first step.

Oh, in case it helps, I'm trying to model a relativistic star. By solving the TOV equations.

In[1]:= n = 1 (*polytropic index*)

Out[1]= 1

In[2]:= \[CapitalGamma] = 1 + 1/n

Out[2]= 2

In[3]:= K1 = 5.3802*10^9 (*polytropic contant in cgs units*)

Out[3]= 5.3802*10^9

In[4]:= Clear[P]

In[5]:= \[Rho]c = 1*10^(15) (*central density of the star, i guessed*)

Out[5]= 1000000000000000

In[6]:= pc = K1 * \[Rho]c^(\[CapitalGamma]) (*central pressure*)

Out[6]= 5.3802*10^39

In[7]:= \[Rho][r_] = (P[r]/K1)^(1/\[CapitalGamma]) (*polytropic equation. different than it usually look because I'm solving it for ρ instead of P. Because I will substitute in ρ for P in the equations below*)

Out[7]= 0.0000136333 Sqrt[P[r]]

In[8]:= M0 = 1.98892*10^33 (*one solar mass in grams*)

Out[8]= 1.98892*10^33

In[9]:= M = 1.4*M0 (*mass of my star*)

Out[9]= 2.78449*10^33

sol = NDSolve[ {D[\[Nu][r], r] == 2*(m[r] + 4*\[Pi]*r^3*P[r])/(r (r - 2*M)),
D[P[r], r] == -(\[Rho][r] + P[r])* ((m[r] + 4*\[Pi]*r^3*P[r])/(r (r - 2*M))),
D[m[r], r] == 4*\[Pi]*r^2*\[Rho][r],
P[1] == pc, m[1] == 0, \[Nu][1] == 11}, {\[Nu], P, m}, {r, 1, 12000},
StartingStepSize → 10]

During evaluation of In[13]:= NDSolve::derlen: The length of the derivative operator Derivative[1] in (\[Nu]^\[Prime])[r] is not the same as the number of arguments. >>

 

[phyzguy]

Your differential equations include 4 functions of r, nu(r), P(r), rho(r), and m(r), but your list of functions in the NDSolve statement only involves 3 of these, nu, P, and m. This is what it is complaining about.

 

[tau1777]

Thanks phyzguy,

I had defined ρ(r) in terms of P(r) , and was hoping that Mathematica would pull that definition but it doesn't. I've tried entering this now:

sol = NDSolve[ {D[\[Nu][r], r] ==
2*(m[r] + 4*\[Pi]*r^3*P[r])/(r (r - 2*M)),
D[P[r], r] == -(0.000013633293674139719` Sqrt[P[r]] +
P[r])* ((m[r] + 4*\[Pi]*r^3*P[r])/(r (r - 2*M))),
D[m[r], r] == 4*\[Pi]*r^2*0.000013633293674139719` Sqrt[P[r]],
P[1] == pc, m[1] == 0, \[Nu][1] == 11}, {\[Nu], P, m}, {r, 1,
12000}, StartingStepSize -> 10]

but I get the error " NDSolve::ndcf: Repeated convergence test failure at r == 1.`; unable to continue."

Any thoughts? Thanks again.

 

[phyzguy]

I think this means that it is unable to find a self-consistent set of values at the first iteration. Are you sure the equations you've entered are self-consistent? Try calculating manually the values of the variables and their derivatives at the first 1 or 2 points.

 

[tau1777]

PLEASE HELP. Re: NDSolve::derlen: - Mathematica error - TOV equations

PLEASE HELP. Ok, I'm about to lose my mind trying to do this Mathematica, I was having trouble implementing my code to solve the TOV equations. So I searched on the internet and found this very helpful pdf here : http://ipnweb.in2p3.fr/~PhT-IPN/membres/docs/Margueron-lecture1-Coimbra.pdf

I wrote the code out exactly the way they have it in page 11 I believe, titled "Solving the TOV equations with Mathematica". However when I go to solve I get the same error about the length of the derivative operator. I went through it and saw that just as I was defining my density function, they were doing that as well. So replaced the density function inside of NDSolve with a function explicitly of pressure. And STILL I get this error about the length of the derivative operator.

If you can find the problem please tell me. Here's my code:

(*constants in units of CGS*)

msol = 1.989*10^(-33) (*1 solar mass*)

Out[1]= 1.989*10^-33

G = 6.67*10^(-8) (*gravitational constant*)

Out[2]= 6.67*10^-8

c = 2.9979*10^(-10) (*speed of light*)

Out[21]= 2.9979*10^-10

rsol = 2*G*msol/c^2 (*Schwarschild radius of the Sun*)

Out[4]= 2.95227*10^-21

rhosol = msol*3/(4*Pi*rsol^3) (*density of the Sun*)

Out[5]= 1.84534*10^28

In[6]:=
(*define polytropic equation and its inverse*)

In[7]:= K1 = 5.38*10^(9)

Out[7]= 5.38*10^9

In[8]:= n1 = 1

Out[8]= 1

In[9]:= \[CapitalGamma] = 1 + 1/n1

Out[9]= 2

prs[d_] := K1*d^\[CapitalGamma] (*polytropic equation*)

rhs[p_] := (1/K1)*
p^(1/\[CapitalGamma]) (*inverse polytropic equation*)(*set critical density, and other initial values*)

rho0 = 1.0 *10^(15) (*I'm guessing at this value*)

Out[12]= 1.*10^15

p0 = prs[rho0] (*should be the critical density*)

Out[13]= 5.38*10^39

In[14]:= rmin = .001

Out[14]= 0.001

In[15]:= rmax = 5

Out[15]= 5

In[16]:= rhs[p[r]*p0]

Out[16]= 1.36335*10^10 Sqrt[p[r]](*solve the equations*)

Clear

In[20]:= s =
NDSolve[{D[m[r], r] == If[p[r] > 0.0, 3*(1.3633547078730297`*^10 Sqrt[p[r]])/rhosol*r^2, 0.0],

D[p[r], r]* r*(m[r] - r) == If[p[r] > 0.0, 0.5*(1.3633547078730297`*^10 Sqrt[p[r]])*c^2/p0*m[r]*(1.0 +
p[r]*p0/(1.3633547078730297`*^10 Sqrt[p[r]]*c^2))*(1.0 + 3*p0/(rhosol*c^2)*p[r]*r^3/m[r]), 0.0],

p[rmin] == 1.0 , m[rmin] == rho0*4.0/3.0*Pi*rmin^3/msol}, {m, p}, {r, rmin, rmax},

Method →{"ExplicitRungeKutta", "DifferenceOrder" → 8} ];

NDSolve::derlen: The length of the derivative operator Derivative[1] in (m^\[Prime])[r] is not the same as the number of arguments. >>

I've been trying to do this modeling for days with no success. Any help is greatly, greatly, appreciated.