Group Abstract

Message Boards

WOLFRAM COMMUNITY

3.9K Views

1 Reply

0 Total Likes

View groups...

Follow this post

Share this post:

GROUPS:

Is it possible to impose a variable to vanish in NDSolve ?

Thibaut Arnoulx de Pirey

Posted 10 years ago

Hi everyone ! I am integrating numerically on Mathematica the TOV equations (the equations for the interior of neutron stars) in some theory of modified gravity. When the density ([Rho]) first goes to 0, that means that you have reached the surface of the star. From this point I would like to keep integrating the Einstein equation (in vacuum then) by imposing the density to remain 0. I join my code below. It tried to do so with two WhenEvent condition and by adding a discrete variable to annul the derivative of the density but it doesn't work... Thank's a lot ! Code : \[Beta] = - 4.5 ; \[CapitalGamma] = 2.34 (dimensionless) ; K = 0.0195 (dimensionless) ; \[Rho]0 = 10^(18) (kg.m^(-3)) ; c = 3 10^8 (m.s^-1) ; G = 6.67 10^(-11) (SI) ; \[Epsilon][r_] = (\[Rho][r] + K \[Rho]0 (\[Rho][r]/\[Rho]0)^\[CapitalGamma] ) c^2 ; p[r_] = K \[Rho]0 (\[Rho][r]/\[Rho]0)^\[CapitalGamma] c^2 (\[CapitalGamma] - 1) ; A[r_] = Exp[(1/2) \[Beta] \[CurlyPhi][r]^2] ; Eqn1 = (\[Mu]'[r] - (4 Pi G/c^2) r^2 (\[Epsilon][r]/c^2) (A[r])^4 - (1/2) r (r - 2 \[Mu][r]) (\[Psi][r])^2 == 0) ; Eqn2 = (\[Nu]'[r] - (8 Pi G/c^4) r^2 (A[r])^4 p[r]/(r - 2 \[Mu][r]) - r (\[Psi][r])^2 + 2 \[Mu][r]/(r (r - 2 \[Mu][r])) == 0) ; Eqn3 = (\[CurlyPhi]'[r] - \[Psi][r] == 0) ; Eqn4 = (\[Psi]'[r] - (4 Pi G/c^4) r (A[r])^4/(r - 2 \[Mu][r]) (\[Beta] \[CurlyPhi][r] (\[Epsilon][r] - 3 p[r]) + r \[Psi][r] (\[Epsilon][r] - p[r])) + 2 (r - \[Mu][r]) \[Psi[r]/(r (r - 2 \[Mu][r])) == 0) ; Eqn5 = (p'[r] + on[r] ((\[Epsilon][r] + p[r]) ((4 Pi G/c^4) r^2 (A[r])^4 p[r]/(r - 2 \[Mu][r]) + (1/2) r (\[Psi][r])^2 + \[Mu][r]/(r (r - 2 \[Mu][r])) + \[Beta]\[CurlyPhi][r] \[Psi][r])) == 0) ; sol = NDSolve[{Eqn1, Eqn2, Eqn3, Eqn4, Eqn5, \[Mu][10^(-8)] == 0, \[Nu][10^(-8)] == 0, \[Rho][10^(-8)] == \[Rho]0, \[Psi][10^(-8)] == 0, \[CurlyPhi][10^(-8)] == 0.5, on[10^(-8)] == 1, WhenEvent[\[Rho][r] < 3 10^(17), on[r] -> 0], WhenEvent[\[Rho][r] < 3 10^(17), \[Rho][r] -> 0]}, {\[Mu][r], \[Nu][r], \[Rho][r], \[Psi][r], \[CurlyPhi][r]}, {r, 10^(-8), 30000}, DiscreteVariables -> {on[r] \[Element] {0, 1}}, Method -> {"StiffnessSwitching"}]

Hi everyone ! I am integrating numerically on Mathematica the TOV equations (the equations for the interior of neutron stars) in some theory of modified gravity. When the density ([Rho]) first goes to 0, that means that you have reached the surface of the star. From this point I would like to keep integrating the Einstein equation (in vacuum then) by imposing the density to remain 0. I join my code below. It tried to do so with two WhenEvent condition and by adding a discrete variable to annul the derivative of the density but it doesn't work... Thank's a lot !

Code :

\[Beta] = - 4.5 ;

\[CapitalGamma] = 2.34 (*dimensionless*) ;

K = 0.0195 (*dimensionless*) ;

\[Rho]0 = 10^(18) (*kg.m^(-3)*) ;

c = 3 10^8 (*m.s^-1*) ;

G = 6.67 10^(-11) (*SI*) ;

\[Epsilon][r_] = (\[Rho][r] +  K \[Rho]0 (\[Rho][r]/\[Rho]0)^\[CapitalGamma] ) c^2 ;

p[r_] = K \[Rho]0 (\[Rho][r]/\[Rho]0)^\[CapitalGamma] c^2 (\[CapitalGamma] - 1) ;

A[r_] = Exp[(1/2) \[Beta] \[CurlyPhi][r]^2] ;

Eqn1 = (\[Mu]'[r] - (4 Pi G/c^2) r^2 (\[Epsilon][r]/c^2) (A[r])^4 - (1/2) r (r - 2 \[Mu][r]) (\[Psi][r])^2 == 0) ;

Eqn2 = (\[Nu]'[r] - (8 Pi G/c^4) r^2 (A[r])^4 p[r]/(r - 2 \[Mu][r]) - r (\[Psi][r])^2 + 2 \[Mu][r]/(r (r - 2 \[Mu][r])) == 0) ;

Eqn3 = (\[CurlyPhi]'[r] - \[Psi][r] == 0) ;

Eqn4 = (\[Psi]'[r] - (4 Pi G/c^4) r (A[r])^4/(r - 2 \[Mu][r]) (\[Beta] \[CurlyPhi][r] (\[Epsilon][r] - 3 p[r]) + r \[Psi][r] (\[Epsilon][r] - p[r])) + 2 (r - \[Mu][r]) \[Psi[r]/(r (r - 2 \[Mu][r])) == 0) ;

Eqn5 = (p'[r] + on[r] ((\[Epsilon][r] + p[r]) ((4 Pi G/c^4) r^2 (A[r])^4 p[r]/(r - 2 \[Mu][r]) + (1/2) r (\[Psi][r])^2 + \[Mu][r]/(r (r - 2 \[Mu][r])) + \[Beta]\[CurlyPhi][r] \[Psi][r])) == 0) ;

sol = NDSolve[{Eqn1, Eqn2, Eqn3, Eqn4, Eqn5, \[Mu][10^(-8)] == 0, \[Nu][10^(-8)] == 0, \[Rho][10^(-8)] == \[Rho]0, \[Psi][10^(-8)] == 0, \[CurlyPhi][10^(-8)] == 0.5, on[10^(-8)] == 1, WhenEvent[\[Rho][r] < 3 10^(17), on[r] -> 0], WhenEvent[\[Rho][r] < 3 10^(17), \[Rho][r] -> 0]}, {\[Mu][r], \[Nu][r], \[Rho][r], \[Psi][r], \[CurlyPhi][r]}, {r, 10^(-8), 30000}, DiscreteVariables -> {on[r] \[Element] {0, 1}}, Method -> {"StiffnessSwitching"}]

POSTED BY: Thibaut Arnoulx de Pirey